@f]
where $\kappa := \lambda + 2/3 \mu$ is the bulk modulus ($\lambda$ and $\mu$ are the Lame parameters)
and $\overline{I}_1 := \textrm{tr}\ \overline{\mathbf{b}}$.
-The function $\mathcal{G}(J)$ is required to be strictly convex and satisfy the condition $\mathcal{G}(1) = 0$.
+The function $\mathcal{G}(J)$ is required to be strictly convex and satisfy the condition $\mathcal{G}(1) = 0$,
+among others, see Holzapfel (2001) for further details.
In this work $\mathcal{G}:=\frac{1}{4} [ J^2 - 1 - 2\textrm{ln}J ]$.
Incompressibility imposes the isochoric constraint that $J=1$ for all motions $\mathbf{\varphi}$.
&=
- \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V
+ \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V
- + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \,.
+ + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \\
+&=
+- \int_{\Omega_0} [\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} :\boldsymbol{\tau}~\mathrm{d}V
++ \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V
+ + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \, ,
@f}
+where
+$[\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} = 1/2[ \mathrm{grad}\delta\mathbf{u} + [\mathrm{grad}\delta\mathbf{u}]^T] $.
We will use an iterative Newton-Raphson method to solve the nonlinear residual equation $R$.
For the sake of simplicity we assume dead loading, i.e. the loading does not change due to the deformation.
A typical screen output generated by running the problem is shown below.
The particular case demonstrated is that of the $Q_2-DGPM_1-DGPM_1$ formulation.
It is clear that, using the Newton-Raphson method, quadratic convergence of the solution is obtained.
-Other than in the first time step, solution convergence is achieved within 5 Newton increments.
+Solution convergence is achieved within 5 Newton increments for all time-steps.
The converged displacement's $L_2$-norm is several orders of magnitude less than the geometry scale.
@code
___________________________________________________________________________________________________________________________________________________________
SOLVER STEP | LIN_IT LIN_RES RES_NORM RES_U RES_P RES_J NU_NORM NU_U NU_P NU_J
___________________________________________________________________________________________________________________________________________________________
- 0 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 516 2.432e-06 1.000e+00 1.000e+00 0.000e+00 0.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
- 1 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 290 4.660e-03 5.446e+00 8.554e-02 2.602e-10 1.437e+01 4.191e+00 2.550e+00 4.191e+00 1.456e+04
- 2 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 339 7.754e-03 2.780e+00 1.901e+00 7.310e-12 5.352e+00 3.597e+00 2.129e+00 3.597e+00 1.316e+03
- 3 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 325 2.931e-03 1.443e+00 1.380e+00 4.568e-12 1.115e+00 9.874e-01 4.266e-01 9.874e-01 1.250e+03
- 4 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 390 2.072e-04 8.894e-02 8.523e-02 2.469e-13 6.711e-02 4.998e-02 2.656e-02 4.998e-02 1.237e+02
- 5 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 506 8.187e-07 4.081e-04 3.921e-04 8.575e-16 2.980e-04 2.522e-04 1.344e-04 2.522e-04 5.493e-01
- 6 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 613 1.665e-11 3.105e-07 3.105e-07 2.613e-20 5.834e-09 3.108e-07 3.428e-08 3.108e-07 1.066e-05
- 7 ASM_R CONVERGED!
+ 0 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 786 2.118e-06 1.000e+00 1.000e+00 0.000e+00 0.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
+ 1 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 552 1.031e-03 8.563e-02 8.563e-02 9.200e-13 3.929e-08 1.060e-01 3.816e-02 1.060e-01 1.060e-01
+ 2 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 667 5.602e-06 2.482e-03 2.482e-03 3.373e-15 2.982e-10 2.936e-03 2.053e-04 2.936e-03 2.936e-03
+ 3 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 856 6.469e-10 2.129e-06 2.129e-06 2.245e-19 1.244e-13 1.887e-06 7.289e-07 1.887e-06 1.887e-06
+ 4 ASM_R CONVERGED!
___________________________________________________________________________________________________________________________________________________________
Relative errors:
-Displacement: 3.428e-08
-Force: 6.310e-12
+Displacement: 7.289e-07
+Force: 2.451e-10
Dilatation: 1.353e-07
v / V_0: 1.000e-09 / 1.000e-09 = 1.000e+00
___________________________________________________________________________________________________________________________________________________________
SOLVER STEP | LIN_IT LIN_RES RES_NORM RES_U RES_P RES_J NU_NORM NU_U NU_P NU_J
___________________________________________________________________________________________________________________________________________________________
- 0 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 556 2.542e-06 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
- 1 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 455 8.017e-05 1.538e-01 1.538e-01 3.522e+11 4.617e+07 6.010e-02 7.396e-02 6.010e-02 6.011e-02
- 2 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 554 4.012e-07 2.896e-03 2.896e-03 2.127e+09 5.117e+04 2.720e-03 1.430e-03 2.720e-03 2.724e-03
- 3 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 620 4.946e-10 1.728e-06 1.728e-06 1.315e+06 9.778e+01 1.150e-06 6.268e-07 1.150e-06 1.150e-06
+ 0 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 874 2.358e-06 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
+ 1 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 658 2.942e-04 1.544e-01 1.544e-01 1.208e+13 1.855e+06 6.014e-02 7.398e-02 6.014e-02 6.014e-02
+ 2 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 790 2.206e-06 2.908e-03 2.908e-03 7.302e+10 2.067e+03 2.716e-03 1.433e-03 2.716e-03 2.717e-03
+ 3 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 893 2.374e-09 1.919e-06 1.919e-06 4.527e+07 4.100e+00 1.672e-06 6.842e-07 1.672e-06 1.672e-06
4 ASM_R CONVERGED!
___________________________________________________________________________________________________________________________________________________________
Relative errors:
-Displacement: 6.268e-07
-Force: 1.873e-10
-Dilatation: 1.525e-06
-v / V_0: 9.999e-10 / 1.000e-09 = 9.999e-01
+Displacement: 6.842e-07
+Force: 8.995e-10
+Dilatation: 1.528e-06
+v / V_0: 1.000e-09 / 1.000e-09 = 1.000e+00
@endcode
<td align="center">
@image html "step-44.Normalised_runtime.png"
<p align="center">
- Runtime on a 4-core machine, normalised against the lowest grid resolution Q1-P0 solution that utilised a SSOR preconditioner.
+ Runtime on a 4-core machine, normalised against the lowest grid resolution $Q_1-DGPM_0-DGPM_0$ solution that utilised a SSOR preconditioner.
</p>
</td>
</tr>
// interchangeable. We choose to increment time linearly using a constant time
// step size.
//
-// We start the function with preprocessing, and then output the initial grid
-// before starting the simulation proper with the first time (and loading)
-// increment:
+// We start the function with preprocessing, setting the initial dilatation
+// values, and then output the initial grid before starting the simulation
+// proper with the first time (and loading)
+// increment.
+//
+// Care must be taken (or at least some thought given) when imposing the
+// constraint $\widetilde{J}=1$ on the initial solution field. The constraint
+// corresponds to the determinant of the deformation gradient in the undeformed
+// configuration, which is the identity tensor.
+// We use FE_DGPMonomial bases to interpolate the dilatation field, thus we can't
+// simply set the corresponding dof to unity as they correspond to the
+// monomial coefficients. Thus we use the VectorTools::project function to do
+// the work for us. The VectorTools::project function requires an argument
+// indicating the hanging node constraints. We have none in this program
+// So we have to create a constraint object. In its original state, constraint
+// objects are unsorted, and have to be sorted (using the ConstraintMatrix::close function)
+// before they can be used. Have a look at step-21 for more information.
+// We only need to enforce the initial condition on the dilatation.
+// In order to do this, we make use of a ComponentSelectFunction which acts
+// as a mask and sets the J_component of n_components to 1. This is exactly what
+// we want. Have a look at its usage in step-20 for more information.
template <int dim>
void Solid<dim>::run()
{
make_grid();
system_setup();
+ {
+ ConstraintMatrix constraints;
+ constraints.close();
+
+ const ComponentSelectFunction<dim>
+ J_mask (J_component, n_components);
+
+ VectorTools::project (dof_handler_ref,
+ constraints,
+ QGauss<dim>(degree+2),
+ J_mask,
+ solution_n);
+ }
output_results();
time.increment();
// this operation (we could, in principle simply create a new task using
// Threads::new_task for each cell) but there is not much harm done to doing
// it this way anyway.
+// Furthermore, should their be different material models associated with a
+// quadrature point, requiring varying levels of computational expense, then
+// the method used here could be advantageous.
template <int dim>
struct Solid<dim>::PerTaskData_UQPH
{
tangent_matrix.reinit(sparsity_pattern);
- // We then set up storage vectors noting
- // that the dilatation is unity
- // (i.e. $\widetilde{J} = 1$) in the
- // undeformed configuration...
+ // We then set up storage vectors
system_rhs.reinit(dofs_per_block);
system_rhs.collect_sizes();
solution_n.reinit(dofs_per_block);
solution_n.collect_sizes();
- solution_n.block(J_dof) = 1.0;
// ...and finally set up the quadrature
// point history: