3.2. Governing Equations

3.2.1. Generalized Conservation Equation

We first introduce a general conservation equation, as a model for the specific equations that Aria solves, demonstrating how the Galerkin finite element method is applied to it, and how the integration by parts is carried out on its individual terms. Following [3], the conservation of a general scalar quantity $\symSlnVar(\x,t)$ , with units of amount-per-unit-volume, at a point $\x$ and time $t$ can be expressed as

(3.1) $\pt{\symSlnVar} + \div(\symSlnVar\v) = - \div\f + \symSlnVar_\Vol$

where $\v$ is the mass average velocity, $\f$ is the diffusive flux of $\symSlnVar$ , and $\symSlnVar_\Vol$ is the volumetric source of $\symSlnVar$ . In Equations Aria Solves, the physical meaning of the source term of each equation is discussed.

The Galerkin FEM (G/FEM) residual form of (3.1) is formed by bringing the right hand side terms to the left, multiplying by the FEM weight function $\symWeightFunc^i$ and integrating over the volume $\Vol$ ,

(3.2) $\symRes_\symSlnVar^i = \int\limits_\Vol\left(\pt{\symSlnVar} + \v\bcdot\grad \symSlnVar + \symSlnVar\div\v + \div\f - B_\Vol\right)\symWeightFunc^i\dV = 0.$

In many applications $\div\v=0$ so we ignore that term from here on. However, it is straight forward to account for this term via the source term $B_\Vol$ . Using the vector identity $(\div\f)\symWeightFunc^i = \div(\f\symWeightFunc^i) - \grad\symWeightFunc^i\bcdot\f$ and using the divergence theorem, (3.2) becomes

(3.3) $\symRes_\symSlnVar^i = \int\limits_\Vol\left[\left(\pt{\symSlnVar} + \v\bcdot\grad \symSlnVar - B_\Vol\right)\symWeightFunc^i - \grad\symWeightFunc^i\bcdot\f\right]\dV + \int\limits_\Surf \n\bcdot\f\symWeightFunc^i\dS = 0.$

Here $\n$ is a unit normal along the boundary $\Surf$ , pointing out of the volume $\Vol$ .

Equation (3.3) embodies the sign convention for sources, fluxes and equation terms used within Aria. For example, scalar flux expressions in Aria provide values for $f_n \equiv \n\bcdot\f$ and should be positive for a flux of $\symSlnVar$ leaving the volume $\Vol$ .

Note also that we have not assigned a units convention to the equation. Any unit system may be employed in the specification of the individual terms in (3.1). However, each term in (3.1) must have overall units of [ $\symSlnVar$ ] / [time], and the overall residual expression has units of [ $\symSlnVar$ ] * [L]³/ [time], where [ $\symSlnVar$ ] are units of the conserved quantity, $\symSlnVar$ , [L] is the unit of the length scale, and [time] is the unit for time.

3.2.2. Statement of the Transient Problem

We are now in a position to state mathematically the initial–boundary value problem described at the beginning of this chapter. Let the domain $\Omega$ consist of $N$ non-overlapping subdomains. Let $\mathcal{T}$ be the set of all subdomains. Then the statement of the boundary value problem becomes, find the solution $\symSlnVar = \symSlnVar(\bfx, t)$ , which satisfies

(3.4) $\dfrac{\partial \symSlnVar}{\partial t} + \v \cdot \grad \symSlnVar - \div \f = \symSlnVar_V && \forall\; \bfx \in \{ \Omega_i \,|\, \Omega_i \in \mathcal{T} \}; t > t_0 \\ \symSlnVar(\bfx,t_0) = \symSlnVar_0(\bfx) && \forall\; \bfx \in \Omega \\ \symSlnVar(\bfx, t) = \symSlnVar_b(\bfx, t) && \forall\; \bfx \in \Gamma_\symSlnVar \\ f_n = f_b(\bfx,t,\symSlnVar) && \forall\; \bfx \in \partial\Omega\setminus \Gamma_\symSlnVar \\ \ljump f_n \rjump = 0 && \forall\; \bfx \in \{\partial\Omega_{i-k} \,|\, \Omega_i,\Omega_k \in \mathcal{T}\}$

Here $f_n$ is the component of flux normal to a surface, where $f_n = \f \cdot \n$ , and $n$ is the outward unit normal vector. The notation $\partial\Omega\setminus \Gamma_\symSlnVar$ indicates the complement of $\Gamma_\symSlnVar$ in $\partial\Omega$ ,

