From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 24 Jun 2019 14:43:36 +0200
Subject: Initialize at openlb-1-3
---
src/functors/analytical/frameChangeF3D.h | 521 +++++++++++++++++++++++++++++++
1 file changed, 521 insertions(+)
create mode 100644 src/functors/analytical/frameChangeF3D.h
(limited to 'src/functors/analytical/frameChangeF3D.h')
diff --git a/src/functors/analytical/frameChangeF3D.h b/src/functors/analytical/frameChangeF3D.h
new file mode 100644
index 0000000..066f4ba
--- /dev/null
+++ b/src/functors/analytical/frameChangeF3D.h
@@ -0,0 +1,521 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2013, 2015 Gilles Zahnd, Mathias J. Krause
+ * Marie-Luise Maier
+ * E-mail contact: info@openlb.net
+ * The most recent release of OpenLB can be downloaded at
+ *
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public
+ * License along with this program; if not, write to the Free
+ * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ * Boston, MA 02110-1301, USA.
+*/
+
+#ifndef FRAME_CHANGE_F_3D_H
+#define FRAME_CHANGE_F_3D_H
+
+#include
+#include
+
+#include "analyticalF.h"
+
+/** \file
+ This file contains two different classes of functors, in the FIRST part
+ - for simulations in a rotating frame
+ - different functors for
+ velocity (3d, RotatingLinear3D),
+ pressure (1d, RotatingQuadratic1D) and
+ force (3d, RotatingForceField3D)
+ The functors return the displacement of a point x in a fixed amount of time.
+
+ The ones in the SECOND part are useful to set Poiseuille velocity profiles on
+ - pipes with round cross-section and
+ - pipes with square-shaped cross-section.
+*/
+
+/** To enable simulations in a rotating frame, the axis is set in the
+ * constructor with axisPoint and axisDirection. The axisPoint can be the
+ * coordinate of any point on the axis. The axisDirection has to be a normed to
+ * 1. The pulse w is in rad/s. It determines the pulse speed by its norm while
+ * the trigonometric or clockwise direction is determined by its sign: When the
+ * axisDirection is pointing "towards you", a positive pulse makes it turn in
+ * the trigonometric way. It has to be noticed that putting both axisDirection
+ * into -axisDirection and w into -w yields an exactly identical situation.
+ */
+
+namespace olb {
+
+
+template class SuperGeometry3D;
+
+// PART 1: /////////////////////////////////////////////////////////////////////
+// Functors for rotating the coordinate system (velocity, pressure, force,...)
+
+/**
+ * This functor gives a linar profile for a given point x as it computes
+ * the distance between x and the axis.
+ *
+ * The field in outcome is the velocity field of q rotating solid
+ */
+
+/// Functor with a linear profile e.g. for rotating velocity fields.
+template
+class RotatingLinear3D final : public AnalyticalF3D {
+protected:
+ std::vector axisPoint;
+ std::vector axisDirection;
+ T w;
+ T scale;
+public:
+ RotatingLinear3D(std::vector axisPoint_, std::vector axisDirection_, T w_, T scale_=1);
+ bool operator()(T output[], const T x[]) override;
+};
+
+/**
+ * This functor gives a linar profile in an annulus for a given point x between the inner and outer radius as it computes
+ * the distance between x and the inner and outer radius.
+ *
+ * The field in outcome is the velocity field of q rotating solid in an annulus
+ */
+
+/// Functor with a linear profile e.g. for rotating velocity fields.
+template
+class RotatingLinearAnnulus3D final : public AnalyticalF3D {
+protected:
+ std::vector axisPoint;
+ std::vector axisDirection;
+ T w;
+ T ri;
+ T ro;
+ T scale;
+public:
+ RotatingLinearAnnulus3D(std::vector axisPoint_, std::vector axisDirection_, T w_, T ri_, T ro_, T scale_=1);
+ bool operator()(T output[], const T x[]);
+};
+
+
+/**
+ * This functor gives a parabolic profile for a given point x as it computes
+ * the distance between x and the axis.
+ *
+ * This field is a scalar field, a vector with one component will be used
+ */
+
+/// Functor with a parabolic profile e.g. for rotating pressure fields.
+template
+class RotatingQuadratic1D final : public AnalyticalF3D {
+protected:
+ std::vector axisPoint;
+ std::vector axisDirection;
+ T w;
+ T scale;
+ T additive;
+public:
+ RotatingQuadratic1D(std::vector axisPoint_, std::vector axisDirection_,
+ T w_, T scale_=1, T additive_=0);
+ bool operator()(T output[], const T x[]) override;
+};
+
+
+// PART 2: /////////////////////////////////////////////////////////////////////
+// Functors for setting velocities on a velocity boundary of a pipe
+
+/**
+ * This functor returns a quadratic Poiseuille profile for use with a pipe with
+ * round cross-section. It uses cylinder coordinates and is valid for the
+ * entire length of the pipe.
+ *
+ * This functor gives a parabolic velocity profile for a given point x as it
+ * computes the distance between x and the axis.
+ *
+ * The axis is set in the input with axisPoint and axisDirection. The axisPoint
+ * can be the coordinate of any point where the axis passes.
+ * axisDirection has to be normed to 1.
+ * Once the axis is set in the middle of the pipe, the radius of the
+ * pipe "radius" and the velocity in the middle of the pipe "maxVelocity"
+ * determine the Poisseuille profile entierly.
+ */
+
+/// Velocity profile for round pipes and power law fluids: u(r)=u_max*(1-(r/R)^((n+1)/n)). The exponent n characterizes the fluid behavior.
+/// n<1: Pseudoplastic, n=1: Newtonian fluid, n>1: Dilatant
+template
+class CirclePowerLaw3D : public AnalyticalF3D {
+protected:
+ std::vector _center;
+ std::vector _normal;
+ T _radius;
+ T _maxVelocity;
+ T _n;
+ T _scale;
+
+public:
+ CirclePowerLaw3D(std::vector axisPoint, std::vector axisDirection, T maxVelocity, T radius, T n, T scale = T(1));
+ CirclePowerLaw3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T n, T scale = T(1));
+ CirclePowerLaw3D(SuperGeometry3D& superGeometry, int material, T maxVelocity, T n, T scale = T(1));
+
+ CirclePowerLaw3D(bool useMeanVelocity, std::vector axisPoint, std::vector axisDirection, T Velocity, T radius, T n, T scale = T(1));
+ CirclePowerLaw3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T n, T scale = T(1));
+ CirclePowerLaw3D(bool useMeanVelocity, SuperGeometry3D& superGeometry, int material, T Velocity, T n, T scale = T(1));
+
+ /// Returns centerpoint vector
+ std::vector getCenter()
+ {
+ return _center;
+ };
+ /// Returns normal vector
+ std::vector getNormal()
+ {
+ return _normal;
+ };
+ /// Returns radi
+ T getRadius()
+ {
+ return _radius;
+ };
+
+ bool operator()(T output[], const T x[]) override;
+};
+
+/// Velocity profile for round pipes and turbulent flows: u(r)=u_max*(1-r/R)^(1/n) The exponent n can be calculated by n = 1.03 * ln(Re) - 3.6
+/// n=7 is used for many flow applications
+template
+class CirclePowerLawTurbulent3D : public CirclePowerLaw3D {
+public:
+ CirclePowerLawTurbulent3D(std::vector axisPoint_, std::vector axisDirection, T maxVelocity, T radius, T n, T scale = T(1));
+ CirclePowerLawTurbulent3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T n, T scale = T(1));
+ CirclePowerLawTurbulent3D(SuperGeometry3D& superGeometry, int material, T maxVelocity, T n, T scale = T(1));
+
+ CirclePowerLawTurbulent3D(bool useMeanVelocity, std::vector axisPoint, std::vector axisDirection, T Velocity, T radius, T n, T scale = T(1));
+ CirclePowerLawTurbulent3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T n, T scale = T(1));
+ CirclePowerLawTurbulent3D(bool useMeanVelocity, SuperGeometry3D& superGeometry, int material, T Velocity, T n, T scale = T(1));
+
+ bool operator()(T output[], const T x[]) override;
+};
+
+/// Velocity profile for round pipes and a laminar flow of a Newtonian fluid: u(r)=u_max*(1-(r/R)^2)
+template
+class CirclePoiseuille3D final : public CirclePowerLaw3D {
+
+public:
+ CirclePoiseuille3D(std::vector axisPoint, std::vector axisDirection, T maxVelocity, T radius, T scale = T(1));
+ CirclePoiseuille3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T scale = T(1));
+ CirclePoiseuille3D(SuperGeometry3D& superGeometry, int material, T maxVelocity, T scale = T(1));
+
+ CirclePoiseuille3D(bool useMeanVelocity, std::vector axisPoint, std::vector axisDirection, T Velocity, T radius, T scale = T(1));
+ CirclePoiseuille3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T scale = T(1));
+ CirclePoiseuille3D(bool useMeanVelocity, SuperGeometry3D& superGeometry, int material, T Velocity, T scale = T(1));
+};
+
+/// Strain rate for round pipes and laminar flow of a Newtonian fluid
+template < typename T,typename DESCRIPTOR>
+class CirclePoiseuilleStrainRate3D : public AnalyticalF3D {
+protected:
+ T lengthY;
+ T lengthZ;
+ T maxVelocity;
+
+public:
+ CirclePoiseuilleStrainRate3D(UnitConverter const& converter, T ly);
+ bool operator()(T output[], const T input[]);
+};
+
+/**
+ This functor returns a Poiseuille profile for use with a pipe with square shaped cross-section.
+*/
+/// velocity profile for pipes with rectangular cross-section.
+template
+class RectanglePoiseuille3D final : public AnalyticalF3D {
+protected:
+ OstreamManager clout;
+ std::vector x0;
+ std::vector x1;
+ std::vector x2;
+ std::vector maxVelocity;
+
+public:
+ RectanglePoiseuille3D(std::vector& x0_, std::vector& x1_, std::vector& x2_, std::vector& maxVelocity_);
+ /// constructor from material numbers
+ /** offsetX,Y,Z is a positive number denoting the distance from border
+ * voxels of material_ to the zerovelocity boundary */
+ RectanglePoiseuille3D(SuperGeometry3D& superGeometry_, int material_, std::vector& maxVelocity_, T offsetX, T offsetY, T offsetZ);
+ bool operator()(T output[], const T x[]) override;
+};
+
+
+/**
+ * This functor returns a quadratic Poiseuille profile for use with a pipe with
+ * elliptic cross-section.
+ *
+ * Ellpise is orthogonal to z-direction.
+ * a is in x-direction
+ * b is in y-direction
+ *
+ */
+/// Velocity profile for round pipes.
+template
+class EllipticPoiseuille3D final : public AnalyticalF3D {
+protected:
+ OstreamManager clout;
+ std::vector _center;
+ T _a2, _b2;
+ T _maxVel;
+
+public:
+ EllipticPoiseuille3D(std::vector center, T a, T b, T maxVel);
+ bool operator()(T output[], const T x[]) override;
+};
+
+
+/// Analytical solution of porous media channel flow with low Reynolds number
+/// See Spaid and Phelan (doi:10.1063/1.869392)
+template
+class AnalyticalPorousVelocity3D : public AnalyticalF3D {
+protected:
+ std::vector center;
+ std::vector normal;
+ T K, mu, gradP, radius;
+ T eps;
+public:
+ AnalyticalPorousVelocity3D(SuperGeometry3D& superGeometry, int material, T K_, T mu_, T gradP_, T radius_, T eps_=T(1));
+ T getPeakVelocity();
+ bool operator()(T output[], const T input[]) override;
+};
+
+
+////////////////Coordinate Transformation//////////////
+
+
+/// This class calculates the angle alpha between vector _referenceVector and
+/// any other vector x.
+/// Vector x has to be turned by alpha in mathematical positive sense depending
+/// to _orientation to lie over vector _referenceVector
+template
+class AngleBetweenVectors3D final : public AnalyticalF3D {
+protected:
+ /// between _referenceVector and vector x, angle is calculated
+ std::vector _referenceVector;
+ /// direction around which x has to be turned with angle is in math pos sense
+ /// to lie over _referenceVector
+ std::vector _orientation;
+public:
+ /// constructor defines _referenceVector and _orientation
+ AngleBetweenVectors3D(std::vector referenceVector,
+ std::vector orientation);
+ /// operator writes angle between x and _referenceVector inro output field
+ bool operator()(T output[], const S x[]) override;
+};
+
+
+/// This class saves coordinates of the resulting point after rotation in the
+/// output field. The resulting point is achieved after rotation of x with
+/// angle _alpha in math positive sense relative to _rotAxisDirection with a
+/// given _origin.
+template
+class RotationRoundAxis3D final : public AnalyticalF3D {
+protected:
+ /// origin, around which x is turned
+ std::vector _origin;
+ /// direction, around which x is turned
+ std::vector _rotAxisDirection;
+ /// angle, by which vector x is rotated around _rotAxisDirection
+ T _alpha;
+public:
+ /// constructor defines _rotAxisDirection, _alpha and _origin
+ RotationRoundAxis3D(std::vector origin, std::vector rotAxisDirection,
+ T alpha);
+ /// operator writes coordinates of the resulting point after rotation
+ /// of x by angle _alpha in math positive sense to _rotAxisDirection
+ /// with _origin into output field
+ bool operator()(T output[], const S x[]) override;
+};
+
+
+////////// Cylinder & Cartesian ///////////////
+
+
+/// This class converts cylindrical of point x (x[0] = radius, x[1] = phi,
+/// x[2] = z) to Cartesian coordinates (wrote into output field).
+/// Initial situation for the Cartesian coordinate system is that angle phi
+/// lies in the x-y-plane and turns round the _cylinderOrigin, the z-axis
+/// direction equals the axis direction of the cylinder.
+template
+class CylinderToCartesian3D final : public AnalyticalF3D {
+protected:
+ std::vector _cylinderOrigin;
+public:
+ /// constructor defines _cylinderOrigin
+ CylinderToCartesian3D(std::vector cylinderOrigin);
+ /// operator writes Cartesian coordinates of cylinder coordinates
+ /// x[0] = radius >= 0, x[1] = phi in [0, 2*Pi), x[2] = z into output field
+ bool operator()(T output[], const S x[]) override;
+};
+
+
+/// This class converts Cartesian coordinates of point x to cylindrical
+/// coordinates wrote into output field
+/// (output[0] = radius, output[1] = phi, output[2] = z).
+/// Initial situation for the cylindrical coordinate system is that angle phi
+/// lies in plane perpendicular to the _axisDirection and turns around the
+/// _cartesianOrigin. The radius is the distance of point x to the
+/// _axisDirection and z the distance to _cartesianOrigin along _axisDirection.
+template
+class CartesianToCylinder3D final : public AnalyticalF3D {
+protected:
+ /// origin of the Cartesian coordinate system
+ std::vector _cartesianOrigin;
+ /// direction of the axis along which the cylindrical coordinates are calculated
+ std::vector _axisDirection;
+ /// direction to know orientation for math positive to obtain angle phi
+ /// of Cartesian point x
+ std::vector _orientation;
+public:
+ CartesianToCylinder3D(std::vector cartesianOrigin, T& angle,
+ std::vector orientation = {T(1),T(),T()});
+ CartesianToCylinder3D(std::vector cartesianOrigin,
+ std::vector axisDirection,
+ std::vector orientation = {T(1),T(),T()});
+ CartesianToCylinder3D(T cartesianOriginX, T cartesianOriginY,
+ T cartesianOriginZ,
+ T axisDirectionX, T axisDirectionY,
+ T axisDirectionZ,
+ T orientationX = T(1), T orientationY = T(),
+ T orientationZ = T());
+ /// operator writes cylindrical coordinates of Cartesian point x into output
+ /// field,
+ /// output[0] = radius ( >= 0), output[1] = phi (in [0, 2Pi)), output[2] = z
+ bool operator()(T output[], const S x[]) override;
+ /// Read only access to _axisDirection
+ std::vector const& getAxisDirection() const;
+ /// Read and write access to _axisDirection
+ std::vector& getAxisDirection();
+ /// Read only access to _catesianOrigin
+ std::vector const& getCartesianOrigin() const;
+ /// Read and write access to _catesianOrigin
+ std::vector& getCartesianOrigin();
+ /// Read and write access to _axisDirection using angle
+ void setAngle(T angle);
+};
+
+
+////////// Spherical & Cartesian ///////////////
+
+
+/// This class converts spherical coordinates of point x
+/// (x[0] = radius, x[1] = phi, x[2] = theta) to Cartesian coordinates wrote
+/// into output field.
+/// Initial situation for the Cartesian coordinate system is that angle phi
+/// (phi in [0, 2Pi)) lies in x-y-plane and turns around the Cartesian origin
+/// (0,0,0). Angle theta (theta in [0, Pi]) lies in y-z-plane. z axis equals
+/// the _orientation of the spherical coordinate system.
+/// The radius is distance of point x to the Cartesian _origin
+template
+class SphericalToCartesian3D final : public AnalyticalF3D {
+ //protected:
+public:
+ SphericalToCartesian3D();
+ /// operator writes spherical coordinates of Cartesian coordinates of x
+ /// (x[0] = radius, x[1] = phi, x[2] = theta) into output field.
+ /// phi is in x-y-plane, theta in x-z-plane, z axis as orientation
+ bool operator()(T output[], const S x[]) override;
+};
+
+
+/// This class converts Cartesian coordinates of point x to spherical
+/// coordinates wrote into output field
+/// (output[0] = radius, output[1] = phi, output[2] = theta).
+/// Initial situation for the spherical coordinate system is that angle phi
+/// lies in x-y-plane and theta in x-z-plane.
+/// The z axis is the orientation for the mathematical positive sense of phi.
+template
+class CartesianToSpherical3D final : public AnalyticalF3D {
+protected:
+ /// origin of the Cartesian coordinate system
+ std::vector _cartesianOrigin;
+ /// direction of the axis along which the spherical coordinates are calculated
+ std::vector _axisDirection;
+public:
+ CartesianToSpherical3D(std::vector cartesianOrigin,
+ std::vector axisDirection);//, std::vector orientation);
+ /// operator writes shperical coordinates of Cartesian point x into output
+ /// field,
+ /// output[0] = radius ( >= 0), output[1] = phi (in [0, 2Pi)),
+ /// output[2] = theta (in [0, Pi])
+ bool operator()(T output[], const S x[]) override;
+};
+
+
+////////// Magnetic Field ///////////////
+
+
+/// Magnetic field that creates magnetization in wire has to be orthogonal to
+/// the wire. To calculate the magnetic force on particles around a cylinder
+/// (J. Lindner et al. / Computers and Chemical Engineering 54 (2013) 111-121)
+/// Fm = mu0*4/3.*PI*radParticle^3*_Mp*radWire^2/r^3 *
+/// [radWire^2/r^2+cos(2*theta), sin(2*theta), 0]
+template
+class MagneticForceFromCylinder3D final : public AnalyticalF3D {
+protected:
+ CartesianToCylinder3D& _car2cyl;
+ /// length of the wire, from origin to _car2cyl.axisDirection
+ T _length;
+ /// wire radius
+ T _radWire;
+ /// maximal distance from wire cutoff/threshold
+ T _cutoff;
+ /// saturation magnetization wire, linear scaling factor
+ T _Mw;
+ /// magnetization particle, linear scaling factor
+ T _Mp;
+ /// particle radius, cubic scaling factor
+ T _radParticle;
+ /// permeability constant, linear scaling factor
+ T _mu0;
+ /// factor = mu0*4/3.*PI*radParticle^3*_Mp*radWire^2/r^3
+ T _factor;
+public:
+ MagneticForceFromCylinder3D(CartesianToCylinder3D& car2cyl, T length,
+ T radWire, T cutoff, T Mw, T Mp = T(1),
+ T radParticle = T(1), T mu0 = 4*3.14159265e-7);
+ /// operator writes the magnetic force in a point x round a cylindrical wire into output field
+ bool operator()(T output[], const S x[]) override;
+};
+
+template
+class MagneticFieldFromCylinder3D final : public AnalyticalF3D {
+protected:
+ CartesianToCylinder3D& _car2cyl;
+ /// length of the wire, from origin to _car2cyl.axisDirection
+ T _length;
+ /// wire radius
+ T _radWire;
+ /// maximal distance from wire cutoff/threshold
+ T _cutoff;
+ /// saturation magnetization wire, linear scaling factor
+ T _Mw;
+ /// factor = mu0*4/3.*PI*radParticle^3*_Mp*radWire^2/r^3
+ T _factor;
+public:
+ MagneticFieldFromCylinder3D(CartesianToCylinder3D& car2cyl,
+ T length,
+ T radWire,
+ T cutoff,
+ T Mw
+ );
+ /// operator writes the magnetic force in a point x round a cylindrical wire into output field
+ bool operator()(T output[], const S x[]);
+};
+
+} // end namespace olb
+
+#endif
--
cgit v1.2.3