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

$\begin{equation} \mathbf{u}(t) = \mathbf{u}_{\textrm e}(t) \qquad \textrm{on } \Gamma_{\textrm e}, \end{equation}$

$\begin{equation} \mathbf{\sigma}(t)\cdot\mathbf{n} = \mathbf{\sigma}_{\textrm n}(t) \qquad \textrm{on } \Gamma_{\textrm n}, \end{equation}$

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>

* computeLHS: Picks up the stresses from a previous 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

$\begin{equation} \mathbf{C} = \begin{bmatrix} 2 \mu + \lambda & \lambda & \lambda & 0 & 0 & 0\\ \lambda & 2 \mu + \lambda & \lambda & 0 & 0 & 0\\ \lambda & \lambda & 2 \mu + \lambda & 0 & 0 & 0\\ 0 & 0 & 0 & \mu & 0 & 0\\ 0 & 0 & 0 & 0 & \mu & 0\\ 0 & 0 & 0 & 0 & 0 & \mu \end{bmatrix}. \end{equation}$

Here, the Voigt notation as well as the definition for $\mu$ and $\lambda$ given by

$\begin{equation} \mu = \frac{E}{2(1+\nu)} \end{equation}$

$\begin{equation} \lambda = \frac{E \nu}{(1+\nu)(1-2\nu)} \end{equation}$

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

$\mathbf{D} = \alpha \mathbf{M} + \beta \mathbf{K}$

of the mass matrix $\mathbf{M}$ and the stiffness matrix $\mathbf{K}$ is added to the FE equation as

$\mathbf{M} \mathbf{\ddot u}(t) + \mathbf{D} \mathbf{\dot u}(t) +\mathbf{K} \mathbf{u}(t) = \mathbf{F}(t) \, .$

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

$\mathrm{tan}\Delta_i = 2 \zeta_i = \frac{\alpha}{\omega_i} + \beta \omega_i,$

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>

Additionally, the Rayleigh damping properties must be defined in the material XML file. This can currently be done in three different ways. First, the

$\alpha$ and

$\beta$ coefficients can be defined explicitly:

<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,

$\begin{bmatrix} \frac{1}{\omega_1} & \omega_1 \\ \frac{1}{\omega_2} & \omega_2 \end{bmatrix} \cdot \begin{bmatrix} \alpha \\ \beta \end{bmatrix} = \begin{bmatrix} \tan \delta_1 \\ \tan \delta_2 \end{bmatrix}.$

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>

or on the two edge frequencies of a frequency band,

<damping>
  <rayleigh>
    <computeFromTanDeltaAtFrequency> 
      <centerFreq> 100 </centerFreq>
      <ratioDeltaF> 1e-1 </ratioDeltaF>
      <lossTangensDelta> 0.01 </lossTangensDelta>
    </computeFromTanDeltaAtFrequency>
  </rayleigh>
</damping>

where 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>

Additionally, the damping property must be specified in the material XML file. Here, two options are available. First, it is possible to pass a constant non-negative float value for the

$\mathrm{tan}\Delta$ :

<damping>
  <adaptedLossTangensDelta value="0.01"/>
</damping>

The second option is to provide a file containing a specified

$\mathrm{tan}\Delta$ and corresponding frequencies. This file can be read by openCFS and interpolated to the simulated frequencies via the mathParser. This is enabled, e.g., by specifying the file using the sample1D interpolator:

<damping>
  <adaptedLossTangensDelta value="sample1D('./frequencyDamping.txt',f,1)"/>
</damping>

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$

    <analysis>
      <transient>
        <numSteps>10</numSteps>
        <deltaT>0.1</deltaT>
      </transient>
    </analysis>

Harmonic: $\partial / \partial t (\cdot) = j \omega (\cdot)$

<analysis>
    <harmonic>
        <numFreq>10</numFreq>
        <startFreq>10</startFreq>
        <stopFreq>3000</stopFreq>
        <sampling>linear</sampling>
    </harmonic>
</analysis>

Eigenfrequency: $\partial / \partial t (\cdot) = j \omega (\cdot)$ and $RHS = 0$

    <analysis>
      <eigenFrequency>
        <isQuadratic>no</isQuadratic>
        <numModes>5</numModes>
        <freqShift>0</freqShift>
        <writeModes>yes</writeModes>        
      </eigenFrequency>
    </analysis>

Postprocessing Results¶

Node Results¶

Mechanic Displacement: Displacement vector for each node. (primary solution)
```
<nodeResult type="mechDisplacement">
```
Mechanic Velocity: Velocity vector for each node. \begin{equation} \mathbf{v} = \frac{\partial \mathbf{u}}{\partial t} \end{equation}
```
<nodeResult type="mechVelocity">
```
Mechanic Acceleration: Acceleration vector for each node. \begin{equation} \mathbf{a} = \frac{\partial^2 \mathbf{u}}{\partial t^2} \end{equation}
```
<nodeResult type="mechAcceleration">
```
Mechanic RHS load: Right hand side source term in the FE formulation as nodal loads
```
<nodeResult type="mechAcceleration">
```

Element Results¶

Mechanic Strain \begin{equation} \mathbf{s} = \frac{1}{2} ( \mathbf{u} + \mathbf{u}^{\textrm{T}}). \end{equation}
```
<elemResult type="mechStrain">
```
Mechanic Stress \begin{equation} \mathbf{\sigma} = \mathbf{C} : \mathbf{s}, \end{equation}
```
<elemResult type="mechStress">
```
Mechanic Kinetic Energy Density \begin{equation} e_\mathrm{kin} = \frac{1}{2} \rho \mathbf{v}^{\textrm{T}} \cdot \mathbf{v} \end{equation}
```
<elemResult type="mechKinEnergyDensity">
```
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}
```
<elemResult type="vonMisesStress">
```

Region Results¶

Mechanic Kinetic Energy \begin{equation} E_\mathrm{kin} = \int_{{\Omega}} e_\mathrm{kin} \, dV \end{equation}
```
<regionResult type="mechKinEnergy">
```
Mechanic Deformation Energy \begin{equation} E_\mathrm{def} = \frac{1}{2} <u,Ku> \end{equation}
```
<regionResult type="mechKinEnergy">
```

Tutorials¶

Cantilever