KratosMultiphysics
KRATOS Multiphysics (Kratos) is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface.
two_fluid_navier_stokes_data.h
Go to the documentation of this file.
1 // | / |
2 // ' / __| _` | __| _ \ __|
3 // . \ | ( | | ( |\__ `
4 // _|\_\_| \__,_|\__|\___/ ____/
5 // Multi-Physics
6 //
7 // License: BSD License
8 // Kratos default license: kratos/license.txt
9 //
10 // Main authors: Daniel Diez, Jordi Cotela
11 //
12 
13 
14 #if !defined(KRATOS_TWO_FLUID_NAVIER_STOKES_DATA_H)
15 #define KRATOS_TWO_FLUID_NAVIER_STOKES_DATA_H
16 
18 
23 
24 namespace Kratos {
25 
28 
31 
32 template< size_t TDim, size_t TNumNodes >
33 class TwoFluidNavierStokesData : public FluidElementData<TDim,TNumNodes, true>
34 {
35 public:
36 
39 
47 
51 
57 
62 
63 double Density;
65 double DeltaTime; // Time increment
66 double DynamicTau; // Dynamic tau considered in ASGS stabilization coefficients
70 double DarcyTerm;
71 double VolumeError;
72 double bdf0;
73 double bdf1;
74 double bdf2;
75 
76 // Auxiliary containers for the symbolically-generated matrices
77 BoundedMatrix<double,TNumNodes*(TDim+1),TNumNodes*(TDim+1)> lhs;
78 array_1d<double,TNumNodes*(TDim+1)> rhs;
79 BoundedMatrix<double, TNumNodes*(TDim + 1), TNumNodes> V;
80 BoundedMatrix<double, TNumNodes, TNumNodes*(TDim + 1)> H;
83 
84 double ElementSize;
85 
90 
93 
96 
99 
102 unsigned int NumberOfDivisions;
103 
104 
108 
109 void Initialize(const Element& rElement, const ProcessInfo& rProcessInfo) override
110 {
111  // Base class Initialize manages constitutive law parameters
113 
114  const Geometry< Node >& r_geometry = rElement.GetGeometry();
115  const Properties& r_properties = rElement.GetProperties();
116  this->FillFromHistoricalNodalData(Velocity,VELOCITY,r_geometry);
117  this->FillFromHistoricalNodalData(Velocity_OldStep1,VELOCITY,r_geometry,1);
118  this->FillFromHistoricalNodalData(Velocity_OldStep2,VELOCITY,r_geometry,2);
119  this->FillFromHistoricalNodalData(Distance, DISTANCE, r_geometry);
120  this->FillFromHistoricalNodalData(MeshVelocity,MESH_VELOCITY,r_geometry);
121  this->FillFromHistoricalNodalData(BodyForce,BODY_FORCE,r_geometry);
122  this->FillFromHistoricalNodalData(Pressure,PRESSURE,r_geometry);
123  this->FillFromHistoricalNodalData(NodalDensity, DENSITY, r_geometry);
124  this->FillFromHistoricalNodalData(NodalDynamicViscosity, DYNAMIC_VISCOSITY, r_geometry);
125  this->FillFromProperties(SmagorinskyConstant, C_SMAGORINSKY, r_properties);
126  this->FillFromProperties(LinearDarcyCoefficient, LIN_DARCY_COEF, r_properties);
127  this->FillFromProperties(NonLinearDarcyCoefficient, NONLIN_DARCY_COEF, r_properties);
128  this->FillFromProcessInfo(DeltaTime,DELTA_TIME,rProcessInfo);
129  this->FillFromProcessInfo(DynamicTau,DYNAMIC_TAU,rProcessInfo);
130  this->FillFromProcessInfo(VolumeError,VOLUME_ERROR,rProcessInfo);
131 
132  const Vector& BDFVector = rProcessInfo[BDF_COEFFICIENTS];
133  bdf0 = BDFVector[0];
134  bdf1 = BDFVector[1];
135  bdf2 = BDFVector[2];
136 
137  noalias(lhs) = ZeroMatrix(TNumNodes*(TDim+1),TNumNodes*(TDim+1));
138  noalias(rhs) = ZeroVector(TNumNodes*(TDim+1));
139  noalias(V) = ZeroMatrix(TNumNodes*(TDim + 1), TNumNodes);
140  noalias(H) = ZeroMatrix(TNumNodes, TNumNodes*(TDim + 1));
141  noalias(Kee) = ZeroMatrix(TNumNodes, TNumNodes);
142  noalias(rhs_ee) = ZeroVector(TNumNodes);
143 
144  NumPositiveNodes = 0;
145  NumNegativeNodes = 0;
146 
147  for (unsigned int i = 0; i < TNumNodes; i++){
148  if(Distance[i] > 0)
150  else
152  }
153 }
154 
156  unsigned int IntegrationPointIndex,
157  double NewWeight,
158  const MatrixRowType& rN,
159  const BoundedMatrix<double, TNumNodes, TDim>& rDN_DX) override
160 {
164 }
165 
167  unsigned int IntegrationPointIndex,
168  double NewWeight,
169  const MatrixRowType& rN,
171  const MatrixRowType& rNenr,
173 {
176  noalias(this->Nenr) = rNenr;
177  noalias(this->DN_DXenr) = rDN_DXenr;
179 }
180 
181 static int Check(const Element& rElement, const ProcessInfo& rProcessInfo)
182 {
183  const Geometry< Node >& r_geometry = rElement.GetGeometry();
184 
185  for (unsigned int i = 0; i < TNumNodes; i++)
186  {
187  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(VELOCITY,r_geometry[i]);
188  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(DISTANCE, r_geometry[i]);
189  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(MESH_VELOCITY,r_geometry[i]);
190  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(BODY_FORCE,r_geometry[i]);
191  KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(PRESSURE,r_geometry[i]);
192  }
193 
194  return 0;
195 }
196 
197 bool IsCut() {
198  return (NumPositiveNodes > 0) && (NumNegativeNodes > 0);
199 }
200 
201 bool IsAir() {
202  return (NumPositiveNodes == TNumNodes);
203 }
204 
206  const unsigned int strain_size = 3 * (TDim - 1);
207 
208  if(this->C.size1() != strain_size)
209  this->C.resize(strain_size,strain_size,false);
210  if(this->ShearStress.size() != strain_size)
211  this->ShearStress.resize(strain_size,false);
212 
213  ComputeStrain();
214 
216 
217  const double mu = this->EffectiveViscosity;
218  const double c1 = 2.0*mu;
219  const double c2 = mu;
220  this->C.clear();
222  Vector& stress = this->ShearStress;
223  Vector& strain = this->StrainRate;
224 
226  this->C = c_mat;
227 
228  if constexpr (TDim == 2)
229  {
230  const double trace = strain[0] + strain[1];
231  const double volumetric_part = trace/2.0; // Note: this should be small for an incompressible fluid (it is basically the incompressibility error)
232 
233  stress[0] = c1 * (strain[0] - volumetric_part);
234  stress[1] = c1 * (strain[1] - volumetric_part);
235  stress[2] = c2 * strain[2];
236  }
237 
238  else if constexpr (TDim == 3)
239  {
240  const double trace = strain[0] + strain[1] + strain[2];
241  const double volumetric_part = trace/3.0; // Note: this should be small for an incompressible fluid (it is basically the incompressibility error)
242 
243  stress[0] = c1*(strain[0] - volumetric_part);
244  stress[1] = c1*(strain[1] - volumetric_part);
245  stress[2] = c1*(strain[2] - volumetric_part);
246  stress[3] = c2*strain[3];
247  stress[4] = c2*strain[4];
248  stress[5] = c2*strain[5];
249  }
250 }
251 
253 {
256 
257  // Compute strain (B*v)
258  // 3D strain computation
259  if constexpr (TDim == 3)
260  {
261  this->StrainRate[0] = DN(0,0)*v(0,0) + DN(1,0)*v(1,0) + DN(2,0)*v(2,0) + DN(3,0)*v(3,0);
262  this->StrainRate[1] = DN(0,1)*v(0,1) + DN(1,1)*v(1,1) + DN(2,1)*v(2,1) + DN(3,1)*v(3,1);
263  this->StrainRate[2] = DN(0,2)*v(0,2) + DN(1,2)*v(1,2) + DN(2,2)*v(2,2) + DN(3,2)*v(3,2);
264  this->StrainRate[3] = DN(0,0)*v(0,1) + DN(0,1)*v(0,0) + DN(1,0)*v(1,1) + DN(1,1)*v(1,0) + DN(2,0)*v(2,1) + DN(2,1)*v(2,0) + DN(3,0)*v(3,1) + DN(3,1)*v(3,0);
265  this->StrainRate[4] = DN(0,1)*v(0,2) + DN(0,2)*v(0,1) + DN(1,1)*v(1,2) + DN(1,2)*v(1,1) + DN(2,1)*v(2,2) + DN(2,2)*v(2,1) + DN(3,1)*v(3,2) + DN(3,2)*v(3,1);
266  this->StrainRate[5] = DN(0,0)*v(0,2) + DN(0,2)*v(0,0) + DN(1,0)*v(1,2) + DN(1,2)*v(1,0) + DN(2,0)*v(2,2) + DN(2,2)*v(2,0) + DN(3,0)*v(3,2) + DN(3,2)*v(3,0);
267  }
268  // 2D strain computation
269  else if constexpr (TDim == 2)
270  {
271  this->StrainRate[0] = DN(0,0)*v(0,0) + DN(1,0)*v(1,0) + DN(2,0)*v(2,0);
272  this->StrainRate[1] = DN(0,1)*v(0,1) + DN(1,1)*v(1,1) + DN(2,1)*v(2,1);
273  this->StrainRate[2] = DN(0,1)*v(0,0) + DN(0,0)*v(0,1) + DN(1,1)*v(1,0) + DN(1,0)*v(1,1) + DN(2,1)*v(2,0) + DN(2,0)*v(2,1);
274  }
275 }
276 
278 {
279  double strain_rate_norm;
280  Vector& S = this->StrainRate;
281  if constexpr (TDim == 3)
282  {
283  strain_rate_norm = std::sqrt(2.*S[0] * S[0] + 2.*S[1] * S[1] + 2.*S[2] * S[2] +
284  S[3] * S[3] + S[4] * S[4] + S[5] * S[5]);
285  }
286 
287  else if constexpr (TDim == 2)
288  {
289  strain_rate_norm = std::sqrt(2.*S[0] * S[0] + 2.*S[1] * S[1] + S[2] * S[2]);
290  }
291  return strain_rate_norm;
292 }
293 
295 {
296  double dist = 0.0;
297  for (unsigned int i = 0; i < TNumNodes; i++)
298  dist += this->N[i] * Distance[i];
299 
300  int navg = 0;
301  double density = 0.0;
302  for (unsigned int i = 0; i < TNumNodes; i++)
303  {
304  if (dist * Distance[i] > 0.0)
305  {
306  navg += 1;
307  density += NodalDensity[i];
308  }
309  }
310 
311  Density = density / navg;
312 }
313 
315 {
316  double dist = 0.0;
317  for (unsigned int i = 0; i < TNumNodes; i++)
318  dist += this->N[i] * Distance[i];
319 
320  int navg = 0;
321  double dynamic_viscosity = 0.0;
322  for (unsigned int i = 0; i < TNumNodes; i++)
323  {
324  if (dist * Distance[i] > 0.0)
325  {
326  navg += 1;
327  dynamic_viscosity += NodalDynamicViscosity[i];
328  }
329  }
330  DynamicViscosity = dynamic_viscosity / navg;
331 
332  if (SmagorinskyConstant > 0.0)
333  {
334  const double strain_rate_norm = ComputeStrainNorm();
335 
336  double length_scale = SmagorinskyConstant*ElementSize;
337  length_scale *= length_scale; // square
338  this->EffectiveViscosity = DynamicViscosity + 2.0*length_scale*strain_rate_norm;
339  }
341 }
342 
344 {
345  array_1d<double, 3> convective_velocity(3, 0.0);
346  for (size_t i = 0; i < TNumNodes; i++) {
347  for (size_t j = 0; j < TDim; j++) {
348  convective_velocity[j] += this->N[i] * (Velocity(i, j) - MeshVelocity(i, j));
349  }
350  }
351  const double convective_velocity_norm = MathUtils<double>::Norm(convective_velocity);
352  DarcyTerm = this->EffectiveViscosity * LinearDarcyCoefficient + Density * NonLinearDarcyCoefficient * convective_velocity_norm;
353 }
355 
356 
357 
358 
359 
360 };
361 
363 
365 
366 }
367 
368 #endif
Base class for all Elements.
Definition: element.h:60
PropertiesType & GetProperties()
Definition: element.h:1024
static double GradientsElementSize(const BoundedMatrix< double, 3, 2 > &rDN_DX)
Element size based on the shape functions gradients. Triangle element version.
Definition: element_size_calculator.cpp:1456
Base class for data containers used within FluidElement and derived types.
Definition: fluid_element_data.h:37
void FillFromProperties(double &rData, const Variable< double > &rVariable, const Properties &rProperties)
Definition: fluid_element_data.cpp:192
void FillFromHistoricalNodalData(NodalScalarData &rData, const Variable< double > &rVariable, const Geometry< Node > &rGeometry)
Definition: fluid_element_data.cpp:65
virtual void Initialize(const Element &rElement, const ProcessInfo &rProcessInfo)
Definition: fluid_element_data.cpp:18
unsigned int IntegrationPointIndex
Definition: fluid_element_data.h:100
ShapeFunctionsType N
Definition: fluid_element_data.h:104
Vector StrainRate
Strain rate (symmetric gradient of velocity) vector in Voigt notation.
Definition: fluid_element_data.h:110
Vector ShearStress
Shear stress vector in Voigt notation.
Definition: fluid_element_data.h:114
virtual void UpdateGeometryValues(unsigned int IntegrationPointIndex, double NewWeight, const MatrixRowType &rN, const ShapeDerivativesType &rDN_DX)
Definition: fluid_element_data.cpp:52
double EffectiveViscosity
Effective viscosity (in dynamic units) produced by the constitutive law.
Definition: fluid_element_data.h:124
ShapeDerivativesType DN_DX
Definition: fluid_element_data.h:106
void FillFromProcessInfo(double &rData, const Variable< double > &rVariable, const ProcessInfo &rProcessInfo)
Definition: fluid_element_data.cpp:153
Matrix C
Constitutive tensor C (expressed as a Matrix).
Definition: fluid_element_data.h:118
static void GetNewtonianConstitutiveMatrix(const double DynamicViscosity, BoundedMatrix< double, VoigtVector2DSize, VoigtVector2DSize > &rConstitutiveMatrix)
Definition: fluid_element_utilities.cpp:58
GeometryType & GetGeometry()
Returns the reference of the geometry.
Definition: geometrical_object.h:158
Geometry base class.
Definition: geometry.h:71
Definition: amatrix_interface.h:41
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
void clear()
Definition: amatrix_interface.h:284
static double Norm(const Vector &a)
Calculates the norm of vector "a".
Definition: math_utils.h:703
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Properties encapsulates data shared by different Elements or Conditions. It can store any type of dat...
Definition: properties.h:69
Definition: two_fluid_navier_stokes_data.h:34
bool IsCut()
Definition: two_fluid_navier_stokes_data.h:197
double DarcyTerm
Definition: two_fluid_navier_stokes_data.h:70
size_t NumNegativeNodes
Definition: two_fluid_navier_stokes_data.h:101
BoundedMatrix< double, TNumNodes, TNumNodes *(TDim+1)> H
Definition: two_fluid_navier_stokes_data.h:80
static int Check(const Element &rElement, const ProcessInfo &rProcessInfo)
Definition: two_fluid_navier_stokes_data.h:181
GeometryType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType
Definition: two_fluid_navier_stokes_data.h:46
typename FluidElementData< TDim, TNumNodes, true >::MatrixRowType MatrixRowType
Definition: two_fluid_navier_stokes_data.h:44
ShapeFunctionsGradientsType DN_DX_pos_side
Definition: two_fluid_navier_stokes_data.h:88
array_1d< double, TNumNodes > rhs_ee
Definition: two_fluid_navier_stokes_data.h:82
double SmagorinskyConstant
Definition: two_fluid_navier_stokes_data.h:67
bool IsAir()
Definition: two_fluid_navier_stokes_data.h:201
double Density
Definition: two_fluid_navier_stokes_data.h:63
unsigned int NumberOfDivisions
Definition: two_fluid_navier_stokes_data.h:102
typename FluidElementData< TDim, TNumNodes, true >::NodalVectorData NodalVectorData
Definition: two_fluid_navier_stokes_data.h:41
Vector w_gauss_neg_side
Definition: two_fluid_navier_stokes_data.h:95
void UpdateGeometryValues(unsigned int IntegrationPointIndex, double NewWeight, const MatrixRowType &rN, const BoundedMatrix< double, TNumNodes, TDim > &rDN_DX, const MatrixRowType &rNenr, const BoundedMatrix< double, TNumNodes, TDim > &rDN_DXenr)
Definition: two_fluid_navier_stokes_data.h:166
ShapeFunctionsGradientsType DN_DX_neg_side
Definition: two_fluid_navier_stokes_data.h:89
void CalculateAirMaterialResponse()
Definition: two_fluid_navier_stokes_data.h:205
void CalculateDensityAtGaussPoint()
Definition: two_fluid_navier_stokes_data.h:294
void ComputeStrain()
Definition: two_fluid_navier_stokes_data.h:252
double ElementSize
Definition: two_fluid_navier_stokes_data.h:84
BoundedMatrix< double, TNumNodes *(TDim+1), TNumNodes > V
Definition: two_fluid_navier_stokes_data.h:79
NodalScalarData Distance
Definition: two_fluid_navier_stokes_data.h:59
Matrix N_pos_side
Definition: two_fluid_navier_stokes_data.h:86
Vector w_gauss_pos_side
Definition: two_fluid_navier_stokes_data.h:94
Matrix N_neg_side
Definition: two_fluid_navier_stokes_data.h:87
double bdf0
Definition: two_fluid_navier_stokes_data.h:72
typename FluidElementData< TDim, TNumNodes, true >::ShapeFunctionsType ShapeFunctionsType
Definition: two_fluid_navier_stokes_data.h:42
NodalScalarData Pressure
Definition: two_fluid_navier_stokes_data.h:58
Geometry< Node > GeometryType
Definition: two_fluid_navier_stokes_data.h:45
BoundedMatrix< double, TNumNodes, TNumNodes > Enr_Neg_Interp
Definition: two_fluid_navier_stokes_data.h:92
NodalScalarData NodalDensity
Definition: two_fluid_navier_stokes_data.h:60
BoundedMatrix< double, TNumNodes *(TDim+1), TNumNodes *(TDim+1)> lhs
Definition: two_fluid_navier_stokes_data.h:77
BoundedMatrix< double, TNumNodes, TNumNodes > Enr_Pos_Interp
Definition: two_fluid_navier_stokes_data.h:91
NodalVectorData MeshVelocity
Definition: two_fluid_navier_stokes_data.h:55
double ComputeStrainNorm()
Definition: two_fluid_navier_stokes_data.h:277
void CalculateEffectiveViscosityAtGaussPoint()
Definition: two_fluid_navier_stokes_data.h:314
double LinearDarcyCoefficient
Definition: two_fluid_navier_stokes_data.h:68
NodalVectorData Velocity_OldStep1
Definition: two_fluid_navier_stokes_data.h:53
void UpdateGeometryValues(unsigned int IntegrationPointIndex, double NewWeight, const MatrixRowType &rN, const BoundedMatrix< double, TNumNodes, TDim > &rDN_DX) override
Definition: two_fluid_navier_stokes_data.h:155
array_1d< double, TNumNodes *(TDim+1)> rhs
Definition: two_fluid_navier_stokes_data.h:78
double DynamicTau
Definition: two_fluid_navier_stokes_data.h:66
BoundedMatrix< double, TNumNodes, TNumNodes > Kee
Definition: two_fluid_navier_stokes_data.h:81
NodalScalarData NodalDynamicViscosity
Definition: two_fluid_navier_stokes_data.h:61
size_t NumPositiveNodes
Definition: two_fluid_navier_stokes_data.h:100
double bdf1
Definition: two_fluid_navier_stokes_data.h:73
NodalVectorData BodyForce
Definition: two_fluid_navier_stokes_data.h:56
NodalVectorData Velocity
Definition: two_fluid_navier_stokes_data.h:52
void Initialize(const Element &rElement, const ProcessInfo &rProcessInfo) override
Definition: two_fluid_navier_stokes_data.h:109
ShapeDerivativesType DN_DXenr
Definition: two_fluid_navier_stokes_data.h:98
NodalVectorData Velocity_OldStep2
Definition: two_fluid_navier_stokes_data.h:54
void ComputeDarcyTerm()
Definition: two_fluid_navier_stokes_data.h:343
typename FluidElementData< TDim, TNumNodes, true >::ShapeDerivativesType ShapeDerivativesType
Definition: two_fluid_navier_stokes_data.h:43
typename FluidElementData< TDim, TNumNodes, true >::NodalScalarData NodalScalarData
Definition: two_fluid_navier_stokes_data.h:40
double NonLinearDarcyCoefficient
Definition: two_fluid_navier_stokes_data.h:69
double VolumeError
Definition: two_fluid_navier_stokes_data.h:71
double DynamicViscosity
Definition: two_fluid_navier_stokes_data.h:64
double bdf2
Definition: two_fluid_navier_stokes_data.h:74
double DeltaTime
Definition: two_fluid_navier_stokes_data.h:65
ShapeFunctionsType Nenr
Definition: two_fluid_navier_stokes_data.h:97
Short class definition.
Definition: array_1d.h:61
#define KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(TheVariable, TheNode)
Definition: checks.h:171
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
TABLE_NUMBER_ANGULAR_VELOCITY TABLE_NUMBER_MOMENT I33 BEAM_INERTIA_ROT_UNIT_LENGHT_Y KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, BEAM_INERTIA_ROT_UNIT_LENGHT_Z) typedef std double
Definition: DEM_application_variables.h:182
float dist
Definition: edgebased_PureConvection.py:89
float density
Definition: face_heat.py:56
v
Definition: generate_convection_diffusion_explicit_element.py:114
DN
Definition: generate_convection_diffusion_explicit_element.py:98
stress
Stress vector definition.
Definition: generate_droplet_dynamics.py:82
mu
Definition: generate_frictional_mortar_condition.py:127
int strain_size
Definition: generate_hyper_elastic_simo_taylor_neo_hookean.py:16
strain
HERE WE MUST SUBSTITUTE EXPRESSIONS.
Definition: generate_stokes_twofluid_element.py:104
S
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:39
int j
Definition: quadrature.py:648
integer i
Definition: TensorModule.f:17