Elastodynamics: 1D Wave Propagation in an Elastic Solid
Overview
The objective of this case study is to verify FEBio’s implicit and explicit structural solver in a dynamic analysis, against an analytical solution for 1D wave propagation in an elastic solid medium under infinitesimal strains. This case study also provides an opportunity to investigate energy conservation in a solid dynamics analysis, as well as to compare the performances of FEBio’s implicit and explicit solvers.
Elastodynamics in FEBio
Governing Equation
As presented in the Theory Manual, in solid dynamics FEBio solves the linear momentum balance equation:
| (1) | $\operatorname{div}\boldsymbol{\sigma}+\rho\mathbf{b}=\rho\mathbf{a}$ |
where $\boldsymbol{\sigma}$ is the Cauchy stress tensor in the medium, $\rho$ is its mass density, $\mathbf{b}$ is the external specific body force and $\mathbf{a}$ is the medium’s acceleration. Our first goal is to solve this equation analytically for a one-dimensional wave propagation problem, then replicate that problem in FEBio and solve it numerically.
1D Wave Propagation Analysis
Consider that the 1D analysis is performed along the $y-$direction, with the unit basis vector $\mathbf{e}_{y}$ oriented upward. For a 1D analysis we consider that the displacement vector has the form $\mathbf{u}=u\mathbf{e}_{y}$, therefore the acceleration vector is $\mathbf{a}=\frac{d^{2}u}{dt^{2}}\mathbf{e}_{y}$. The specific body force is assumed to be gravity, pointing along the negative $y-$direction, so that $\mathbf{b}=-g\mathbf{e}_{y}$. The stress is assumed to be one-dimensional as well (i.e., we neglect the lateral expansion/contraction of our material domain relative to its axial deformation, which is equivalent to assuming that Poisson’s ratio is zero), thus $\boldsymbol{\sigma}=\sigma\mathbf{e}_{y}\otimes\mathbf{e}_{y}$. Finally, we use Hooke’s law to relate the stress to the strain according to $\sigma=E\frac{du}{dy}$ where $E$ is Young’s modulus. Substituting these relations into eq.(1) produces the classical wave equation,
| (2) | $\frac{d^{2}u}{dy^{2}}-\frac{1}{c^{2}}\frac{d^{2}u}{dt^{2}}=\frac{1}{c^{2}}g$ |
where
| (3) | $c=\sqrt{\frac{E}{\rho}}$ |
is the speed of sound in the elastic solid medium.
Essentially, we are solving the problem of a prismatic bar whose long axis is along the $y-$ direction, which is suddenly subjected to the gravitational force, producing an oscillatory response in that bar. Let the height of the bar be given by $h$. Consider that the bar is prevented from moving along $y$ at $y=0$, while the stress is zero, $\sigma=0$, at the top of the bar, $y=h$. Therefore, the boundary conditions for the partial differential equation in (2) are
| (4) | $u\left(0,t\right)=0\,,\quad\sigma\left(h,t\right)=E\frac{du}{dy}\left(h,t\right)=0$ |
The initial conditions for this problem are that the displacement and velocity are zero,
| $u\left(y,0\right)=0\,,\quad\frac{du}{dt}\left(y,0\right)=0$ | (5) |
We can solve the PDE of eq.(2), subject to the boundary conditions of eq.(4) and initial conditions of eq.(5), using the method of Laplace transforms. In Laplace transform space, the PDE reduces to an ordinary differential equation of the form
| (6) | $\frac{d^{2}\bar{u}}{dy^{2}}-\frac{s^{2}}{c^{2}}\bar{u}=\frac{g}{sc^{2}}$ |
where $\bar{u}\left(y,s\right)=\mathcal{L}\left\{ u\left(y,t\right)\right\}$, subject to the boundary conditions
| (7) | $\bar{u}\left(0,s\right)=0\,,\quad\frac{d\bar{u}}{ds}\left(h,s\right)=0$ |
The solution to this ODE is
| (8) | $\bar{u}\left(y,s\right)=-\frac{g}{s^{3}}+\frac{\cosh\frac{s}{c}\left(h-y\right)}{\cosh\frac{s}{c}h}\frac{g}{s^{3}}$ |
Using tables of Laplace transforms, we can invert this expression to produce the final analytical solution for the displacement in the bar,
| (9) | $\frac{u\left(y,t\right)}{g\tau^{2}}=\frac{1}{2}\frac{y}{h}\left(\frac{y}{h}-2\right)-\frac{2}{\pi^{3}}\sum_{n=1}^{\infty}\frac{\left(-1\right)^{n}}{\left(n-\frac{1}{2}\right)^{3}}\cos\left(n-\frac{1}{2}\right)\pi\left(1-\frac{y}{h}\right)\cos\left(n-\frac{1}{2}\right)\pi\frac{t}{\tau}$ |
where
| (10) | $\tau=\frac{h}{c}$ |
is the time required for a wave to propagate through the entire height $h$ of the bar. Note that the form of eq.(9) embodies the natural normalization of the variables in this problem. Thus,
| (11) | $\hat{y}\equiv\frac{y}{h}\,,\quad\hat{t}\equiv\frac{t}{\tau}\,,\quad\hat{u}\equiv\frac{u}{g\tau^{2}}$ |
respectively represent the normalized spatial coordinate, normalized time, and normalized axial displacement.
FEBio Simulation
Geometry, Mesh, Material Properties and Boundary Conditions
We analyze this problem using a cube with sides $h=1\,\text{m}$. The cube is meshed with 100 uniform hex8 elements along the $y-$direction, and one element along $x-$ and $z-$directions. The elastic solid is taken to be neo-Hookean, with $E=10^6\,\text{Pa}$, $\nu=0$, and $\rho=1000\,\text{kg}/\text{m}^3$. The displacements are fixed on the left ($u_x=0$), bottom ($u_y=0$) and back ($u_z=0$) faces. Based on these properties, we calculate the speed of sound in this solid using eq.(3) as $c=31.62\,\text{m}/\text{s}$, from which we can evaluate the time $\tau$ using eq.(10) as $\tau=0.03162\,\text{s}$. We plan for a total analysis time of $\sim20\tau=0.64\,\text{s}$. For simplicity, we set $g=10\,\text{m}/\text{s}^2$ and we assign a step load curve to this body load. (Alternatively, the load curve assigned by default to the body force can be removed.)
Implicit Solver
The standard FEBio solver for structural mechanics is an implicit solver, implying that it factors the inertia term $\rho\mathbf{a}$ into the stiffness matrix. It allows us to use relatively large time increments without significantly compromising the accuracy of an analysis. For this 1D wave propagation problem we select a time increment $\Delta t\sim\tau/20=0.0016\,\text{s}$, such that the total number of steps is $400$.
As explained in the Theory Manual, the standard solution scheme in finite element analyses of structural problems only solves the momentum balance equation, as given in eq.(1). However, when performing elastodynamic analyses, it is expected that ideal elastic bodies should not dissipate any energy. In practice however, the choice of numerical discretization scheme may lead to lesser or greater amounts of numerical energy dissipation. To overcome this issue, one may choose to employ an energy-momentum conservation scheme, such as the one described in the Theory Manual. This scheme is automatically employed when setting the solver parameter called the spectral radius ($\rho_\infty$) to $1$. In the Generalized-$\alpha$ scheme, this parameter varies in the range $0\le\rho_\infty\le1$. In FEBio, we can also set $\rho_\infty=-1$ for the backward Euler integration scheme, or $\rho_\infty=-2$ to recover the trapezoidal integration scheme, which used to be the only available option in older versions of FEBio.
For this case study we set $\rho_\infty=1$, although it can be verified by trial and error that some of the other available integration schemes can perform equally well for this particular problem. The solver settings for this analysis are summarized below:
<Module type="solid"/>
<Control>
<analysis>DYNAMIC</analysis>
<time_steps>400</time_steps>
<step_size>0.0016</step_size>
<solver>
<max_refs>25</max_refs>
<max_ups>0</max_ups>
<diverge_reform>1</diverge_reform>
<reform_each_time_step>1</reform_each_time_step>
<dtol>0.001</dtol>
<etol>0.01</etol>
<rtol>0.001</rtol>
<lstol>0.9</lstol>
<min_residual>1e-20</min_residual>
<qnmethod>BFGS</qnmethod>
<rhoi>1</rhoi>
</solver>
<time_stepper>
<dtmin>1.6e-05</dtmin>
<dtmax>0.0016</dtmax>
<max_retries>5</max_retries>
<opt_iter>10</opt_iter>
</time_stepper>
</Control>We then compare the FEBio solution to the analytical solution of eq.(9), using the normalization of variables presented in eq.(11).
Results shown in Figure 1 confirm that FEBio evaluates the solution to this problem accurately. As expected, the FEBio solution starts to slightly deviate from the analytical solution with increasing time, since temporal discretization errors accumulate over time.
Explicit Solver
This case study also serves as a convenient testbed for the structural mechanics explicit solver implemented in FEBio. In an explicit solver, the inertia term $\rho\mathbf{a}$ is typically evaluated by lumping its value at the nodes of the mesh, then discretizing the acceleration using a forward difference scheme. This makes it possible to explicitly evaluate the unknown nodal displacements at the current time/iteration using the solutions at the previous time/iteration (i.e., without having to solve a system of equations simultaneously for all unknowns). Explicit schemes do not converge unconditionally. Thus, for a given mesh, it is important to select a sufficiently small time increments to produce a converged (and reasonably accurate) solution.
If we tried to use the same time step size $\Delta t$ as for the implicit solver, we would find that the explicit solver would fail at the very first time step, with a negative jacobian error. For this problem, the critical time step can be evaluated from,
| (10) | $\tau_c=\frac{h_e}{c}$ |
Here, $h_e$ is the element size, which in this case is 0.01m. From this, we can calculate that the critical time step is, $\tau_c=0.0003162\,\text{s}$.
Therefore, for this problem, we choose to employ a time step size $\Delta t\sim\tau/200=1.6\times 10^{-4}$, requiring a total of $4000$ steps. In effect, we reduce the step size by a factor of 10 relative to that employed with the implicit solver. However, the solution time of the explicit solver remains shorter than that of the implicit solver, since the explicit solver does not solve a system of equations simultaneously. The settings for the explicit solver are shown below.
<Module type="explicit-solid"/>
<Control>
<analysis>DYNAMIC</analysis>
<time_steps>4000</time_steps>
<step_size>0.00016</step_size>
<plot_zero_state>0</plot_zero_state>
<plot_range>0,-1</plot_range>
<plot_level>PLOT_MAJOR_ITRS</plot_level>
<output_level>OUTPUT_MAJOR_ITRS</output_level>
<plot_stride>1</plot_stride>
<adaptor_re_solve>1</adaptor_re_solve>
<solver>
<symmetric_stiffness>symmetric</symmetric_stiffness>
<equation_scheme>staggered</equation_scheme>
<equation_order>default</equation_order>
<optimize_bw>0</optimize_bw>
<mass_lumping>1</mass_lumping>
<dyn_damping>1</dyn_damping>
</solver>
</Control>
The explicit solver solution is compared to the analytical solution in Figure 2. Here again, we find that the FEBio solution agrees with the analytical solution, thus verifying the accuracy of the FEBio explicit solver.
Summary
The objective of this case study was to verify FEBio’s implicit and explicit structural solver in a dynamic analysis, against an analytical solution for 1D wave propagation in an elastic solid medium under infinitesimal strains. An exact solution was obtained for a 1D elastodynamic problem under infinitesimal strains, as presented in eq.(9). Results presented in Figures 1 and 2 verify that implicit and explicit solvers implemented in FEBio produce accurate solutions. This case study illustrates FEBio’s ability to model elastodynamic problems.
Suggested Exercise
It is not necessary to discretize the mesh with 100 elements along the $y-$axis. On your own, remesh this model using four hex27 elements along $y$, and show that the implicit and explicit solutions can produce essentially the same answers as the examples above. Keep in mind that when changing the element size, the critical time step for the explicit solver is also changed.