$\partial\Omega\setminus \Gamma_\symSlnVar = \left\{\bfx \in \partial\Omega \mid \bfx \notin \Gamma_\symSlnVar\right\},$

or the boundary of $\Omega$ excluding the surface $\Gamma_\symSlnVar$ . In order for the problem to be well posed, it is not possible to specify both the flux and the temperature at the same location. Note that the initial condition $\symSlnVar_0$ , the boundary condition $\symSlnVar$ , and the flux boundary condition $f_b$ are usually specified in a piecewise manner over the subdomains and their boundaries. Various forms of the specified flux function $f_b(\bfx, t, \symSlnVar)$ are possible and will be described later in this chapter. For example, an adiabatic condition is specified if $f_b(\bfx, t, \symSlnVar) = 0$ . The notation $\ljump f_n \rjump$ represents the jump in the normal flux across a surface between two subdomains.

3.2.3. Statement of the Stationary Problem

For stationary or steady-state problems (3.4) may be simplified since the time derivative vanishes by definition. Accordingly, we obtain

(3.5) $\v \cdot \grad \symSlnVar - \div \f = \symSlnVar_V && \forall\; \bfx \in \{ \Omega_i \,|\, \Omega_i \in \mathcal{T} \}; t > t_0 \\ \symSlnVar(\bfx) = \symSlnVar_b(\bfx) && \forall\; \bfx \in \Gamma_\symSlnVar \\ f_n = f_b(\bfx,\symSlnVar) && \forall\; \bfx \in \partial\Omega\setminus \Gamma_\symSlnVar \\ \ljump f_n \rjump = 0 && \forall\; \bfx \in \{\partial\Omega_{i-k} \,|\, \Omega_i,\Omega_k \in \mathcal{T}\}$

3.2.4. Equations Aria Solves

Below, a more detailed description is provided for each equation currently supported by Aria:

3.2.5. Notes on Solid Mechanics

Some of the standard references on solid mechanics include [4], [5], [6] and [7]. As is often the case, the mathematical notion used through-out these texts is different in many cases and this is often a source of confusion. Here, we’ll lay out some basic definitions in our notation and, when possible, give the notation used in these other texts.

In what follows, $\x$ is the position vector of a material particle in the deformed or current spatial configuration and $\X$ is the position vector of a material particle in the undeformed or initial or reference configuration. The displacement vector $\d$ is the difference between these to states 1 viz. $\x = \X + \d$ .

We will make extensive use of the gradients of these fields and so it is important to distinguish between gradients in the reference configuration and the current configuration. Gradients in the current configuration are denoted $\grad$ in Gibbs notation or $\partial\ /\partial x_i$ in index notation. Gradients in the reference configuration are denoted $\GRAD$ in Gibbs notation or $\partial\ /\partial X_i$ in index notation 2.

Next, we define the deformation gradient $\F$

$\F &\equiv \GRAD\x^t \\ &= \GRAD\X^t + \GRAD\d^t \\ &= \I + \GRAD\d^t$

where the superscript $^t$ denotes the transpose operator 3. The inverse deformation gradient 4, $\F^{-1}$ , is also useful and can be computed as

$\F^{-1} &\equiv \grad\X^t \\ &= \grad\x^t - \grad\d^t \\ &= \I - \grad\d^t.$

The determinants of $\F$ and $\F^{-1}$ are denoted $J$ and $J^{-1}$ respectively and are often used in transformations between different stress definitions 5.

It’s worth noting at this point that in Aria both gradient operators, $\grad$ and $\GRAD$ , are available as Expression objects as are $\F$ , $\F^{-1}$ , $J$ and $J^{-1}$ .

The Green or Green-Lagrange strain tensor $\E$ is defined 6 as

$\E &\equiv \frac{1}{2}\left( \F^t\bcdot\F - \I \right) \\ &= \frac{1}{2}\left(\GRAD\d + \GRAD\d^t + \GRAD\d\bcdot\GRAD\d^t\right).$

The Green strain is a strain measure in the reference configuration and is suitable for large deformations and large rotations. The analogous Eulerian (or Almansi’s) strain tensor $\E^*$ is defined 7 as

