3. Digital integration#

3.1. Discretized equations#

The time discretization of behavioral equations is based on an implicit Euler schema, i.e. the various variables of the problem are expressed at the final moment of the time step. The expression of stress at the end of the time step thus has the same relationship as that in the table, with deformation and damage also being expressed at the end of the time step. There is still the delicate part which consists in calculating the evolution of the damage.

Integrating the damage evolution equation is certainly a scalar but not linear problem. First of all, by following the resolution approach proposed in [Lorentz and Godard, 2011] and repeated in the fascicle [R5.04.01], a term affine in damage replaces the term in Laplacian for the damage field in the threshold function \(g\):

(3.1)#\[ \ stackrel {} {g} (a) =-A\ text {'} (a)\ Gamma (\ epsilon) +\ lambda +r (\ alpha -a) -k\]

At the stage of integration of the law of behavior, the quantities \(k\) and \(\alpha\) are considered known (as well as the positive penalty coefficient \(r\)), like the deformations at the end of the time step because they are expressed as a function of the nodal unknowns of the problem. Likewise, the value of the function \(\Gamma (\epsilon )\) is also known, even if its calculation also requires the resolution of a non-linear equation to which we will come back later. That’s why the above algorithmic threshold function only depends on \(a\). In the end, the problem to be solved is to find \(a\le 1\) such as:

(3.2)#\[ \ mathrm {\ Delta} a\ ge 0\ phantom {\ rule {3em} {0ex}};\ phantom {\ rule {3em} {0ex}}\ stackrel {~} {g} (a)\ le 0\ phantom {\ rule {3em} {0ex}};\ phantom {\ rule {3em} {0ex}}}\ mathrm {\ Delta} (a)\ le 0\ phantom {\ rule {3em} {0ex}} (a)\ le 0\ phantom {\ rule {3em} {0ex}} (a)\ le 0\ phantom {\ rule {3em} {0ex}} (a)\ le 0\ phantom {\ rule {3em} {0ex}} (a)\ le 0\ phantom {\ rule {3em} {\ stackrel {~} {g} (a) =0\]

where \(\Delta a=a–{a}_{n}\) refers to the damage increment during the time step. And if there is no solution (less than \(1\)), then \(a=1\) (saturated damage).

3.2. Solving the algorithmic evolution equation#

Since \(A\) is convex, \(\Gamma (\epsilon )\) is positive, and \(r\) is positive, the algorithmic threshold function is decreasing. There is therefore existence and uniqueness of the solution (important for the robustness of digital integration). In practice, the problem is solved in three steps:

  • If \(\stackrel{~}{g}({a}_{n})\le 0\) then \(a={a}_{n}\)

  • If \(\stackrel{~}{g}(1)\ge 0\) then \(a=1\)

  • Otherwise solve \(\stackrel{~}{g}(a)=0\) with \({a}_{n}\le a\le 1\) using a Newton method with controlled terminals

3.3. Effective calculation of the gamma function#

According to the table, the calculation of \(\Gamma (\epsilon )\) refers to that of \(\chi (\epsilon )\) which is only defined implicitly: it is therefore necessary to solve a non-linear equation that we know has a unique solution, cf. [Lorentz, 2016]. To do this, we start by introducing some notations and a change of variable in which \(x\) is the new unknown:

(3.3)#\[ {\ mathrm {\ sigma}} ^ {e} =\ mathrm {\ lambda}\ mathit {tr} (\ mathrm {\ epsilon})\ mathit {Id} +2\ mathrm {\ id} +2\ mathrm {\ epsilon} +2\ mathrm {\ epsilon}\]

The space for tensor \(s\) is generated by:

\[\]

: label: eq-4

mathit {span} (s) ={{s} _ {1}, {s} _ {2}, {s} _ {3}}

With:

\[\]

: label: eq-5

{s} _ {beta} =parallel s+ s+ ({beta} _ {0} -frac {1} {3})mathit {tr} (s)mathit {Id}parallel

The equation to be solved to determine \(x\) and therefore \({\rm X}\) is written as follows:

(3.4)#\[ {s} _ {\ beta} x+\ sqrt {\ sum _ {i=1} ^ {3}\ mathrm {exp} (2 {s} _ {i} x)}} - {\ gamma} _ {0} =0\]

