Mechanic PDE¶
Governing Equations¶
The governing equation can easily be derived from the conservation equation for the momentum density (\rho \dot{\mathbf{u}}). Assuming that the force density \mathbf{f} can be represented as \mathbf{f} = \nabla \cdot \mathbf{\sigma} + \mathbf{g}, we obtain \begin{equation} \frac{\partial^2 \rho \mathbf{u}}{\partial t^2} - \nabla \cdot \mathbf{\sigma} = \mathbf{g}, \end{equation} where \rho denotes the density, \mathbf{u} the displacement, \mathbf{\sigma} the stress tensor (rank 2) and \mathbf{g} the force density. It has to be noted that \mathbf{f} is a body force density which is exerted from the outside. We assume linear, elastic material behaviour which may be anisotropic and express it via Hooks Law which is given by \begin{equation} \mathbf{\sigma} = \mathbf{C} : \mathbf{s}, \end{equation} where \mathbf{C} denotes the stiffness tensor (rank 4) and \mathbf{s} the strain tensor (rank 2). Furthermore, we assume small strains in order to use the linearized strain displacement relationship given as \begin{equation} \mathbf{s} = \frac{1}{2} ( \nabla\mathbf{u} + (\nabla\mathbf{u})^{\textrm{T}}). \end{equation} Multiplying (1) with a the test function \mathbf{u}^\prime and inserting (2) and (3) finally leads to the weak form of the PDE given by \begin{equation} \int_{\Omega} \mathbf{u}^\prime \cdot \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} \mathrm{d}\Omega + \frac{1}{2} \int_{\Omega} (\nabla \mathbf{u}^\prime) : \mathbf{C} : \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^{\textrm{T}} \right) \mathrm{d}\Omega - \int_{\Gamma} (\mathbf{u}^\prime \cdot \mathbf{\sigma}) \cdot \mathbf{n} \mathrm{d}\Gamma = \int_{\Omega} \mathbf{u}^\prime \cdot \mathbf{g} \, \mathrm{d}\Omega. \end{equation}
Boundary conditions¶
The boundary conditions can be given as
with \mathbf{u}_{\textrm e} as the Dirichlet (essential) BC-value, \Gamma_{\textrm e} the essential boundary, \mathbf{n} the normal vector on the Neumann (natural) boundary \Gamma_{\textrm n} and \mathbf{\sigma}_{\textrm n} the traction at the boundary. Please note that all surfaces that are not manually defined will automatically be free of traction. The boundary conditions can be defined in the xml-file within the <bcsAndLoads>
tag.
<bcsAndLoads>
<fix name="">
<comp dof=""/>
</fix>
<displacement name="">
<comp dof="" value=""/>
</displacement>
<pressure name="" value=""/>
<normalStiffness name="" volumeRegion="" value=""/>
<thermalStrain name="" value=""/>
</bcsAndLoads>
- fix: Fixes the nodes of the given geometry where only the direction (e.g. "x", "y", "z") specified by the tag
<comp dof ="">
is locked. This allows to fix the geometry in one or more directions (homogeneous Dirichlet value \mathbf{u} = \mathbf{0}). - displacement: Similar to
<fix>
, but instead of setting the value of the displacement to zero a certain value can be specified (inhomogeneous Dirichlet value \mathbf{u} = \mathbf{u}_{\textrm e}). - pressure: Prescribes a value for a uniform pressure applied to a surface (inhomogeneous Neumann value \mathbf{\sigma}\cdot\mathbf{n} = \mathbf{\sigma}_{\textrm n}).
- normalStiffness: Prescribes a stiffness acting in normal direction to a surface.
- thermalStrain: Used to prescribe a thermal strain on a volume.
Sources¶
Besides providing values for the boundary conditions the RHS can be used. For the mechanic PDE we can provide the given sources via the xml within the <bcsAndLoads>
tag.
<bcsAndLoads>
<force name="">
<comp dof="" value=""/>
</force>
<forceDensity name="">
<comp dof="" value=""/>
</forceDensity>
<volumeAcceleration name="">
<comp dof="" value=""/>
</volumeAcceleration>
</bcsAndLoads>
- force: Used to prescribe a force on a volume, where the direction is specified by the tag
<comp dof ="">
(e.g. "x", "y", "z"). - forceDensity: Similar to
<force>
but with a force-density instead of an absolute force. - volumeAcceleration: Applies a RHS based on a specified acceleration vector. Can be used, for example, to apply gravity in a static simulation by specifying the gravitational constant g = 9.80665 m/s^2 in the desired direction.
Pre-Stress¶
There is also a possibility to prescribe a pre-stress to your computation. There are two ways in the <bcsAndLoads>
to prescribe the pre-stress:
<bcsAndLoads>
<preStress region="">
<computeLHS sequenceStep=""/>
</preStress>
<preStress>
<prescribedLHS value="s0 s1 s2 s3 s4 s5"/>
</preStress>
</bcsAndLoads>
sequenceStep
.
* prescribedLHS: Allows to specify an individual stress tensor, where the s0 ... s5 variables have to be replaced by the numbers defining the respective components of the tensor, assuming Voigt notation.
Material¶
In the general case the material used for the simulaiton is anisotropic which can be defined in the material-xml as
<material name="material2">
<mechanical>
<density>
<linear>
<real> 1.0e+3 </real>
</linear>
</density>
<elasticity>
<linear>
<tensor dim1="6" dim2="6">
<real>
1.34615E+03 5.76923E+02 5.76923E+02 0.00000E+00 0.00000E+00 0.00000E+00
5.76923E+02 1.34615E+03 5.76923E+02 0.00000E+00 0.00000E+00 0.00000E+00
5.76923E+02 5.76923E+02 1.34615E+03 0.00000E+00 0.00000E+00 0.00000E+00
0.00000E+00 0.00000E+00 0.00000E+00 3.84615E+02 0.00000E+00 0.00000E+00
0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 3.84615E+02 0.00000E+00
0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 3.84615E+02
</real>
</tensor>
</linear>
</elasticity>
</mechanical>
</material>
For this example, all 21 individual components of the elasticity tensor can be defined. For a simpler case, e.g. an isotropic material, the definition of the elasticity reduces to
<elasticity>
<linear>
<isotropic>
<elasticityModulus>
<real> 3e7 </real>
</elasticityModulus>
<poissonNumber>
<real> 0.2 </real>
</poissonNumber>
</isotropic>
</linear>
</elasticity>
where the stiffness tensor is defined via
Here, the Voigt notation as well as the definition for \mu and \lambda given by
has been used. In these equations, E enotes the elasticity modulus and \nu the poisson ration.
Damping¶
Rayleigh Damping¶
A common approach to consider damping is the Rayleigh damping model, where a velocity proportional damping term with a damping matrix \mathbf{D} defined as a linear combination
of the mass matrix \mathbf{M} and the stiffness matrix \mathbf{K} is added to the FE equation as
A Rayleigh type damping can be added to the system in multiple ways. The 'classic' Rayleigh damping assumes constant coefficients \alpha and \beta, wich results in a frequency-dependent damping curve. Defining the so-called loss tangens delta function \mathrm{tan}\Delta, which is two times the damping ratio, the damping of a mode i can be shown according to follow
where \omega_i is the resonance frequency of the damped mode.
To employ the classic Rayleigh damping, at first a dampingId
has to be added to the damped region in the region list and the corresponding damping must be defined in a damping list below, e.g.,
<regionList>
<region name="solid_damped" dampingId="Rayl"/>
</regionList>
<dampingList>
<rayleigh id="Rayl"/>
</dampingList>
<damping>
<rayleigh>
<coefficients>
<alpha> 30 </alpha>
<beta> 0.00002 </beta>
</coefficients>
</rayleigh>
</damping>
The other two options to define the material damping are by specifying a loss-tangens delta \mathrm{tan} \Delta at specific frequencies. The two Rayleigh coefficients can in this case be determined by solving a system of two equations,
The \mathrm{tan}\Delta can either be specified either on two distinct frequency points,
<damping>
<rayleigh>
<computeFromTwoPoints>
<point1 frequency="90" lossTangensDelta="0.01"/>
<point2 frequency="110" lossTangensDelta="0.01"/>
</computeFromTwoPoints>
</rayleigh>
</damping>
<damping>
<rayleigh>
<computeFromTanDeltaAtFrequency>
<centerFreq> 100 </centerFreq>
<ratioDeltaF> 1e-1 </ratioDeltaF>
<lossTangensDelta> 0.01 </lossTangensDelta>
</computeFromTanDeltaAtFrequency>
</rayleigh>
</damping>
ratioDeltaF
is the percentage of variation from the center frequency centerFreq
for giving the upper/lower bounds of the band. If the ratioDeltaF
is kept sufficiently small it is thus also possible to specify an approximate damping at only a single frequency.
Frequency-Adapted Rayleigh Damping¶
For time-harmonic computations, the Rayleigh coefficients can easily be adjusted in each step to obtain a desired damping for the corresponding frequencies. In this case, the \alpha and \beta coefficients are re-computed internally at each frequency to satisfy a specified \mathrm{tan}\Delta. To enable this type of damping, specify an adaptedLossTangensDelta
damping in the simulation XML, e.g.,
<regionList>
<region name="domain1" dampingId="adaptedLoss"/>
</regionList>
<dampingList>
<adaptedLossTangensDelta id="adaptedLoss"/>
</dampingList>
mathParser
. This is enabled, e.g., by specifying the file using the sample1D
interpolator:
The file containing the damping information might look like
where the first column contains the frequencies and the second the corresponding \mathrm{tan}\Delta values.
Note: For Eigenfrequency analysis, the quadratic Eigenvalue problem must be solved, if damping should be considered (see definition below). For detailed information see [^Kaltenbacher2010].
Analysis Types¶
Since in general we are dealing with a time-dependent PDE, we can distinguish three analysis types:
- Transient: \frac{\partial }{\partial t}\neq 0
-
Harmonic: \partial / \partial t (\cdot) = j \omega (\cdot)
-
Eigenfrequency: \partial / \partial t (\cdot) = j \omega (\cdot) and RHS = 0
Postprocessing Results¶
Node Results¶
-
Mechanic Displacement: Displacement vector for each node. (primary solution)
-
Mechanic Velocity: Velocity vector for each node. \begin{equation} \mathbf{v} = \frac{\partial \mathbf{u}}{\partial t} \end{equation}
-
Mechanic Acceleration: Acceleration vector for each node. \begin{equation} \mathbf{a} = \frac{\partial^2 \mathbf{u}}{\partial t^2} \end{equation}
-
Mechanic RHS load: Right hand side source term in the FE formulation as nodal loads
Element Results¶
-
Mechanic Strain \begin{equation} \mathbf{s} = \frac{1}{2} ( \mathbf{u} + \mathbf{u}^{\textrm{T}}). \end{equation}
-
Mechanic Stress \begin{equation} \mathbf{\sigma} = \mathbf{C} : \mathbf{s}, \end{equation}
-
Mechanic Kinetic Energy Density \begin{equation} e_\mathrm{kin} = \frac{1}{2} \rho \mathbf{v}^{\textrm{T}} \cdot \mathbf{v} \end{equation}
-
Von Mises stress: Equivalent stress according to von Mises. \begin{equation} \sigma_\mathrm{v} = \sqrt{\frac{1}{2} \left( (\sigma_{11}-\sigma_{22})^2 + (\sigma_{22}-\sigma_{33})^2 + (\sigma_{33}-\sigma_{11})^2 \right) + 3 (\sigma_{12}^2 + \sigma_{23}^2 + \sigma_{31}^2)} \end{equation}
Region Results¶
-
Mechanic Kinetic Energy \begin{equation} E_\mathrm{kin} = \int_{{\Omega}} e_\mathrm{kin} \, dV \end{equation}
-
Mechanic Deformation Energy \begin{equation} E_\mathrm{def} = \frac{1}{2} <u,Ku> \end{equation}