2. Benchmark solution#
In this test we manufacture the reference solution by estimating the energy return rate in two ways.
The first solution is obtained by the classical G-theta method (operator CALC_G), i.e. by estimating the surface integral on a crown placed around the crack.
In order to find this same value by another means, we first notice that in the case of linear elasticity, it is possible to estimate \(G\) using the Irwin formula linking the global quantities with the local quantities that are the stress intensity factors \({K}_{I};{K}_{\mathit{II}}\):
And:
Stress intensity factors can be estimated by analyzing the singularity of stress fields or displacement fields close to the peak. In the case of hyper-elasticity, the energy return rate is still an invariant of the integration contour, which allows us to imagine a similar procedure.
First we can theoretically contract the integration contour indefinitely towards the crack tip while keeping the same value of \(G\). So even for a hyper-elastic law we can establish a link between the global \(G\) and local quantities \(K\). However, the generalization of the concept of stress intensity factor is not always easy for all types of hyper-elastic laws. On the other hand, for a bilinear hyper-elastic law, this link can be established without too much difficulty. For this law, in the vicinity of the peak, the behavior will be linear regardless of the external load applied. We will therefore find the classical singularity of stress and deformation in the form of \(\sqrt{r}\).
In this test case we therefore use the non-linear elastic law of the Hencky type, which has the advantage of reproducing the von Mises plasticity law for any monotonic radial load:
where \(3K=E/(1-2\mathrm{\nu })\) is the compressibility module (in English: bulk modulus) and \(\mathrm{\mu }({\mathrm{\epsilon }}^{\mathit{eq}})\) is a Lamé coefficient, which is assumed to depend exclusively on the equivalent von Mises deformation: \({\mathrm{\epsilon }}^{\mathit{eq}}=\sqrt{\frac{3}{2}\cdot {\underline{\underline{\mathrm{\epsilon }}}}^{D}\mathrm{:}{\underline{\underline{\mathrm{\epsilon }}}}^{D}}\).
In the general form of Hencky’s nonlinear law, the dependence of the Lamé \(\mathrm{\mu }({\mathrm{\epsilon }}^{\mathit{eq}})\) coefficient on deformation is identified from the uniaxial tension curve. For the bilinear case the law is written as follows in the initial linear elasticity domain:
: label: eq-4
underline {underline {sigma}} (underline {underline {epsilon}}) =Kmathit {Tr} (underline {underline {epsilon}})underline {epsilon}})underline {epsilon}})underline {epsilon}})underline {epsilon}})underline {underline {epsilon}})underline {underline {epsilon}})underline {underline {epsilon}})underline {underline {epsilon}})underline {underline {epsilon}})
And in the field of non-linear elasticity:
: label: eq-5
underline {underline {sigma}} (underline {underline {epsilon}}) =Kmathit {Tr} (underline {underline {epsilon}})underline {underline {I}}} (underline {I}}}) = {sigma} _ {c} (1- {mu} _ {mathrm {infty}}})underline {underline {sigma}}}/mu)underline {underline {{epsilon} ^ {D}}}}/sqrt {{underline {epsilon}}} ^ {D}mathrm {:} {underline {underline {epsilon}}}}} ^ {D}}} +2 {mu} _ {mathrm {infty}} {infty}} {underline {infty}} {underline {infty}} {underline {epsilon}}}} ^ {D}} +2 {mu} _ {mathrm {infty}}}underline {infty}} {underline {infty}} {underline {epsilon}}}underline {epsilon Underline {{epsilon} ^ {D}}}
with \(3K=E/(1-2\mathrm{\nu })\) and \(2\mathrm{\mu }=E/(1+\mathrm{\nu })\).
The Lamé coefficient \(2{\mathrm{\mu }}_{\mathrm{\infty }}\) for loads exceeding the linear elasticity threshold is linked to the second slope of the law, which is noted \({E}_{\mathrm{\infty }}\) (keyword D_ SIGM_EPSI):
Finally, for large loads \({\mathrm{\epsilon }}^{\mathit{eq}}\gg \frac{{\mathrm{\sigma }}_{c}}{2\mathrm{\mu }}\), this law becomes asymptotically linear with the following linear characteristics:
With \(E\to {E}_{\mathrm{\infty }}\) and \(2{\mathrm{\nu }}_{\mathrm{\infty }}=1-(1-2\mathrm{\nu })\cdot {E}_{\mathrm{\infty }}/E\).
Analyzing the singularity of the displacement jump with these new values of elasticity coefficients with « infinite » loads allows us to identify the values of stress intensity factors using the operator already available in code_aster: POST_K1_K2_K3, then use the classical Irwin formula:, then use the classical Irwin formula:
Thus, two tests are carried out on the calculated energy return rates.
First test:
We test directly the ratio between the G from CALC_G and the G from POST_K1_K2_K3. This ratio should be equal to 1.
This is a reference of type ANALYTIQUE.
Second test:
The test case was run, for each modeling, on a mesh that was ten times finer at the bottom of the crack and then on a mesh that was a hundred times finer (size of the cells connected to the bottom of the crack respectively \({10}^{-4}\mathit{mm}\) and \({10}^{-5}\mathit{mm}\)), with a loading covering a wider range than that of the test case.
At each calculation moment, the values of \(G\) from the method 1 and 3 for calculating POST_K1_K2_K3 (see documentation U4.82.05) were extracted.
The comparison of the quantities calculated for each mesh shows the convergence (in the sense of the mesh) of these quantities.
The reference solution adopted is the value of \(G\) obtained with method 3 with the finest mesh (size of cells connected to the crack bottom of \({10}^{-5}\mathit{mm}\)). The relative difference between the value of \(G\) from methods 1 and 3 gives an estimate of the precision of this reference value.
This is a reference of type AUTRE_ASTER.