Skip to content

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 coe´Čâcient \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 1.

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