The first limb is convex. We estimate a major \({x}_{0}\) of the solution then we solve the equation by a Newton method starting with \({x}_{0}\), which leads to a series of decreasing values and which converge towards the solution of the equation due to the predicted convexity. Note that the majority is obtained by one of the following two relationships, where \({s}_{M}\) designates the largest of the three eigenvalues of \(s\):

(3.5)#\[ {s} _ {\ beta} x+\ sqrt {\ sum _ {i=1}} ^ {3}\ mathrm {exp} (2 {s} _ {i} x)}\ ge {s} _ {\ beta} _ {\ beta} _ {\ beta} x\ Rightarrow x\ le\ frac {{\ gamma}} _ {\ beta}}\]

Or:

(3.6)#\[ {s} _ {\ mathrm {\ beta}} x+\ sqrt {\ sum _ {\ sum _ {i=1} ^ {3}\ mathrm {exp} (2 {s} _ {i} x)}\ ge\ mathrm {exp}}\ ge\ mathrm {exp}}\ ge\ mathrm {exp}}\ ge\ mathrm {exp}} ({s} _ {exp}}\ ge\ mathrm {exp}} ({s}}\ mathrm {exp}}\ ge\ mathrm {exp}} ({s})\ ge\ mathrm {exp}} ({s}}\ mathrm {exp}}\ ge\ mathrm {exp}} ({s} {2em} {0ex}} x\ le\ frac {\ mathrm {ln} {\ mathrm {\ gamma}} _ {0}} {{s}} {{s} _ {M}}\]

We will take the smaller of the two majors for the sake of performance (or the only one if \({s}_{M}\le 0\)).

3.4. Tangent matrix#

Determining the tangent matrix consistent with the chosen integration scheme does not raise any particular difficulty (derivation of compound functions and derivation of implicit functions), except that the eigenvalues of the deformation are involved in the calculation of \(\Gamma (\epsilon )\) and in that of the elastic law when distinguishing tension from compression. Here we examine each of these points more precisely, in order to avoid the derivative of the associated eigenvectors which, as we recall, are not defined when the eigenvalues are not distinct.

For the derivative of \(\Gamma (\epsilon )\), only the term that with \(x\) associates the norm of \(\mathrm{exp}(xs)\) requires special treatment. Gold:

(3.7)#\[ \ delta\ parallel\ mathrm {exp} (xs)\ parallel =\ frac {\ delta\ mathrm {exp} (xs)\ mathrm {:}\ mathrm {exp} (xs)} {\ parallel\ mathrm {xs) {exp} (xs)\ parallel} (xs)\ parallel\ mathrm {exp} (xs)\ parallel\ x+\ frac {\ mathrm {exp} (xs)} {\ parallel\ mathrm {exp} (xs)\ parallel}\ mathrm {:}\ mathrm {exp}\ text {'} (xs) (xs)\ mathrm {:} (x\ delta s)\]

At this stage, even if the derivative of the exponential (fourth-order tensor) appears formally in the expression, it is not useful to actually proceed with the derivative. Indeed, the derivative of the exponential is symmetric because it is the derivative of the following isotropic scalar function (cf. [Silhavy, 2000] for differentiability properties):

(3.8)#\[ f (S) =\ sum _ {i=1} ^ {3} {e} ^ {{A} _ {i}}\]

And so:

(3.9)#\[ \ mathrm {exp} (A)\ mathrm {:}\ mathrm {:}\ mathrm {exp}\ mathrm {exp}\ text {'} (A)\ mathrm {:}\ mathrm {:}\ mathrm {exp} (A) =\ mathrm {exp} (2A)\]

This last result is easily shown by noticing that \(A\) and \(\mathrm{exp}(A)\) switch. In the end, the calculation of the tangent matrix will involve that of the tensor exponential, according to the relationship above, without the need to derive the eigenvectors. And for the calculation of the tensor exponential, we will rely on the so-called « scaling and squaring » method, see for example [Higham, 2005], which does away with the calculation of eigenvectors.

As far as the stress-strain relationship is concerned, it shows a trace term that does not pose any particular derivation problem and an infinitely differentiable functional relationship between the natural stress \({\sigma }_{i}\) and the natural deformation \({\epsilon }_{i}\). For this one, we enter the framework of the derivative of tensor functions based on a scalar function of eigenvalues: the article [Carlson and Hoger, 1986] then provides the expressions of the (tensor) derivative whether the eigenvalues are distinct or not.