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>
        </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.

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

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) \, .$$

For doing so, the damping-ID has to be added to the damped region in the region list and the corresponding a damping list including the damping-ID needs to be added as follows:

        <regionList>
          <region  dampingId="Rayl" name="solid_damped"/>
        </regionList>   
        <dampingList>
          <rayleigh id="Rayl"/>
        </dampingList>

The mass proportional damping coefficient $\alpha$ and the stiffness proportional damping coefficient $\beta$ must be provided in a dedicated tag damping in the material definition (XML) in addition to the density and elasticity tags as

                <rayleigh>
                    <alpha>4</alpha>
                    <beta>8</beta>
                    <measuredFreq>0</measuredFreq>
                </rayleigh>

, where the measured frequency ($f_M$) can be used in a harmonic analysis (for transient analysis it is ignored) when the damping coefficients should be adjusted depending on the frequency $f$ as $\alpha' = \alpha \frac{f}{f_M}$ and $\beta' = \beta \frac{f_M}{f}$. This is must be defined in the damping list in in the simulation-XML as:

        </dampingList>
          <rayleigh id="Rayl">
            <adjustDamping>yes</adjustDamping>
            <ratioDeltaF>2</ratioDeltaF>
          </rayleigh>
        </dampingList>

An alternative way to directly providing the damping coefficients is to is to provide the loss factor $\tan\delta$ and the measured frequency $f_M$ in the material definition (XML) as

                <rayleigh>
                    <lossTangensDelta>10</lossTangensDelta>
                    <measuredFreq>2</measuredFreq>
                </rayleigh>

, where the loss factor $\tan \delta_{i}$ of the mode $i$ relates to the damping ratio $\xi_{i}$ as

$$\tan \delta_{i}=2 \xi_{i}=\frac{\alpha + \beta \omega_{i}^{2}}{\omega_{i}} \, ,$$

resulting from

$$\alpha + \beta \omega_{i}^{2}=2 \omega_{i} \xi_{i} \ .$$

The Rayleigh damping parameters $\alpha$ and $\beta$ are defined from the following equations:

$$\alpha + \beta (\omega_{i} + \Delta \omega)^{2}=2 (\omega_{i} + \Delta\omega) \xi_{i} \, ,$$

$$\alpha + \beta ( \omega_{i} - \Delta \omega)^{2}=2 ( \omega_{i} - \Delta\omega) \xi_{i} \, .$$

$\Delta\omega$ is kept small to meet prescribed $\xi_{i}$ and $\omega_{i}$.

Note: For Eigenfrequency analysis, the quadratic Eigenvalue problem must be solved, if damping should be considered (see definition below). For detailed information see ¹.

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

---

1. M. Kaltenbacher. Numerical Simulation of Mechatronic Sensors and Actuators: Finite Elements for Computational Multiphysics. Springer Berlin Heidelberg, 2015. ISBN 978-3-642-40169-5. URL: https://www.springer.com/de/book/9783642401695. ↩