2. Formulation

In this chapter, we present the various equations governing the shell deformation problem in the framework of a theory of large transformations.

2.1. Geometry of solid shell elements

The volume shell is represented by volume \(\Omega\) (set of points \(Q({\xi }_{3}\mathrm{\ne }0)\)) built around the mean surface \(\omega\) (set of points \(P({\xi }_{3}\mathrm{=}0)\)). At any point \(Q\) of \(\Omega\), we build a local orthonormal coordinate system \(\left[{t}_{1}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{:}{t}_{2}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{:}n({\xi }_{1},{\xi }_{2})\right]\). Vector \(n({\xi }_{1},{\xi }_{2})\) represents the normal to surface \(\omega\).

Figure 2.1-a: Solid shell. Local reference configuration guidelines

In the initial configuration, the position of any point \(Q\) normal to the mean surface can be expressed, as a function of the position of the center of gravity \(P\) of the normal fiber, in the following way:

\({x}_{Q}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{=}{\xi }_{P}({\xi }_{1},{\xi }_{2})+{\xi }_{3}\frac{h}{2}n({\xi }_{1},{\xi }_{2})\)

2.2. Kinematics of solid shells

Figure 2.2-a: Volumic shell.Major transformations of a fiber that was initially normal to the average surface

In the deformed configuration, the position of point \(Q\) can also be expressed as a function of the position of point \(P\):

\({x}_{Q}^{\varphi }({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{=}{x}_{P}^{\varphi }({\xi }_{1},{\xi }_{2})+{\xi }_{3}\frac{h}{2}{n}^{\varphi }({\xi }_{1},{\xi }_{2})\)

where \({n}^{\varphi }\) is the unit vector obtained by large rotation of normal \(n\).

The vector \({n}^{\varphi }\) is not necessarily normal to the deformed mean surface, due to the transverse shear deformation. It is linked to the initial normal vector by the relationship:

\({n}^{\varphi }\mathrm{=}\Lambda ({\xi }_{1},{\xi }_{2})n\)

\(\Lambda\) is the orthogonal operator of the large rotation around the vector \(\Theta\), with an angle \(\theta\), undergone by the fiber, which was initially normal to the mean surface, the expression of which is given by:

\(\Lambda \mathrm{=}\text{exp}\mathrm{[}\Theta \mathrm{\times }\mathrm{]}\mathrm{=}\text{cos}\theta \mathrm{[}I\mathrm{]}+\frac{\text{sin}\theta }{\theta }\mathrm{[}\Theta \mathrm{\times }\mathrm{]}+\frac{1\mathrm{-}\text{cos}\theta }{{\theta }^{2}}\mathrm{[}\Theta \mathrm{\otimes }\Theta \mathrm{]}\)

where \(\mathrm{[}\Theta \mathrm{\times }\mathrm{]}\) is the antisymmetric operator for the total rotation vector \(\Theta\) whose matrix expression is:

\(\mathrm{[}\Theta \mathrm{\times }\mathrm{]}\mathrm{=}\left[\begin{array}{ccc}0& \mathrm{-}{\Theta }_{z}& {\Theta }_{y}\\ {\Theta }_{z}& 0& \mathrm{-}{\Theta }_{x}\\ \mathrm{-}{\Theta }_{y}& {\Theta }_{x}& 0\end{array}\right]\)

and \(\mathrm{[}\Theta \mathrm{\otimes }\Theta \mathrm{]}\) is the symmetric operator given by \(\mathrm{[}\Theta \mathrm{\otimes }\Theta \mathrm{]}\mathrm{=}\Theta {\Theta }^{T}\).

More details about large rotations and their numerical processing can be found in [bib1] or [R5.03.40]. You can also write:

\(\begin{array}{c}{t}_{1}^{\varphi }\mathrm{=}\Lambda ({\xi }_{1},{\xi }_{2}){t}_{1}\\ \\ {t}_{2}^{\varphi }\mathrm{=}\Lambda ({\xi }_{1},{\xi }_{2}){t}_{2}\end{array}\)

The virtual variation of the high rotation operator can be expressed in the form:

\(\delta \Lambda \mathrm{=}\left[\delta w\mathrm{\times }\right]\Lambda\)

where \(\mathrm{[}\delta w\mathrm{\times }\mathrm{]}\) is the antisymmetric operator of the spatial virtual rotation vector \(\delta w\) which is also the rotation part of the test functions:

\(\left[\delta \mathrm{w}\mathrm{\times }\right]\mathrm{b}\mathrm{=}\delta \mathrm{w}\mathrm{\wedge }\mathrm{b}\mathrm{\forall }b\mathrm{\in }{\mathrm{ℝ}}^{3}\)

Its matrix expression is:

\(\left[\delta w\mathrm{\times }\right]\mathrm{=}\left[\begin{array}{ccc}0& \mathrm{-}\delta {w}_{z}& \delta {w}_{y}\\ \delta {w}_{z}& 0& \mathrm{-}\delta {w}_{x}\\ \mathrm{-}\delta {w}_{y}& \delta {w}_{x}& 0\end{array}\right]\)

We can also express the iterative variation of the large rotation operator in the form:

\(\Delta \Lambda \mathrm{=}\left[\Delta w\mathrm{\times }\right]\Lambda\)

where \(\Delta w\) is the spatial iterative rotation vector, which is also the rotation part of the solution to the system of linearized equations.

This vector can be linked to the total iterative rotation vector. So we have the relationships:

\(\Delta w\mathrm{=}T(\Theta )\Delta \Theta\) and \(\delta w\mathrm{=}T(\Theta )\delta \Theta\)

where \(T(\Theta )\) is the differential rotation operator, whose expression in terms of the total rotation vector is given by:

\(T(\Theta )\mathrm{=}\frac{\text{sin}\theta }{\theta }\mathrm{[}I\mathrm{]}\mathrm{-}\frac{1\mathrm{-}\text{cos}\theta }{{\theta }^{2}}\mathrm{[}\Theta \mathrm{\times }\mathrm{]}+\frac{\theta \mathrm{-}\text{sin}\theta }{{\theta }^{3}}\mathrm{[}\Theta \mathrm{\otimes }\Theta \mathrm{]}\)

This matrix has the same values and eigenvectors as the matrix \(\Lambda\) and verifies the relationship:

\(T(\Theta )\mathrm{=}\Lambda (\Theta ){T}^{T}(\Theta )\)

Moreover, the iterative variation of the virtual rotation matrix can be put in the form of:

\(\Delta \delta \Lambda \mathrm{=}\left[\delta w\mathrm{\times }\right]\left[\Delta w\mathrm{\times }\right]\Lambda\)

The total displacement of point \(Q\) on the fiber can be linked to the displacement of the center of gravity \(P\):

\({u}_{Q}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{=}{u}_{P}({\xi }_{1},{\xi }_{2})+{\xi }_{3}\frac{h}{2}({n}^{\varphi }({\xi }_{1},{\xi }_{2})\mathrm{-}n({\xi }_{1},{\xi }_{2}))\)

In order to arrive at a system of linearized equations, obtained from the weak form of equilibrium, we need to calculate different differential variations of this total displacement. Virtual displacement has the following expression:

\(\delta {u}_{Q}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{=}\delta {u}_{P}({\xi }_{1},{\xi }_{2})+{\xi }_{3}\frac{h}{2}\delta w({\xi }_{1},{\xi }_{2})\mathrm{\wedge }{n}^{\varphi }({\xi }_{1},{\xi }_{2});\delta n\mathrm{=}0\)

The expression for iterative displacement is:

\(\Delta {u}_{Q}({\xi }_{1},{\xi }_{2},{x}_{3})\mathrm{=}\Delta {u}_{P}({\xi }_{1},{\xi }_{2})+{x}_{3}\frac{h}{2}\Delta w({\xi }_{1},{\xi }_{2})\mathrm{\wedge }{n}^{\varphi }({\xi }_{1},{\xi }_{2});\Delta n\mathrm{=}0\)

The iterative variation of virtual displacement is expressed as:

\(\Delta \delta {u}_{Q}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{=}{\xi }_{3}\frac{h}{2}\delta w({\xi }_{1},{\xi }_{2})\mathrm{\wedge }(\Delta w({\xi }_{1},{\xi }_{2})\mathrm{\wedge }{n}^{\varphi }({\xi }_{1},{\xi }_{2}))\)

Note: The proposed formulation is still limited to rotations less than \(2\pi\). This limit is due to the particular choice of updating the large rotations implemented in Code_Aster. This is due to the non-bijection between the total rotation vector and the orthogonal rotation matrix.

2.3. Law of behavior

We consider a linear hyperelastic law of behavior: the local Piola‑Kirchhoff constraints of the second kind are proportional to the local Green-Lagrange deformations:

\(\tilde{S}\mathrm{=}D\tilde{E}\)

Hereinafter, the symbol \(\stackrel{}{˜}\) designates the quantities expressed in the \(\left[{t}_{1}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{:}{t}_{2}({\xi }_{1},{\xi }_{2},{\xi }_{3})\mathrm{:}n({\xi }_{1},{\xi }_{2})\right]\) orthonormal coordinate system.

The linear elastic behavior matrix under plane stresses is written as follows:

\(D\mathrm{=}\left[\begin{array}{ccccc}\frac{E}{1\mathrm{-}{\nu }^{2}}& \frac{\nu E}{1\mathrm{-}{\nu }^{2}}& 0& 0& 0\\ & \frac{E}{1\mathrm{-}{\nu }^{2}}& 0& 0& 0\\ & & \frac{E}{2(1+\nu )}& 0& 0\\ & \text{sym}& & \frac{\text{Ek}}{2(1+\nu )}& 0\\ & & & & \frac{\text{Ek}}{2(1+\nu )}\end{array}\right]\)

\(E\) being the Young’s modulus, \(\nu\) the Poisson’s ratio, and \(k\) the transverse shear correction coefficient.

In the local coordinate system, the second type of Piola-Kirchhoff stress state is plane \(({\tilde{S}}_{\text{nn}}\mathrm{=}0)\) and can be characterized by a 5-component vector:

\(\tilde{S}\mathrm{=}(\begin{array}{c}{\tilde{S}}_{{t}_{1}{t}_{1}}\\ {\tilde{S}}_{{t}_{2}{t}_{2}}\\ {\tilde{S}}_{{t}_{1}{t}_{2}}\\ \\ {\tilde{S}}_{{t}_{1}n}\\ {\tilde{S}}_{{t}_{2}n}\end{array})\)

The vector of Green-Lagrange deformations is also expressed in the local coordinate system by a vector with 5 components:

\(\tilde{E}\mathrm{=}(\begin{array}{c}{\tilde{E}}_{{t}_{1}{t}_{1}}\\ {\tilde{E}}_{{t}_{2}{t}_{2}}\\ {\tilde{\gamma }}_{{t}_{1}{t}_{2}}\\ \\ {\tilde{\gamma }}_{{t}_{1}n}\\ {\tilde{\gamma }}_{{t}_{2}n}\end{array})\)

Here, we ignored the term \({\tilde{E}}_{\text{nn}}\), which is normal on the average surface and which is not necessarily zero. This is a consequence of the hypothesis of plane stresses.

2.3.1. Taking into account transverse shear

The correction of the transverse shear stress is carried out by extending the energy equivalences determined in the case of small deformations and small displacements [R3.07.03].