2. Formulation of the behavior model#
The use of the theory of thin plates and shells makes it possible to effectively describe the mechanical behavior of reinforced concrete structures, which are generally slender.
In a first step in the construction of the model, it is assumed the existence of a homogenized environment with the same mechanical behavior as the reinforced concrete structure, in which we are interested. To simplify, the hypothesis is made that this medium is isotropic and that the structural element studied is symmetric with respect to its middle sheet. These hypotheses are not essential for the formulation, but were made to simplify the process. In addition, it is estimated that their impact on behavior is lower compared to cracking, which is the focus of the model.
It should be noted that the use of this model is associated with the use of a plate element. It can only be used in the context of finite elements DKT (supported modeling: DKTG), corresponding to the**Love-Kirchhoff theory, where any transverse distortion in the thickness of the plate is neglected.
2.1. Free energy#
For an isotropic homogeneous continuous medium with linear elastic behavior we can write the volume density of free energy as:
\({\Phi }_{e}(\varepsilon )=\frac{\lambda }{2}\mathrm{tr}{(\varepsilon )}^{2}+\mu \sum _{i=1}^{3}{{\tilde{\varepsilon }}_{i}}^{2}\)
where \(\lambda ,\mu\) are the Lamé coefficients, \(\varepsilon\) the deformation tensor and \({\tilde{\varepsilon }}_{i}\) its eigenvalues. As for model ENDO_ISOT_BETON, damage is introduced by a function \(\xi (\mathrm{.},{d}_{i})\), \({d}_{i}\) being a damage variable. So, for a damaging medium, free energy is written as:
with the function \((x,d)\in {ℝ}^{2}\to \xi (x,d)\le 1\) verifying a priory \(\frac{\partial \xi }{\partial d}\le 0\) to represent the loss of stiffness linked to the damage, and \(\frac{\partial \xi }{\partial x}=0\mathrm{pour}x\in \left]-\infty \mathrm{,0}\right[\cup \left]\mathrm{0,}+\infty \right[\), the jump to \(0\) making it possible to represent the discontinuity in behavior between traction and compression.
The equation [eq] is valid for a continuous medium and we will apply it to a Love-Kirchhoff plate \(\omega \times \left]\frac{-h}{2};\frac{h}{2}\right[\), of thickness \(h\) (we note \(z={x}_{3}\)), where we make Hencky-Mindlin kinematic hypotheses (see [bib9]):
\((\begin{array}{}{U}_{1}({x}_{\mathrm{1,}}{x}_{\mathrm{2,}}z)\\ {U}_{2}({x}_{\mathrm{1,}}{x}_{\mathrm{2,}}z)\\ {U}_{z}({x}_{\mathrm{1,}}{x}_{\mathrm{2,}}z)\end{array})=\underset{\begin{array}{}\mathrm{cinématique}\mathrm{de}\mathrm{plaque}\\ {u}^{s}\in {V}_{s}\end{array}}{\underset{\underbrace{}}{(\begin{array}{}{u}_{1}({x}_{\mathrm{1,}}{x}_{2})\\ {u}_{2}({x}_{\mathrm{1,}}{x}_{2})\\ {u}_{z}({x}_{\mathrm{1,}}{x}_{2})\end{array})+z(\begin{array}{}{\theta }_{2}({x}_{\mathrm{1,}}{x}_{2})\\ -{\theta }_{1}({x}_{\mathrm{1,}}{x}_{2})\\ 0\end{array})}}+\underset{\begin{array}{}\mathrm{déplacement}\mathrm{complémentaire}\\ {u}^{c}\in {V}_{c}\end{array}}{\underset{\underbrace{}}{(\begin{array}{}{u}_{1}^{c}({x}_{\mathrm{1,}}{x}_{\mathrm{2,}}z)\\ {u}_{2}^{c}({x}_{\mathrm{1,}}{x}_{\mathrm{2,}}z)\\ {u}_{z}^{c}({x}_{\mathrm{1,}}{x}_{\mathrm{2,}}z)\end{array})}}\)
where \(U={({U}_{1}{U}_{2}{U}_{z})}^{T}\) is the 3D displacement field, \(u={({u}_{1}{u}_{2}{u}_{z})}^{T}\) the displacement of the middle sheet and \({\theta }_{x}\), \({\theta }_{y}\) its rotations. Thus, the strain tensor, defined as:
\({\varepsilon }_{\mathrm{ij}}=\frac{1}{2}(\frac{\partial {U}_{i}}{\partial {x}_{j}}+\frac{\partial {U}_{j}}{\partial {x}_{i}}),i,j=\mathrm{1..3}\)
is also written like:
where \(\epsilon\) is the membrane extension tensor, defined in the plane:
\({\epsilon }_{\mathrm{ij}}=\frac{1}{2}(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}),i,j=\mathrm{1..2}\)
and \(\kappa\) the curvature variation tensor, defined in the plane:
\({\kappa }_{11}=\frac{\partial {\theta }_{2}}{\partial {x}_{1}}\), \({\kappa }_{22}=\frac{-\partial {\theta }_{1}}{\partial {x}_{2}}\), \({\kappa }_{12}=\frac{1}{2}(\frac{\partial {\theta }_{2}}{\partial {x}_{2}}-\frac{\partial {\theta }_{1}}{\partial {x}_{1}})\)
relationships to which is added the hypothesis of plane constraints \({\sigma }_{\mathrm{zz}}=0\), \({\sigma }_{\mathrm{1z}}=0\), \({\sigma }_{\mathrm{2z}}=0\) which will determine the complementary displacement field \({u}^{c}\in {V}_{c}\). In the theory used here, only two rotation components \({\theta }_{x}\) and \({\theta }_{y}\) are introduced, which implies that the curvature variation tensor is 2D and has only 3 independent components.
By introducing these kinematic hypotheses, cf. [eq.], into the expression for free energy, [eq.], we can determine the eigenvalues of deformation \({\tilde{\varepsilon }}_{i}=\tilde{{(\epsilon +z\kappa )}_{i}}\). Since these eigenvalues are, in general, non-polynomial functions of the coordinate \(z\), the integral \(\int {\phi }_{\mathrm{ed}}\mathrm{dz}\) cannot be calculated analytically. This formulation would therefore not be suitable for application to structural elements.
Instead of the formulation [eq.], we will rather use a free energy formulation defined directly in generalized deformations, \(\epsilon\) and \(\kappa\):
where \(z(\mathrm{.})\) contains all the terms coupling \(\epsilon\) and \(\kappa\), which disappear after integration on \(z\), if we assume that the plate/beam is symmetric with respect to the middle sheet. The surface density of free energy is thus obtained:
where \(h\) is the thickness of the plate and \({\lambda }_{m},{\lambda }_{f},{\mu }_{m},{\mu }_{f}\) is defined in section § 3.2. The surface density of free energy is strictly convex for \({d}_{j}=0\) (having correctly chosen the elastic coefficients).
In [eq] and [eq] we also made the assumption that the damage is not affected by the extension in \(z\), which results in the absence of \({\varepsilon }_{\mathrm{zz}}\) from the arguments of the damage indicator \({\xi }_{m}\). This is justified by the objective of this model, which is to represent cracking perpendicular to the middle sheet, which is triggered either by membrane stress or by flexure, but never in extension by \({\varepsilon }_{\mathrm{zz}}\). Moreover, from the point of view of numerical resolution discussed below, this hypothesis facilitates the local calculation of damage variables \({d}_{j}\) and the satisfaction of the plane stress hypothesis \({\sigma }_{\mathrm{zz}}=0\).
Note:
In [eq. ], we observe that the deformation \({\varepsilon }_{\mathrm{zz}}\) is only explicitly introduced on the membrane energy term. The effect of plane stresses on membrane behavior will therefore be affected by the damage. This is not the case in flexure. One could imagine that the plane stress condition could be written by integrating a flexion-membrane coupling, which in turn depends on the damage, but this was not chosen here.
We will see section § 3.2 how to determine the parameters of the model from its answers on simple cases.
As for the damage variable, it is described by two components, one of the two components representing overall the damage on the side of the upper face of the plate and the other representing the damage on the side of the lower face of the plate:
\(d(z)=\{\begin{array}{ccc}{d}_{1}& \text{si}& z\ge 0\\ {d}_{2}& \text{si}& z<0\end{array}\)
It remains to define the characteristic functions of damage \({\xi }_{m}(\mathrm{.},{d}_{j})\) and \({\xi }_{f}(\mathrm{.},{d}_{j})\), for the formulation, [eq], to be complete:
\({\xi }_{m}(x,{d}_{\mathrm{1,}}{d}_{2})=\frac{1}{2}((\frac{1+{\gamma }_{\mathrm{mt}}{d}_{1}}{1+{d}_{1}}+\frac{1+{\gamma }_{\mathrm{mt}}{d}_{2}}{1+{d}_{2}})H(x)+(\frac{{\alpha }_{c}+{\gamma }_{\mathrm{mc}}{d}_{1}}{{\alpha }_{c}+{d}_{1}}+\frac{{\alpha }_{c}+{\gamma }_{\mathrm{mc}}{d}_{2}}{{\alpha }_{c}+{d}_{2}})H(-x))\in [\mathrm{0,1}]\)
and
\({\xi }_{f}(x,{d}_{\mathrm{1,}}{d}_{2})=\frac{\alpha +{\gamma }_{f}{d}_{1}}{\alpha +{d}_{1}}H(x)+\frac{\alpha +{\gamma }_{f}{d}_{2}}{\alpha +{d}_{2}}H(-x)\in [\mathrm{0,1}]\)
where \(H(\mathrm{.})\) is the Heaviside function.
The functions \((x,d)\in {ℝ}^{2}\to \xi (x,d)\le 1\), which are expressed in the same way as that chosen in law ENDO_ISOT_BETON, cf. [R7.01.04], by changing the variable on \({d}_{i}\), offer the advantage of giving a constant slope in the phases of damage. They are decreasing \(\frac{\partial \xi }{\partial d}\le 0\) and convex to represent the loss of stiffness due to damage, and \(\frac{\partial \xi }{\partial x}=0\mathrm{pour}x\in \left]-\infty \mathrm{,0}\right[\cup \left]\mathrm{0,}+\infty \right[\), the jump to \(0\) allowing to represent the change in behavior between traction and compression, however without introducing discontinuity. Damage variables \({d}_{i}\) increase to \(+\infty\).
The characteristic functions of damage \({\xi }_{m}(\mathrm{.},{d}_{j})\) and \({\xi }_{f}(\mathrm{.},{d}_{j})\) vary from \(1\) to \({\gamma }_{\mathrm{mt}}\) or \({\gamma }_{\mathrm{mc}}\) respectively, and \({\gamma }_{f}\), for \({d}_{j}\to +\infty\). This model therefore describes a partial but not total step damage to the material.
The damage parameters \({\gamma }_{\mathrm{mt}}\) for membrane traction, \({\gamma }_{\mathrm{mc}}\) for membrane compression, and for membrane compression, and \({\gamma }_{f}\) for flexure, may have values in \([\mathrm{0,1}]\), so that the model is not softening, which would cause difficulties in dependence on spatial discretization and convergence. We will choose \({\gamma }_{\mathrm{.}}\approx 0\) when the corresponding phenomenon will have more impact on the damage and \({\gamma }_{\mathrm{.}}\approx 1\) when this one is negligible. So, for reinforced concrete, we expect \({\gamma }_{\mathrm{mc}}\approx 1\) and \({\gamma }_{\mathrm{mt}}\approx 0\). As for the parameter, \(\alpha\), it allows you to adjust the contribution of flexure to the damage threshold. The parameter \({\alpha }_{c}\) makes it possible to modulate the evolution of compression damage.
2.2. Generalized constraints#
According to the usual procedure, generalized stresses (normal forces and moments) are defined by the derivatives of the free energy density with respect to the generalized deformations, \(\epsilon\) and \(\kappa\):
\(N=\frac{\partial {\Phi }_{\mathrm{ed}}^{S}}{\partial \epsilon }\); \(M=\frac{\partial {\Phi }_{\mathrm{ed}}^{S}}{\partial \kappa }\)
In our application, the generalized stresses are calculated in the coordinate system of the eigenvectors of the generalized deformations and they are written as (see the calculation in the annex):
It is easy to verify (see annex) that the generalized constraints \(N\) and \(M\), defined from the expressions [eq], respectively share the same eigenvectors as the generalized deformations \(\epsilon\) and \(\kappa\). If these eigenvectors are designated as \({Q}_{m}\) and \({Q}_{f}\) respectively, we can write:
\(N={Q}_{m}\tilde{N}{Q}_{m}^{T}\) and \(M={Q}_{f}\tilde{M}{Q}_{f}^{T}\) |
(2.2-2) |
where \(\tilde{N}\) and \(\tilde{M}\) are the diagonal matrices composed of the eigenvalues defined in [eq]. It is important to note that the eigenvectors for the membrane portion \({Q}_{m}\) and the flexure portion and \({Q}_{f}\) are completely independent.
Similarly, pinch constraints are defined as the dual variable in \({\varepsilon }_{\mathrm{zz}}\):
to which we will impose the condition of plane constraints: \({\sigma }_{\mathrm{zz}}=0\).
2.3. Thresholds and evolution of damage#
To be able to define a damage threshold within the framework of the hypothesis of a generalized standard material (see [bib 1], [bib 7]), the thermodynamic forces associated with the variables \({d}_{1}\) and \({d}_{2}\) are introduced:
where
\({Y}_{j}^{m}(\epsilon ,{d}_{j},{\varepsilon }_{\mathrm{zz}})=\frac{1}{{(1+{d}_{j})}^{2}}(\frac{{\lambda }_{m}}{4}{(\mathrm{tr}(\epsilon )+{\varepsilon }_{\mathrm{zz}})}^{2}\mathrm{.}{G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+\frac{{\mu }_{m}}{2}(\sum _{i=1}^{2}{\tilde{\epsilon }}_{i}^{2}\mathrm{.}{G}_{m}({\tilde{\epsilon }}_{i},{d}_{j})))\) by noting \({G}_{m}(x,{d}_{j})=(1-{\gamma }_{\mathrm{mt}})H(x)+(\frac{{\alpha }_{c}(1-{\gamma }_{\mathrm{mc}}){(1+{d}_{j})}^{2}}{{({\alpha }_{c}+{d}_{j})}^{2}})H(-x)\in [\mathrm{0,1}]\) so that: \(\frac{\partial {\xi }_{m}(x,{d}_{\mathrm{1,}}{d}_{2})}{\partial {d}_{j}}=-\frac{{G}_{m}(x,{d}_{j})}{2{(1+{d}_{j})}^{2}}\), and \({Y}_{j}^{f}(\kappa ,{d}_{j})=\frac{\alpha }{{(\alpha +{d}_{j})}^{2}}{Y}_{j}^{f\mathrm{,0}}(\kappa )\) with \({Y}_{j}^{f\mathrm{,0}}(\kappa )=(1-{\gamma }_{f})(\frac{{\lambda }_{f}}{2}\mathrm{tr}{(\kappa )}^{2}H({(-1)}^{(j+1)}\mathrm{.}\mathrm{tr}(\kappa ))+{\mu }_{f}\sum _{i=1}^{2}{\tilde{\kappa }}_{i}^{2}H({(-1)}^{(j+1)}\mathrm{.}{\tilde{\kappa }}_{i}))\) |
(2.3-2) |
The damage thresholds are defined by:
where \({k}_{{0}_{j}}\) are threshold constants. These thresholds define the convex domain of reversibility in the \((\epsilon ,\kappa )\) space.
In principle the threshold constants \({k}_{{0}_{1}}\) and \({k}_{{0}_{2}}\) could be different, but depending on the hypothesis we make about the symmetry of the plate with respect to the middle sheet, the two values are the same: \({k}_{{0}_{1}}={k}_{{0}_{2}}={k}_{0}\) .
We can see in [eq] that the parameter \(\alpha\) regulates the contribution of flexure to the initial damage threshold, since:
\({f}_{{d}_{j}}({d}_{j}=0)={Y}_{j}^{m}(\epsilon ,{d}_{j}=0)+\frac{1}{\alpha }{Y}_{j}^{f\mathrm{,0}}(\kappa )-{k}_{0}\)
The law of evolution of the damage variables \({d}_{1}\) and \({d}_{2}\) is defined by the rule of normality at the thresholds [eq], for which we can define the pseudo-dissipation potential \(D(\delta )\):
The damage values \({d}_{1}\) and \({d}_{2}\) are perfectly determined by the following conditions:
The evolution of the damage variables is therefore obtained using the coherence condition, the functions \(\xi (x,d)\) being convex, the work-hardening modules \(-{f}_{,{d}_{j}}\) are positive (the coefficients verifying \({\gamma }_{\mathrm{.}}\in [\mathrm{0,1}]\)):
It can be seen that in a damaging pure uniaxial membrane load (\(\dot{{d}_{1}}=\dot{{d}_{2}}\ge 0\)), the thermodynamic force is expressed: \({Y}^{m}(d)=-{E}_{\mathrm{eq}}^{m}h\frac{{\varepsilon }^{2}\mathrm{.}{\xi }_{m,d}(d)}{2}\), \({E}_{\mathrm{eq}}^{m}h\) being the uniaxial membrane elastic stiffness. The consistency condition is then written as:
Hence: \(\dot{{d}_{m}}=-\frac{2\dot{\epsilon }{\xi }_{m,d}}{\epsilon {\xi }_{m,\mathrm{dd}}}\) and therefore the law of pure uniaxial damaging membrane charge is:
This makes it possible to interpret the role of the parameter \({\gamma }_{m}\): the slope being constant, which is the justification for the algebraic form of the function \({\xi }_{m}\).
2.4. Digital integration#
Unlike the majority of nonlinear anelastic laws of behavior, the one presented here does not require the discretization of the equations for the evolution of internal variables. A direct discretization method is adopted that is implicit in time.
Variables \({d}_{j}\) can be calculated directly from the consistency condition, [eq]. The only time we refer to the « damage rate » is to verify that it is positive \({\dot{d}}_{j}\ge 0\). For an incremental calculation this condition results in \({d}_{j}^{n}\ge {d}_{j}^{n-1}\) at the \(n>1\) time step, therefore without reference to a particular temporal integration scheme.
Let’s put \({t}_{n}\) on the loading path at a given point in time. First, an elastic prediction step is carried out (elasticity tensor evaluated with the damage variables \({d}_{j}^{n-1}\) fixed in the previous step), hence \(({\epsilon }^{n},{\kappa }^{n})\). We then calculate a new \({\varepsilon }_{\mathrm{zz}}^{\mathrm{n0}}\):
\({\sigma }_{\mathrm{zz}}^{n}=0\Rightarrow {\varepsilon }_{\mathrm{zz}}^{\mathrm{n0}}({\epsilon }^{n},{d}_{1}^{n-1},{d}_{2}^{n-1})=-\frac{{\lambda }_{m}{\xi }_{m}(\mathrm{tr}({\epsilon }^{n}),{d}_{1}^{n-1},{d}_{2}^{n-1})}{2{\mu }_{m}+{\lambda }_{m}{\xi }_{m}(\mathrm{tr}({\epsilon }^{n}),{d}_{1}^{n-1},{d}_{2}^{n-1})}\mathrm{tr}({\epsilon }^{n})\)
see [eq].
We then calculate \({f}_{{d}_{j}}^{\mathrm{n0}}({\epsilon }^{n},{\kappa }^{n},{d}_{j}^{n-1},{\varepsilon }_{\mathrm{zz}}^{\mathrm{n0}})\) to check the damage thresholds. If \({f}_{{d}_{j}}^{\mathrm{n0}}\le 0\), the damage does not evolve: \({d}_{j}^{n}={d}_{j}^{n-1}\) and \({\varepsilon }_{\mathrm{zz}}^{n}={\varepsilon }_{\mathrm{zz}}^{\mathrm{n0}}\), and the generalized constraints are calculated according to [eq].
When \({f}_{{d}_{j}}^{\mathrm{n0}}>0\), the damage can evolve and we must solve the equations:
\({f}_{{d}_{j}}^{n}={Y}_{j}({\epsilon }^{n},{\kappa }^{n},{d}_{j}^{n},{\varepsilon }_{\mathrm{zz}}^{n})-{k}_{0}=0\)
Which corresponds to the solution of the nonlinear equations, to \(\epsilon ,\kappa\) given:
The equation [eq] must be solved by also taking into account the plane stress condition, which allows \({\varepsilon }_{\mathrm{zz}}(\epsilon ,{d}_{j})\) to be expressed, cf. [eq]:
The equations [eq] and [eq] are solved by Newton’s method. We start with a prediction phase, of the explicit Euler type:
\({d}_{j}^{n(0)}={d}_{j}^{n-1}+{(\frac{{\mathrm{dR}}_{{d}_{j}}}{d({d}_{j})})\mid }_{{d}_{j}^{n-1}}^{-1}\mathrm{.}({k}_{0}-{Y}_{j}({\epsilon }^{n},{\kappa }^{n},{d}_{j}^{n-1},{\varepsilon }_{\mathrm{zz}}^{\mathrm{n0}}))\)
Then we treat the correction phase, in iteration \((m)>1\):
\(\Delta {d}_{j}^{n(m)}={(\frac{{\mathrm{dR}}_{{d}_{j}}}{d({d}_{j})})\mid }_{{d}_{j}^{n(m-1)}}^{-1}\mathrm{.}({k}_{0}-{Y}_{j}({\epsilon }^{n},{\kappa }^{n},{d}_{j}^{n(m-1)},{\varepsilon }_{\mathrm{zz}}^{n(m-1)}))\)
\({d}_{j}^{n(m)}={d}_{j}^{n(m-1)}+\Delta {d}_{j}^{n(m)}\)
\({\varepsilon }_{\mathrm{zz}}^{n(m)}={\varepsilon }_{\mathrm{zz}}({\epsilon }^{n},{d}_{1}^{n(m)},{d}_{2}^{n(m)})\)
This correction phase is completed when the convergence criterion expressed in energy terms is reached:
\(\Delta {d}_{j}^{n(m)}\mathrm{.}{R}_{{d}_{j}}^{n(m)}<{\eta }_{\mathrm{tolerance}}\mathrm{.}{k}_{0}\)
The tangent operator for this nonlinear system is defined as:
where:
and:
We check that we have \(\frac{\partial {\xi }_{m}}{\partial {d}_{j}}(x,{d}_{\mathrm{1,}}{d}_{2})<0\) as expected.
2.5. Tangent stiffness operator#
As the main objective of a global model is to offer a simplified approach to the modeling of a complex material, such as reinforced concrete, it is essential that its numerical performance is optimal. Thus, to make the model suitable for calculations with schemes that are implicit in time, either in quasistatics or in transient dynamics, the calculation of the coherent tangent stiffness becomes essential in order to have a quadratic and robust convergence of the global iterative Newton process.
Most of the calculation of the model is carried out in the coordinate system of the eigenvectors of the generalized deformation tensors (and of the generalized constraints, cf. [eq]), the tangent stiffness therefore also being expressed in the same coordinate system. The transformation required to be able to then use it in the assembly of the overall stiffness matrix is specified in the next chapter.
In order to simplify the writing, we define a vector of generalized constraints (membrane force, bending moment), \(\Sigma\), and a vector of generalized deformations (extension, curvature), \(E\), as:
\(\Sigma ={({\tilde{N}}_{1}{\tilde{N}}_{2}{\tilde{M}}_{1}{\tilde{M}}_{2})}^{T}\)
\(E={({\tilde{\epsilon }}_{1}{\tilde{\epsilon }}_{2}{\tilde{\kappa }}_{1}{\tilde{\kappa }}_{2})}^{T}\)
The tangent stiffness operator \(C\) is defined by the relationship in the real evolution:
\(d\Sigma =C\cdot dE\)
It can be calculated as the sum of two contributions, the one that corresponds to a non-evolution of the damage and the one that is due to the evolution of the damage. These contributions can be called: elastic contribution and dissipative contribution:
In addition, we take into account the structure of the \(C\) tensor composed of the contributions normal effort-extension, moment-curvature and their couplings. More particularly, the tensors \({C}_{e}\) and \({C}_{d}\) are of the following form:
We can see from [eq] that the coupling moments/extension and membrane/curvature forces is only introduced through the dissipative part. This coupling has a physical justification, since any cracking perpendicular to the middle sheet of the plate affects both membrane and flexural behavior.
The sub-matrices \({C}_{e}^{\mathrm{mm}}\), \({C}_{e}^{\mathrm{ff}}\), \({C}_{d}^{\mathrm{mm}}\), \({C}_{d}^{\mathrm{mf}}\), \({C}_{d}^{\mathrm{ff}}\) are given in the proper coordinate system by the following expressions:
where
\(\frac{\partial {\tilde{N}}_{i}}{\partial {d}_{k}}={\lambda }_{m}(\mathrm{tr}(\epsilon )+{\varepsilon }_{\mathrm{zz}})\frac{\partial {\xi }_{m}(\mathrm{tr}(\epsilon ),{d}_{j})}{\partial {d}_{k}}+{\lambda }_{m}{\xi }_{m}(\mathrm{tr}(\epsilon ),{d}_{j})\frac{\partial {\varepsilon }_{\mathrm{zz}}}{\partial {d}_{k}}+2{\mu }_{m}{\tilde{\epsilon }}_{i}\frac{\partial {\xi }_{m}({\tilde{\epsilon }}_{i},{d}_{j})}{\partial {d}_{k}}\) and \(\frac{\partial {\tilde{M}}_{i}}{\partial {d}_{k}}={\lambda }_{f}\mathrm{tr}(\kappa )\frac{\partial {\xi }_{f}(\mathrm{tr}(\kappa ),{d}_{\mathrm{1,}}{d}_{2})}{\partial {d}_{k}}+2{\mu }_{f}{\tilde{\kappa }}_{i}\frac{\partial {\xi }_{f}({\tilde{\kappa }}_{i},{d}_{\mathrm{1,}}{d}_{2})}{\partial {d}_{k}}\) with \(\frac{\partial {\xi }_{f}}{\partial {d}_{k}}(x,{d}_{\mathrm{1,}}{d}_{2})=-\alpha \frac{(1-{\gamma }_{f})}{{(\alpha +{d}_{k})}^{2}}H({(-1)}^{(k+1)}x)\) |
(2.5-5) |
In addition to the expressions [eq], we use the equations [eq] to find \(\frac{\partial {\xi }_{m}}{\partial {d}_{k}}(x,{d}_{\mathrm{1,}}{d}_{2})\).
The derivatives \({\frac{d({d}_{k})}{d{\tilde{\epsilon }}_{i}}\mid }_{{f}_{{d}_{k}}=0}\) and \({\frac{d({d}_{k})}{d{\tilde{\kappa }}_{i}}\mid }_{{f}_{{d}_{k}}=0}\) are calculated by differentiating equation \({R}_{{d}_{k}}=0\) respectively with respect to \({d}_{k}\), \({\tilde{\epsilon }}_{i}\) and \({\tilde{\kappa }}_{i}\), (see [eq.]). When both damage mechanisms are activated, the following systems need to be resolved:
with \({A}_{\mathrm{jk}}=\frac{\partial {Y}_{j}}{\partial {d}_{k}}\), \({B}_{k}^{m}=\frac{-\partial {Y}_{k}}{\partial \epsilon }=\frac{\partial \tilde{N}}{\partial {d}_{k}}\), \({B}_{k}^{f}=\frac{-\partial {Y}_{k}}{\partial \kappa }=\frac{\partial \tilde{M}}{\partial {d}_{k}}\);
So:
\(\frac{d({d}_{1})}{d{\tilde{\epsilon }}_{i}}=\frac{{A}_{22}{B}_{1}^{m}-{A}_{12}{B}_{2}^{m}}{{A}_{22}{A}_{11}-{A}_{12}{A}_{21}}\) and \(\frac{d({d}_{2})}{d{\tilde{\epsilon }}_{i}}=\frac{{A}_{11}{B}_{2}^{m}-{A}_{21}{B}_{1}^{m}}{{A}_{11}{A}_{22}-{A}_{21}{A}_{12}}\)
\(\frac{d({d}_{1})}{d{\tilde{\kappa }}_{i}}=\frac{{A}_{22}{B}_{1}^{f}-{A}_{12}{B}_{2}^{f}}{{A}_{22}{A}_{11}-{A}_{12}{A}_{21}}\) and \(\frac{d({d}_{2})}{d{\tilde{\kappa }}_{i}}=\frac{{A}_{11}{B}_{2}^{f}-{A}_{21}{B}_{1}^{f}}{{A}_{11}{A}_{22}-{A}_{21}{A}_{12}}\)
where:
with:
and \(\frac{\partial {\varepsilon }_{\mathrm{zz}}}{\partial {d}_{k}}\) given by [eq.]. The matrix \(({A}_{\mathrm{jk}})\) is indeed invertible, see section § 2.4.
2.6. Change of frame of reference#
The coordinate system change approach is identical to that developed for model ENDO_ISOT_BETON (see section [2.4.4.1] of [R7.01.04]) with the only difference that it applies to generalized stresses and deformations. We thus obtain the components of:
\(C=(\begin{array}{cc}{C}^{\mathrm{mm}}& {C}^{\mathrm{mf}}\\ {({C}^{\mathrm{mf}})}^{T}& {C}^{\mathrm{ff}}\end{array})\)
as being:
\(\begin{array}{ccccc}{({C}^{\mathrm{mm}})}_{\mathrm{ijkl}}& =& \frac{\partial {N}_{\mathrm{ij}}}{\partial {\epsilon }_{\mathrm{kl}}}& =& \sum _{m,n}{Q}_{\text{im}}^{m}{Q}_{\mathrm{jm}}^{m}{Q}_{\mathrm{kn}}^{m}{Q}_{\text{ln}}^{m}\cdot \frac{\partial {\tilde{N}}_{m}}{\partial {\tilde{\epsilon }}_{n}}\\ & & & +& \frac{1}{2}\sum _{\begin{array}{}m,n\\ m\ne n\end{array}}(\frac{({Q}_{\mathrm{km}}^{m}{Q}_{\text{ln}}^{m}+{Q}_{\mathrm{lm}}^{m}{Q}_{\mathrm{kn}}^{m})({Q}_{\text{in}}^{m}{Q}_{\mathrm{jm}}^{m}+{Q}_{\mathrm{jn}}^{m}{Q}_{\text{im}}^{m})}{{\tilde{\epsilon }}_{n}-{\tilde{\epsilon }}_{m}}){\tilde{N}}_{m}\end{array}\)
\(\begin{array}{ccccc}{({C}^{\mathrm{ff}})}_{\mathrm{ijkl}}& =& \frac{\partial {M}_{\mathrm{ij}}}{\partial {\kappa }_{\mathrm{kl}}}& =& \sum _{m,n}{Q}_{\text{im}}^{f}{Q}_{\mathrm{jm}}^{f}{Q}_{\mathrm{kn}}^{f}{Q}_{\text{ln}}^{f}\cdot \frac{\partial {\tilde{M}}_{m}}{\partial {\tilde{\kappa }}_{n}}\\ & & & +& \frac{1}{2}\sum _{\begin{array}{}m,n\\ m\ne n\end{array}}(\frac{({Q}_{\mathrm{km}}^{f}{Q}_{\text{ln}}^{f}+{Q}_{\mathrm{lm}}^{f}{Q}_{\mathrm{kn}}^{f})({Q}_{\text{in}}^{f}{Q}_{\mathrm{jm}}^{f}+{Q}_{\mathrm{jn}}^{f}{Q}_{\text{im}}^{f})}{{\tilde{\kappa }}_{n}-{\tilde{\kappa }}_{m}}){\tilde{M}}_{m}\end{array}\)
\({({C}^{\mathrm{mf}})}_{\mathrm{ijkl}}=\frac{\partial {N}_{\mathrm{ij}}}{\partial {\kappa }_{\mathrm{kl}}}=\sum _{m,n}{Q}_{\text{im}}^{m}{Q}_{\mathrm{jm}}^{m}{Q}_{\mathrm{kn}}^{f}{Q}_{\text{ln}}^{f}\cdot \frac{\partial {\tilde{N}}_{m}}{\partial {\tilde{\kappa }}_{n}}\)
For their part, generalized constraints are written as, cf. [eq]:
\(N={Q}_{m}\tilde{N}{Q}_{m}^{T}\)
\(M={Q}_{f}\tilde{M}{Q}_{f}^{T}\)
2.7. Calculation of the dissipation#
By definition, the dissipation power density during damage is equal to:
\(\dot{D}=\sigma :\dot{\varepsilon }-{\dot{\Phi }}_{\mathrm{ed}}^{S}=-\sum _{j=\mathrm{1,2}}\frac{\partial {\Phi }_{\mathrm{ed}}^{S}}{\partial {d}_{j}}{\dot{d}}_{j}=\sum _{j=\mathrm{1,2}}{Y}_{j}{\dot{d}}_{j}\)
In this expression, we used the definition of \({Y}_{j}\) [eq.]. In the damaging phase, the threshold functions always satisfy \({f}_{{d}_{j}}={Y}_{j}-{k}_{0}\equiv 0\). Therefore, cumulative dissipation can be calculated as:
It has been shown that the cumulative dissipation of the damage process is directly related to internal variables. All you have to do is sum the two contributions and multiply it by the threshold constant \({k}_{0}\).
The dissipation calculation is carried out by the options DISS_ELGA and DISS_ELNO of CALC_CHAMP. These fields have a single component named ENDO.
2.8. Internal variables of the model#
The model requires two internal variables, \({d}_{1}\) and \({d}_{2}\) (corresponding to Code_Aster variables \(\mathrm{V1}\) and \(\mathrm{V2}\)), which represent damage on the upper side and on the lower side, respectively. The distinction between the upper and lower faces is made through the orientation of the local coordinate system of each Gauss point. Thus, the lower face is damaged with \(\dot{{d}_{2}}>0\) for positive curvatures, and the upper face with \(\dot{{d}_{1}}>0\) for negative curvatures, which is directly deduced from the definition of \({\xi }_{f}\), (cf. section 2.1).
In any case, the choice of the orientation of the coordinate system at a Gauss point does not affect the final result in terms of movements and rotations. There can be an impact on the interpretation of damage and generalized constraints if local orientations are not consistent within a structure. It is strongly recommended to ensure this consistency by entering the keyword ANGL_REP with the nautical angles of the local coordinate system, (see [U4.42.01]):
AFFE_CARA_ELEM (COQUE = _F (ANGL_REP = \((a,b)\)))
In addition to the variables \(\mathrm{V1}\) and \(\mathrm{V2}\), \(\mathrm{V3}\) and \(\mathrm{V4}\), with binary values (\(0\) or \(1\)), are also introduced, which indicate the instantaneous evolution of \(\mathrm{V1}\) and \(\mathrm{V2}\). More precisely, \(\mathrm{V3}\) is \(1\) when \(\mathrm{V1}\) evolves and \(0\) otherwise. Likewise, \(\mathrm{V4}\) is 1 when \(\mathrm{V2}\) evolves and \(0\) otherwise.
We also introduce \(\mathit{V5}\), \(\mathit{V6}\) and \(\mathit{V7}\), whose role is to measure the relative stiffness reduction of the reinforced concrete slab in a rational manner, for example by visualization at each material point:
\(\mathrm{V5}=1-\frac{1}{2}(\frac{1+{\gamma }_{\mathrm{mt}}{d}_{1}}{1+{d}_{1}}+\frac{1+{\gamma }_{\mathrm{mt}}{d}_{2}}{1+{d}_{2}})\) \(\mathrm{V6}=1-\frac{1}{2}(\frac{{\alpha }_{c}+{\gamma }_{\mathrm{mc}}{d}_{1}}{{\alpha }_{c}+{d}_{1}}+\frac{{\alpha }_{c}+{\gamma }_{\mathrm{mc}}{d}_{2}}{{\alpha }_{c}+{d}_{2}})\) and \(\mathrm{V7}=1-\mathrm{Max}(\frac{1+{\gamma }_{f}{d}_{1}}{1+{d}_{1}},\frac{1+{\gamma }_{f}{d}_{2}}{1+{d}_{2}})\)
respectively in tension, in compression and in flexure. These variables will always be between \(0\) and the respective \(1-{\gamma }_{\mathrm{.}}\), and will always be increasing, being zero in the absence of damage. These variables are more « meaningful » than variables \(\mathrm{V1}\) and \(\mathrm{V2}\).
The internal variables \(\mathit{V8}\), \(\mathit{V9}\), \(\mathit{V10}\) and \(\mathit{V11}\) are introduced in order to facilitate the post-treatments of the reinforcements. These variables depend on limit deformations EPSI_ELS (limit to ELS) and EPSI_LIM (limit to ELU), which are provided via DEFI_MATERIAU/GLRC_DM.
\(V8\): ACIXELS: ratio between the deformation of the steel in the X direction (maximum between the lower and upper sheet) and the deformation EPSI_ELS
\(V9\): ACIXELU: ratio between the deformation of the steel in the X direction (maximum between the lower and upper sheet) and the deformation EPSI_ELU. It should be noted that model GLRC_DM does not model the plasticization of steels and is therefore not predictive in this field.
\(V10\): ACIYELS: ratio between the deformation of the steel in the Y direction (maximum between the lower and upper sheet) and the deformation EPSI_ELS
\(V11\): ACIYELU: ratio between the deformation of the steel in the Y direction (maximum between the lower and upper sheet) and the deformation EPSI_ELU. It should be noted that model GLRC_DMne does not model the plasticization of steels and is therefore not predictive in this field.
For concrete, the internal variables \(\mathit{V12}\) and \(\mathit{V13}\) are introduced.
\(V12\): BETSUP: ratio between the lowest main deformation (in compression) of the concrete on the upper side and the limiting deformation of the concrete in compression EPSI_C.
\(V13\): BETINF: ratio between the lowest main deformation (in compression) of concrete on the lower side and the limiting deformation of concrete in compression EPSI_C.
To plot the maximum deformations, we introduce \(\mathit{V14}\), \(\mathit{V15}\) and \(\mathit{V16}\):
\(V14\): TRAMAX: maximum temporal deformation under traction
\(V15\): COMMAX: maximum temporal deformation in compression
\(V16\): FLEMAX: maximum temporal deformation during bending
To assess the validity of the parameters, the calculation of the following errors is proposed:
\(V17\): ERRCOM: error expressed as a percentage between the areas of the curve under theoretical compression (defined by EPSI_C and FCJ) and the approximate curve (NYC, GAMMA_C) between 0 and COMMAX
\(V18\): ERRFLE: error expressed as a percentage between the areas of the curve in theoretical bending (calculation of reinforced concrete section) and the approximate curve (MFY, GAMMA_F) between 0 and FLEMAX
2.9. Use of the GLRC_DM model in thermomechanical situations#
It is assumed that the temperature distribution is refined in thickness: this is admissible for a plate in a stationary thermal situation. The thermal expansion coefficient is assumed to be that of concrete, which is clearly the majority in the cross section: it is therefore directly used in membrane and plate flexure situations.
It is also accepted that the material coefficients and parameters do not depend on the temperature in the range studied.
This fine temperature distribution in the thickness therefore results in membrane deformation and a variation in thermal curvature, resulting in a simple shift in membrane strains and curvatures as a function of the temperatures at the wall of the plate.
Thus, once this discrepancy has occurred, there is no direct impact on the expression of the law of behavior or on the calculation of irreversible changes. Two tests verify the results obtained, cf. § 4.
Model GLRC_DM can therefore be used in stationary thermomechanical situations without modifications, within this framework of hypotheses.