$\E^* &\equiv \frac{1}{2}\left( \I - \F^{-t}\bcdot\F^{-1}\right) \\ &= \frac{1}{2}\left(\grad\d + \grad\d^t - \grad\d\bcdot\grad\d^t\right).$

The Eulerian strain is also a suitable strain measure for large deformations and rotations but is defined in the current configuration.

The Cauchy stress, $\tensor{\sigma}$ , is a stress measure defined in the current configuration as

$\tensor{\sigma} = \lambda E_{kk} \I + 2 \mu \E$

where $E_{kk}$ denotes the trace of $\E$ and $\lambda$ and $\mu$ are the Lamé coefficients. This constitutive equation may also be augmented with some initial residual stress or a thermal stress

$\tensor{\sigma} = \lambda E_{kk} \I + 2 \mu \E - \beta\left(T - T_\circ\right)\I + \tensor{\sigma}_r.$

Here $\beta$ is related to the coefficient of thermal expansion, $T$ is the temperature, $T_\circ$ is the reference temperature of the solid and $\tensor{\sigma}_r$ is the residual stress.

For large deformations and large rotations Aria uses the second Piola-Kirchhoff stress which is defined in the reference configuration and is related to the Cauchy stress as

$\tensor{S} = J \F^{-1}\bcdot\tensor{\sigma}\bcdot\F^{-t}.$

The reverse transformation is readily given by

$\tensor{\sigma} = J^{-1} \F\bcdot\tensor{S}\bcdot\F^t.$

Mathematically, the Cauchy stress $\tensor{\sigma}$ is most conveniently expressed in terms of the Lamé coefficients $\lambda$ , $\mu$ and $\beta$ . In practice, however, it is more common to measure and report a different but related set of parameters: the Young’s modulus $E$ , the Poisson’s ratio $\nu$ and the coefficient of thermal expansion $\alpha$ . (Note, the shear modulus $G = \mu$ .) The relationship between these two sets of parameters is

$2 \mu & = \frac{E}{\left(1+\nu\right)} \\ \lambda & = \frac{\nu E}{\left(1+\nu\right)\left(1-2\nu\right)} \;\; = \;\; 2 \mu \frac{\nu}{\left(1-2\nu\right)} \\ \beta & = \frac{\alpha E}{\left(1-2\nu\right)} \;\; = \;\; \alpha \left( 3 \lambda + 2 \mu \right)$

In Aria, there are separate Expression_Names for the Cauchy stress and the second Piola-Kirchhoff stress expressions. In the input files, users provide their choice of constitutive relations in the material model specification, e.g.,

Begin Aria Material The_Material
   Density     = Constant rho = 2.33e-15
   Mesh Lambda = Constant lambda = 52810.30445
   Mesh Two Mu = Constant two_mu = 134426.2295
   Mesh Stress = Nonlinear_Elastic Reference_Frame=Moving
   Mesh Stress = Residual Sx=-11 Sy=-11
   Mesh Stress = Isothermal T=800 T_ref=450
End

In this Example, there will be three contributions to the stress in the mesh stress: nonlinear elasticity, a planar residual stress and an isotropic linear thermal stress. Internally, Aria contains a separate expression for transforming these Cauchy stresses into Piola-Kirchhoff stresses, viz.

(3.92) $\tensor{S} = J \F^{-1}\bcdot\left(\sum\limits_i\tensor{\sigma}_i\right)\bcdot\F^{-t}.$

Note that second Piola-Kirchhoff stresses are not specified in the input file – only Cauchy stresses are. Aria automatically creates an Expression to compute the transform in equation (3.92).

Footnotes

1: $\d$ is denoted $\vector{u}$ in [4], [5], [6], and [7].
2: In [7] $\GRAD$ is denoted $\grad_0$ .
3: In [5] this is called the conjugate dyadic and is denoted with a subscript $c$ . In [4] the quantity $\GRAD\x^t$ is denoted $\x\stackrel{\leftharpoonup}{\grad}_X$ where the arrow over the gradient operator denotes the direction of the operation.
4: In [5] $\F^{-1}$ is denoted $\tensor{H}$ .
5: Sometimes $J$ is expressed as a ratio of the densities between the reference and current configurations, $\rho_\circ/\rho$ .
6: In [5] $\E$ is denoted $\tensor{L}_G$ and is called the Lagrangian strain.
7: In [5] $\E^*$ is denoted $\E_A$ .