3. Digital integration#
With regard to the integration of the law of behavior, i.e. the calculation of internal variables and stresses with a given (mechanical) deformation history, there is no longer any need to distinguish whether it is a local formulation or a gradient formulation. In fact, the impact of this last one is limited to modifying the work hardening function according to (). Here we will continue to rate it \(R(\kappa )\) without distinction.
On the other hand, we will distinguish the plastic and viscoplastic cases, even if it can be shown that the latter is also limited to a correction of the threshold function in the discretized equations.
3.1. Discretized equations in the case without viscosity#
The time discretization of behavioral equations is based on a staged diagram: for fixed damage, an implicit Euler scheme is applied (even in viscoplasticity) to calculate the plastic variables (plastic deformation, work-hardening variable) as well as the stress, i.e. the different variables of the problem are expressed at the final moment of the time step in question. Then the evolution of the various contributions to damage is evaluated by knowing the evolution of the plastic variables. Refer to [Zhang et al., 2018] for more information on time discretization and the resolution of the corresponding algebraic system.
We note \({Q}_{n}\) the value of a quantity \(Q\) at the start of the time step, \(\mathrm{\Delta }Q\) its increment during the time step and (simply) \(Q\) its value at the end of the time step. The mechanical state at the start of time step \(({E}_{n},{{E}^{p}}_{n},{\kappa }_{n},{f}_{n}^{\mathit{nu}},{f}_{n}^{\mathit{gr}},{d}_{n}^{\mathit{co}},{d}_{n})\) is assumed to be known as well as the deformation increment \(\mathrm{\Delta }E\) (and therefore also the deformation \(E\)). It is then a question of calculating the increments of internal variables as well as the constraint at the end of time step \(T\).
Regarding the damage value chosen to integrate plasticity, [Zhang et al., 2018] proposed adopting the one at the beginning of the time step, i.e. \({d}_{\mathit{ex}}={d}_{n}\). For reasons of precision of integration (and therefore of the size of the time step), we enrich this diagram by extrapolating the damage during a time step by relying on its speed of evolution at the previous step. More precisely, we define:
where the value of \(\theta\) is set by the user under the COMPORTEMENT/PARM_THETA keyword; the value \(\theta =1/2\) is recommended.
Now knowing the extrapolated damage value, we can proceed with the implicit integration of plasticity. If \({d}_{\mathit{ex}}\ge {d}_{f}\), we simply switch to the broken point regime and \(T=0\). Otherwise, the discretization of the stress-strain relationship () is written as follows:
: label: eq-30
T=Kmathrm {tr} (E- {{E} ^ {e} ^ {p}}} _ {n} -mathrm {Delta} {E} ^ {p})mathit {Id} +2mathrm {mu}mathrm {mu}}mathrm {mu}}mathrm {mu}mathrm {mu}}mathrm {dev} (E- {{E}} ^ {n})mathit {Id} +2mathrm {mu}}mathrm {mu}}mathrm {mu}mathrm {mu}mathrm {mu}mathrm {dev} (E- {{E}} ^ {n}) -mathrm {Delta} {p})
We choose to show an elastic stress (or elastic test stress), noted \({T}^{e}\) and which is known:
: label: eq-31
T= {T} ^ {e} -Kmathrm {tr} (mathrm {Delta} {E} ^ {p}) -2mathrm {mu}mathrm {dev} (mathrm {Delta} {tr} (mathrm {Delta} {tr}) (mathrm {Delta} {tr} (mathrm {Delta} {tr} (mathrm {Delta} {tr} {E}) (mathrm {Delta} {tr} (mathrm {Delta} {tr} {E}) (mathrm {Delta} {tr} (mathrm {Delta} {tr}} (mathrm {Delta} {tr}) (mathrm {Delta}} ^ {p}} _ {n})mathit {Id} +2mathrm {mu}mathrm {dev} (E- {{E} ^ {p}}} _ {n})
The discretization of the threshold function is simply written:
And the consistency condition is equal to:
In particular, in plastic load, \(F=0\), so that it is possible to derive/to establish the equivalent stress of GTN as a function of the level of work hardening:
In the case of a regular flow, the temporal discretization of () is written as:
The following functions have been introduced:
By substituting () and () in (), we deduce:
This last expression indicates that \(T\) has the same direction and meaning as \({T}^{e}\). In addition, by substituting \(\mathrm{\Theta }\) with its expression from (), we get the value of \({T}_{\mathit{eq}}\) according to \({T}_{H}\) and \({T}_{\text{*}}\):
Given the definition of the equivalent constraint of GTN (), we can express an implicit link between \({T}_{H}\) and \({T}_{\text{*}}\) by substituting () for it:
The plastic multiplier, equal to \(\mathrm{\Delta }\kappa\), is deduced from ():
When the flow is singular, \({T}_{\text{*}}=0\) so \(T=0\), which makes it possible to calculate \(\mathrm{\Delta }{E}^{p}\) via (). It then remains to verify that the flow does indeed belong to the cone prescribed by () i.e., after discretization:
Knowing the evolution of plasticity \(\mathrm{\Delta }{E}^{p}\) and \(\mathrm{\Delta }\kappa\), we can now proceed with the integration of porosity. Germination porosity \({f}^{\mathit{nu}}\) is simply calculated using the functions () and (); the cumulative plastic deformation is evaluated as follows:
Concerning the growth porosity, a theta-scheme is applied, with the same value of \(\theta\) as during the prediction phase. The corresponding discretization is written as:
From this we deduce:
Finally, coalescence damage is simply defined via:
3.2. Solving the discretized problem without viscosity#
Preliminary note: the equations are solved to a certain degree of precision. If we write \(ϵ\) the precision provided by the user via the RESI_INTE_RELA keyword, we can enter the following details, respectively in terms of stress, deformation and without units for equation ():
If \(1-{d}_{\mathit{ex}}<{ϵ}^{1/2}\), we consider that the point is broken: \(T=0\), \(\mathrm{\Delta }\kappa =0\), \(f={f}_{F}\) (even if \({d}_{\mathit{ex}}<{d}_{f}\), because at the numerical precision requested, the « broken point » solution is actually acceptable) .Otherwise, the discretized problem must be solved.
The latter admits a unique solution. We start by determining which regime it falls under: elastic response, plastic response to singular flow or plastic response to regular flow. The solution is elastic if:
In this case, \(\mathrm{\Delta }\kappa =0\), \(\mathrm{\Delta }{E}^{p}=0\), and \(T={T}^{e}\).
If the solution is not elastic, the case of a singular plastic flow is examined. In this case, \({T}_{\text{*}}=0\), \(T=0\), \(\mathrm{\Delta }{E}^{p}\) is obtained via () and the corresponding value of \({\mathrm{\Pi }}_{N}\) is denoted by \({\mathrm{\Pi }}_{N}^{s}\) according to (). We also calculate \({\kappa }^{z}\) such that (resolution by a Newton method):
: label: eq-50
widehat {T} ({kappa} ^ {z}) =0
The flow is singular if \({\kappa }^{z}\ge {\kappa }_{n}\) and \(J\phantom{\rule{0.5em}{0ex}}{\mathrm{\Pi }}_{N}^{s}\le \mathrm{\Delta }{\kappa }^{z}\) according to (); we then have \(\kappa ={\kappa }^{z}\).
Finally, in the opposite case, the plastic flow is regular and in particular we have \(0<{T}_{\text{*}}<{T}_{\text{*}}^{e}\). It can be noted that at \({T}_{\text{*}}\) given, equation () makes it possible to go back to the value of \({T}_{H}\). By noting \({T}_{H}=p\phantom{\rule{0.5em}{0ex}}{T}_{H}^{e}\), the \(\widehat{g}\mathrm{:}p\to \widehat{G}({T}_{H},{T}_{\text{*}})\) function is increasing and the solution checks the \(0\le p\le 1\) framework. To facilitate the resolution, we can calculate a minimum bound \({p}_{\mathit{min}}>0\) by increasing the function \(\widehat{g}\). Noting \(r\equiv {T}_{H}/{T}_{\text{*}}\) and \(\alpha \equiv \mu /K\), we have:
We can analytically reduce the denominator of the second factor and we denote it \({D}_{\mathit{min}}\). As in addition :math:``, it comes:
: label: eq-52
widehat {g} (p)le {left (frac {{T}} _ {mathit {eq}}} ^ {e}} {{T} _ {text {*}}}}right)}right)}right)}}right)}} ^ {2} {2} {left} {frac {{T}}}right)} ^ {2} + {d} _ {mathit {ex}}}phantom {rule {0.5em} {0ex}} 2mathrm {sh} (frac {3} {2} {2} {q} {q} {q}} {q}} _ {2} r) - {(1- {d}} _ {mathit {ex}})} ^ {2}})} ^ {2} {2}
The root of the right-hand side can be obtained analytically and allows us to go back to the line: math: {p} _ {mathit {min}}}. A resolution of the equation:math: widehat {g} (p) =0 by a Newton method in the range:math:` {p} _ {mathit {min}}le ple 1` makes it possible to find the solution:math: p ({T} _ {text {*}}}) `. The precision is checked via the convergence condition:math: `left|widehat {g}right|le {} _ {G}.
We can thus reduce the incremental problem to the only unknown \(\mathrm{\Delta }\kappa\). Indeed, at a given \(\kappa\), we can calculate \({T}_{\text{*}}=\widehat{T}(\kappa )\). From there, we deduce \(p({T}_{\text{*}})\) by solving the previous equation. Knowing \(p\) therefore \({T}_{H}\), we can calculate \(\lambda =\widehat{\lambda }({T}_{H},{T}_{\text{*}})\) via (), which brings the problem back to solving the scalar equation \(\mathrm{\Delta }\kappa =\lambda\).
Thanks to the convexity properties of the equivalent stress \({T}_{\text{*}}\) compared to \(T\), the analysis of the « return mapping » scheme used shows that at \(\lambda\) increasing from 0 to \(J\phantom{\rule{0.5em}{0ex}}{\mathrm{\Pi }}_{N}^{s}\), \({T}_{\text{*}}\) is decreasing from \({T}_{\text{*}}^{e}\) to 0. Solving the problem therefore comes down to finding the intersection of the increasing curve \(\mathrm{\Delta }\kappa \to {T}_{\text{*}}\) with the decreasing curve \({T}_{\text{*}}\to \lambda\). A minimum bound on \(\mathrm{\Delta }\kappa\) is provided by the maximum of 0 and \(\mathrm{\Delta }{\kappa }_{z}\), because \({T}_{\text{*}}\) solution must be positive and \(\mathrm{\Delta }\kappa\) as well. A maximum bound is provided by the minimum of \(J\phantom{\rule{0.5em}{0ex}}{\mathrm{\Pi }}_{N}^{s}\) and \(\mathrm{\Delta }{\kappa }^{e}\) such as \(\widehat{T}({\kappa }^{e})={T}_{\text{*}}^{e}\) (if it exists), because \({T}_{\text{*}}\) solution must be positive and less than \({T}_{\text{*}}^{e}\).
The equation is solved by a Newton method between these two terminals. We check the precision of the solution by considering that the algorithm has converged as soon as one of the conditions:math: left|mathrm {Delta}kappa -lambdaright|le {} _ {kappa} _ {kappa} or:math: {kappa} or:math:left|left|widehat {T} (mathrm {Delta}right|le {}} _ {kappa} _ {kappa} or:math:left|left|widehat {T} (mathrm {Delta}right|le {} (lambda)right|the {}} _ {sigma} has been reached. We choose to initialize Newton’s method with the solution:math: `mathrm {Delta} {Delta} {kappa} _ {mathit {VM}} of the equivalent von Mises problem (i.e. by fixing:math: {T} _ {H} =0) when it respects the limits or otherwise by a string method with respect to the values at the terminals (i.e. by fixing: math:` {T} _ {H} =0`) when it respects the limits or otherwise by a string method with respect to the values at the terminals.
Finally, we go back to the constraint via the following relationship:
3.3. Solving the discretized problem in the viscoplastic case#
In the presence of viscosity, the consistency condition () is replaced by Norton’s law of evolution (). Once discretized, it is written:
And by reversing this relationship:
Like \(n>1\), \(0<q<1\).
It can still be expressed equivalently in the form of a Kuhn and Tucker relationship:
After discretization, the consideration of viscosity is reduced to the addition of a term (positive and increasing) in the work hardening function.
One might be tempted to apply the same resolution method as before but, unfortunately, this new term is not differentiable in \(\mathrm{\Delta }\kappa =0\) (infinite slope). This is why it is proposed to distinguish between two regimes, one with a steep slope, the other with a slight slope. The « border » slope is provided by \(R\text{'}(0)\) .More precisely, \(V(\mathrm{\Delta }\kappa )={C}_{\mathrm{\Delta }t}{\mathrm{\Delta }\kappa }^{q}\) is the term viscous. The border value \(\mathrm{\Delta }{\kappa }_{v}\) corresponds to:
If \(\mathrm{\Delta }\kappa <\mathrm{\Delta }{\kappa }_{v}\), we are in the regime with a steep slope, i.e. dominated by viscosity. Conversely, if \(\mathrm{\Delta }\kappa >\mathrm{\Delta }{\kappa }_{v}\), we are in the low-slope regime, dominated by plasticity. Given the monotony of the scalar equation \(\mathrm{\Delta }\kappa -\widehat{\lambda }=0\) to be solved with respect to \(\mathrm{\Delta }\kappa\), we can determine from the start the regime of the solution and reduce the search interval if necessary.
In the regime dominated by plasticity, the equation is solved in \(\mathrm{\Delta }\kappa\) as in the previous paragraph; the term viscous does not present any particular numerical difficulties.
On the contrary, in the regime dominated by viscosity, a Newton method involving \(\mathrm{\Delta }\kappa\) is not effective (too steep slopes therefore very slow convergence). We then change the variable \(v=V(\mathrm{\Delta }\kappa )\) (which eliminates the steep slopes in the vicinity of 0) and we apply the Newton method with respect to \(v\). The convergence criteria remain unchanged.
3.4. Tangent matrix#
In the local case, it is a question of determining the tensor of order 4 \({H}_{\mathit{TE}}=\delta T/\delta E\). In the non-local case, it is also necessary to build the 2nd order tensors \({H}_{T\varphi }=\delta T/\delta \varphi\) and \({H}_{\kappa E}=\delta \kappa /\delta E\) as well as the scalar \({H}_{\kappa \varphi }=\delta \kappa /\delta \varphi\).
The tangent operator depends on the flow regime. In the elastic case, it is reduced to:
In the case of the singular regime, the expression is also very simple:
where the function \(\widehat{T}\) contains the terms related to the non local and to the viscosity where appropriate.
The regular flow regime leads to a more complex tangent operator. Here, we simply describe the main lines of the method for calculating it. This regime is characterized by solving the equation system \(\widehat{G}({T}_{H},{T}_{\text{*}};{T}_{H}^{e},{T}_{\mathit{eq}}^{e},J)=0\) and \(\mathrm{\Delta }\kappa -\widehat{\lambda }({T}_{H},{T}_{\text{*}};{T}_{H}^{e},{T}_{\mathit{eq}}^{e},J)=0\) with the relationship \({T}_{\text{*}}=\widehat{T}(\mathrm{\Delta }\kappa ;J,\varphi )\). The derivative of the implicit functions makes it possible to express the variations of the unknowns \(\delta {T}_{H}\) and \(\delta \kappa\) as a function of the variations of the system’s input data (2x2 linear system to be solved). From there, we can go back to the variation in stress by differentiating ().
Finally, it remains to link the variations in the input data to those of deformation \(E\) and \(\varphi\). The only difficulty (minimal) comes from the von Mises constraint. To do this, we recall the following relationship:
3.5. Variables calculated during post-processing#
Four of the internal variables calculated during post-processing still need to be defined at this stage:
The part of the constraints linked to work hardening:
: label: eq-62
text {SIEQ_ ECR =} J (1- {d} _ {mathit {ex}}) R (mathrm {kappa})
The part of the constraints linked to viscosity:
The part of the constraints linked to the non local:
The extrapolation error expressed in constraints: