Understanding the strengths in shells

1 Introduction

The finite element method permits to solve the governing equations of a structural model by using numerical methods over one model discretization. When the structure is thin (i.e. it is a three-dimensional solid whose thickness is very small when compared with other dimensions) it is usual to model it using a shell finite element model. The results of this analysis are the displacements and rotations in every node of the mesh and the strengths (i.e. the values resulting of the integration of the stress tensor components through the shell thickness) in some interior element positions.

In order to get valuable information of the analysis, it is necessary a good understanding and interpretation of the strengths. This document will explain what the strengths in shells are and how can they be used for dimensioning or evaluating one structural model.

We are going to define mathematically the strengths in shells as they are used in the standard shell Reissner-Mindlin theory. The set of strength values describes fully the stress tensor for every point in the shell mid plane.

Image 1. Shell is represented by its mid plane

2 Definitions

One structure can be considered a Shell if it is thin (have one dimension much smaller than the other two).

Some definitions:

Shell

A shell is a finite element model of a three-dimensional structure whose thickness is very small when compared with other dimensions, and in structural terms, is characterized by some simplifying hypotheses like the Reissner-Mindlin hypotheses. The geometry of the shell model is completely defined by its middle surface (plane or curved surface in the 3D) and the thickness and normal in every point on this middle surface

Thickness

for every point in the structure, the thickness is the width of the model through the small dimension

Local axes

for every point in the structure it is necessary to define an attached local axes system X’, Y’ and Z’ where axe Z’ will always be orthogonal to the Shell. X’, Y’ can be any axes but must be defined for every point in the Shell prior to the calculation. All strengths will be expressed on these axes

Section

is the part of the structure which is the intersection between the structure and a plane orthogonal to the Shell mid plane

and then, the Reissner-Mindlin hypotheses are defined as follows:

• One planar section, after deformation will continue to be planar. The deformed section does not generally remain orthogonal to the mid-surface
• The normal stress through the thickness is ignored (${\textstyle \sigma _{z'}=0}$)
• There is a linear variation of displacement across the plate thickness and the plate thickness does not change during deformation

3 Definition of the strengths in the shell

In any internal point of the structure the stress tensor represents the stresses in that point. As we have simplified the shell structure by the mid plane, it is necessary to define some values that are calculated as the integrals of the stress tensor through the shell thickness so as they represent completely, under the accepted hypotheses, the stresses in the mid plane.

Image 2. Strengths in the shell

The strengths are three axial forces in the plane, three bending/twisting moments and two shears

Table 1. Definition of strengths

Axial force in X’	Axial force in Y’	Axial force in XY’
${\textstyle N_{x}=\int _{-t/2}^{t/2}{\sigma _{x}dz}}$	${\textstyle N_{y}=\int _{-t/2}^{t/2}{\sigma _{y}dz}}$	${\textstyle N_{xy}=\int _{-t/2}^{t/2}{\tau _{xy}dz}}$
Bending moment X’	Bending moment Y’	Twisting moment XY’
${\textstyle M_{x}=\int _{-t/2}^{t/2}{\sigma _{x}zdz}}$	${\textstyle M_{y}=\int _{-t/2}^{t/2}{\sigma _{y}zdz}}$	${\textstyle M_{xy}=\int _{-t/2}^{t/2}{\tau _{xy}zdz}}$
Shear force in X’	Shear force in Y’
${\textstyle Q_{x}=\int _{-t/2}^{t/2}{\tau _{xz}dz}}$	${\textstyle Q_{y}=\int _{-t/2}^{t/2}{\tau _{yz}dz}}$

All the integrals are in Z’ axe from ${\textstyle -t/2}$ to ${\textstyle t/2}$, being ${\textstyle t}$ the thickness of the shell in that point. Z’ is the local axis defined as parallel to the normal to the surface at the point.

We need to consider that the real ${\textstyle Q_{x}}$ and ${\textstyle Q_{y}}$ have a parabolic shape through the shell thickness with null values in the top and bottom layers. However, we accept the simplification of considering the shear as constant through all the thickness.

Image 3. Real parabolic shear stress is simplified as constant through the thickness

4 Finite element analysis

The Finite element analysis is a method to solve numerically the differential equation that governs the structural behavior of the structure. To do so, it is necessary to discretize the model into a mesh, which is a collection of elements connected between them by nodes. The variables to solve in the problem based on shell elements are the displacements and rotations of every node of the mesh. From these values, it is possible to interpolate the value of the displacements for every point over the structure and it is also possible to calculate the strengths as a derived result of the displacements.

Image 4. The discretization is the triangle mesh

Once discretized, the finite element problem becomes:

$$Ka=f$$ (1)

where ${\textstyle K}$ is the stiffness matrix, ${\textstyle f}$ is the external forces vector and ${\textstyle a}$ is the displacements vector. As every node has 6 degrees of freedom (${\textstyle d_{x},d_{y},d_{z},\theta _{x},\theta _{y},\theta _{z}}$), the dimension of ${\textstyle a}$ is 6 x number of nodes.

If the problem is linear, ${\textstyle K}$ and ${\textstyle f}$ do not depend on ${\textstyle a}$ and the problem can be solved by simply solving the linear system once.

If the problem is non linear, ${\textstyle K}$ or ${\textstyle F}$ depend on ${\textstyle a}$ and then it is necessary to solve the problem by iterating through a non-linear solver.

If the problem is dynamic, the equation to solve becomes:

$$M{\ddot {a}}+C{\dot {a}}+Ka=f$$ (2)

see [1] for more details.

5 Results of the finite element analysis

The finite element analysis gets reduced to the resolution of a system of equations where the variables are the displacements and rotations defined in every node of the mesh. As these variables are defined in the nodes that are shared between neighboring elements, the displacements and rotations fields are continuum through the elements by the natural definition of the problem.

As described before, the strengths are the integration of the stress tensor through the thickness of the shell and the stresses are obtained from the derivatives of the displacements. Both the stresses and the strength are calculated inside every element. Specifically, they are calculated in the integration points used to compute the numerical integration inside the element during the analysis. These points are the so-called Gauss points which are used by the quadrature rule to increase the precision of the numerical integration.

As the strengths are calculated inside every element, this leads to two effects that are important to consider:

1. The strengths fields are discontinuous between elements
2. When displaying the strengths in a postprocessing program, the software linearly extrapolates the values calculated internally to new values in every node of the element boundary. For relative small differences of the strength values on the Gauss points of the element, the lineal extrapolation process is likely to calculate invalid nodal values (e.g. too large values, or a negative value for a variable that is known to be positive).

6 Smoothing

As the strengths fields are discontinuous between elements, in some cases it might be convenient to calculate another field, similar to the original one, that is continuous. This is only possible between elements that share equal local axes definition. The smoothed field is calculated by comparing the values in every node from each of the neighboring elements. Some kind of mean value is obtained and used as the value for the node.

This correction is generally good, as it provides a value that is closer, in most cases, to the real value. However, in extreme cases and specially near the shell boundary or in points over external elements or reinforcements (beams), it can lead to incorrect values. This effect can be maximized by the addition of the two correction effects: the smoothing and the linearly extrapolation.

It is important to visualize and compare the non-smoothed results and the smoothed ones in order to understand better the strengths field.

Table 2. Moments in coarse mesh. Gauss points results

Table 3. Moments in coarse mesh. Smoothed results

Table 4. Moments in fine mesh. Gauss points results

Table 5. Moments in fine mesh. Smoothed results

The results of the examples displayed in the figures above are summarized in Table 6. We can see that the error, specially for the coarse mesh, is quite high.

Table 6. Some results for the different meshes

Mesh	Gauss points/Smooth	${\textstyle M_{x,border}}$	${\textstyle M_{x,max}}$	${\textstyle M_{x,min}}$
coarse	GP	-2.30	20.0	-77
coarse	SM	-2.30	13.0	-75
fine	GP	-0.17	5.8	-74
fine	SM	-0.17	4.6	-73

7 Main strengths and main stresses

For any tensor defined in a local axes system, there exists another local axes system (the so-called principal axes) for which the tensor is diagonal (i.e. its non-diagonal values are zero). In this system, the values in the diagonal are maximized and minimized.

The principal axes are defined by the eigenvalues and eigenvectors of the matrix.

The local axes systems for the main moments, main axial strengths and main stresses need not to be equivalent and typically will be different between them.

It is important to remind that the stress and the corresponding strength tensors are symmetric (i.e ${\textstyle N_{xy}=N_{yx}}$, ${\textstyle M_{xy}=M_{yx}}$, ${\textstyle \tau _{xy}=\tau _{yx}}$, etc.).

Table 7. Strengths and stresses in principal axes

Strengths in local axes ${\textstyle LA_{0}}$	Main strengths
Axial strengths as membrane components
${\textstyle {\begin{bmatrix}N_{x}&N_{xy}\\N_{yx}&N_{y}\end{bmatrix}}}$	${\textstyle {\begin{bmatrix}N_{11}&0\\0&N_{22}\end{bmatrix}}}$	Local axes axial (${\textstyle LA_{A}}$)
Bending and twisting moments
${\textstyle {\begin{bmatrix}M_{x}&M_{xy}\\M_{yx}&M_{y}\end{bmatrix}}}$	${\textstyle {\begin{bmatrix}M_{11}&0\\0&M_{22}\end{bmatrix}}}$	Local axes moments (${\textstyle LA_{M}}$)
Stress tensor
${\textstyle {\begin{bmatrix}\sigma _{x}&\tau _{xy}&\tau _{xz}\\\tau _{yx}&\sigma _{y}&\tau _{yz}\\\tau _{zx}&\tau _{zy}&\sigma _{z}\end{bmatrix}}}$	${\textstyle {\begin{bmatrix}\sigma _{11}&0&0\\0&\sigma _{22}&0\\0&0&\sigma _{33}\end{bmatrix}}}$	Local axes stresses (${\textstyle LA_{S}}$)

Where ${\textstyle LA_{0}}$, ${\textstyle LA_{A}}$, ${\textstyle LA_{M}}$ and ${\textstyle LA_{S}}$ will represent typically different local axes systems with a common Z’ axe.

Image 5. Main strengths in ${\textstyle LA_{M}}$ and ${\textstyle LA_{M}}$

8 Stresses in top and bottom layers

In the shell model, we assume that stresses vary linearly through the thickness of the Shell. As a consequence, both the maximum and the minimum of any stress component will be located in either the top or the bottom layer for any shell point.

Note: shear stresses in y’ and z’ have a real parabolic shape through the Shell thickness, with a maximum in the middle plane of the shell. However, in our simplified model, we assume a constant value for the shear stresses in the thickness. So, maximum and minimum are also in the external layers.

Top and bottom stresses are calculated in the following way:

$$\sigma _{x}={\frac {N_{x}}{A_{b}}}\pm {\frac {M_{x}}{W_{b}}}\qquad \sigma _{y}={\frac {N_{y}}{A_{b}}}\pm {\frac {M_{y}}{W_{b}}}\qquad \sigma _{z}=0$$ (3)

$$\tau _{xy}={\frac {N_{xy}}{A_{b}}}\pm {\frac {M_{xy}}{W_{b}}}\qquad \tau _{xz}={\frac {Q_{x}}{A_{b}}}\qquad \tau _{yz}={\frac {Q_{y}}{A_{b}}}$$ (4)

where:

$$W_{b}={\frac {I_{z}}{z_{max}}}={\frac {bh^{2}}{6}}$$ (5)

$$A_{b}=bh$$ (6)

the values for ${\textstyle b}$ and ${\textstyle h}$ for shells are: ${\textstyle b=1.0}$ and ${\textstyle h=t}$

9 Dimensioning metallic materials

Based on the stresses tensor calculated before, it is possible to calculate the Von Mises stress for the top and bottom layers following the formula:

$$\sigma _{vm}={\sqrt {\frac {(\sigma _{I}-\sigma _{II})^{2}+(\sigma _{II}-\sigma _{III})^{2}+(\sigma _{III}-\sigma _{I})^{2}}{2}}}$$ (7)

where ${\textstyle \sigma _{I},\sigma _{II},\sigma _{III}}$ are the main stresses associated to the stress tensor ${\textstyle \sigma _{x},\sigma _{y},\sigma _{z},\tau _{xy},\tau _{xz},\tau _{yz}}$.

The Von Mises stress as defined in (7) can be used as a yield criteria for the material. This is commonly used as yield criteria for metallic materials.

10 Dimensioning structural concrete

Structural concrete is a composite material made of concrete and steel bars. Several normatives, like Eurocode and ACI and other national regulations, offer guidelines to dimension and check the concrete cross section and the amount of steel.

For a model where the main axial forces axes ${\textstyle LA_{N}}$ and the main momentum axes ${\textstyle LA_{M}}$ are coincident with the element local axes ${\textstyle L_{0}}$, the strengths in the section will be ${\textstyle N_{x'}}$, ${\textstyle M_{x'}}$ and ${\textstyle Q_{x'}}$ in the ${\textstyle x'}$ direction and ${\textstyle N_{y'}}$, ${\textstyle M_{y'}}$ and ${\textstyle Q_{y'}}$ in the ${\textstyle y'}$ direction. With these strengths, it is possible to dimension the structure by using the regulation standards.

In the general case, there is an additional twisting moment ${\textstyle M_{xy}}$ that is different to zero. It is not in the security side to dimension the section without considering this twisting moment in the verification of the main steel bars in both directions. It is necessary to use some method, that satisfies the equilibrium and geometric compatibility in the section and that has some bearing strengths greater that the external actions.

One commonly used method is the Wood and Armer. In this method, the bearing ${\textstyle m_{ux}}$ and ${\textstyle m_{uy}}$ are increased based on the value of ${\textstyle m_{xy}}$. The resulting values are used for dimensioning the structure.

The CEB-FIP Model Code proposes and alternative method where the thickness is divided into three layers and the strengths are converted into pairs of forces and shear forces for dimensioning.

10.1 Wood and Armer method

The Wood and Armer method is based on the normal moment yield criterion (Johansen’s yield criterion) and proposes that the dimensioning of the maximum bear moments for some given directions ${\textstyle x}$ and ${\textstyle y}$ can be calculated as follows:

${\textstyle m_{ux}=m_{x}\pm \left|m_{xy}\right|}$ (8)

${\textstyle m_{uy}=m_{y}\pm \left|m_{xy}\right|}$ (9)

if ${\textstyle m_{x}}$ and ${\textstyle m_{y}}$ have different sign, and ${\textstyle m_{x}}$ is bigger then:

${\textstyle m_{ux}=m_{x}\pm \left|{\frac {m_{xy}^{2}}{m_{y}}}\right|}$ (10)

${\textstyle m_{uy}=0}$ (11)

10.1.1 Theoretical justification of the Wood and Armer method

Given external moments ${\textstyle m_{x},m_{y}}$ and ${\textstyle m_{xy}}$, the moments in direction ${\textstyle n}$ where the angle between ${\textstyle x}$ and ${\textstyle n}$ is ${\textstyle \theta }$ are the following:

${\textstyle m_{n}=m_{x}\cos ^{2}\theta +m_{y}\sin ^{2}\theta +2m_{xy}\cos \theta \sin \theta }$ (12)

${\textstyle m_{nt}=(m_{y}-m_{x})\cos \theta \sin \theta +m_{xy}(\cos ^{2}\theta -sin^{2}\theta )}$ (13)

where ${\textstyle m_{n}}$ and ${\textstyle m_{nt}}$ are the bending and twisting moments in the ${\textstyle n}$ direction.

On the other side, if ${\textstyle m_{ux}}$ and ${\textstyle m_{uy}}$ are the moment capacities in axes x,y then, the bearing capacity of the section in direction ${\textstyle n}$ where the angle between ${\textstyle x}$ and ${\textstyle n}$ is ${\textstyle \theta }$ is the following:

${\textstyle m_{un}=m_{ux}\cos ^{2}\theta +m_{uy}\sin ^{2}\theta }$ (14)

Imposing that ${\textstyle m_{un}}$ in (14) must be bigger that ${\textstyle m_{n}}$ in (12) for any direction, it is possible to justify that the Wood and Armer proposal is an upper bound of the capacity for any direction.

The justifications and details of the method are described in [2] and [3].

10.2 CEB-FIP Model Code

The slab model is decomposed into three layers, each with ${\textstyle 1/3}$ of the slab thickness. Bending moments are decomposed into pairs of forces and applied to the layers as normal forces. Twisting moments are equally decomposed into pairs of force that are applied as shear forces in the two external layers. Finally, the shear forces are applied to the central layer for verification.

11 Dimensioning composite laminate material

A composite laminate material or fiber reinforced plastic (FRP) is built by adding several layers of fiber that get glued together by the use of a resin. In order to dimension or check its bearing capacity, it is typically necessary to evaluate a failure criterion for every one of its constitutive layers or plies.

Every single layer can be defined by its components of fiber and resin and are normally orthotropic, having different resistance values for every one of the orthogonal directions. This anisotropic composite materials will also have different resistance capacities in tension and compression.

It is common to use the Tsai–Wu failure criterion [4] in order to check the capacity of each single layer.

Image 6. Analysis of a laminate ply

11.1 Tsai–Wu failure criterion

The criterion is defined as follows:

$$\sum _{i=1}^{6}F_{i}\sigma _{i}+\sum _{i=1}^{6}\sum _{j=1}^{6}F_{ij}\sigma _{i}\sigma _{j}\leq 1$$ (15)

where ${\textstyle \sigma _{i}}$ is one of ${\textstyle (\sigma _{xx},\sigma _{yy},\sigma _{zz},\sigma _{yz},\sigma _{xz},\sigma _{xy})}$ and ${\textstyle F_{i}}$ and ${\textstyle F_{ij}}$ are some coefficients calculated from laboratory results.

For plane stress, where ${\textstyle \sigma _{z}=\sigma _{yz}=\sigma _{xz}=0}$, the criteria is reduced to:

${\textstyle F_{1}\sigma _{xx}+F_{2}\sigma _{yy}+F_{11}\sigma _{xx}^{2}+}$ ${\textstyle F_{22}\sigma _{yy}^{2}+F_{66}\tau _{xy}^{2}+2F_{12}\sigma _{xx}\sigma _{yy}\leq 1}$ (16)

where:

${\textstyle F_{1}={\frac {1}{\sigma _{1t}}}-{\frac {1}{\sigma _{1c}}}}$	${\textstyle F_{2}={\frac {1}{\sigma _{2t}}}-{\frac {1}{\sigma _{2c}}}}$
${\textstyle F_{11}={\frac {1}{\sigma _{1c}\sigma _{1t}}}}$	${\textstyle F_{22}={\frac {1}{\sigma _{2c}\sigma _{2t}}}}$
${\textstyle F_{66}={\frac {1}{\tau _{12}^{2}}}}$

As it is difficult to obtain the value of ${\textstyle F_{12}}$ from a laboratory test, some references suggest to use the value ${\textstyle F_{12}=0}$. Tsai recommends to obtain its value based on the following relationship:

$$F_{12}=-0.5(F_{11}F_{22})^{0.5}$$ (17)

Where ${\textstyle \sigma _{1t}}$ and ${\textstyle \sigma _{1c}}$ are the tensile and compressive allowable stress in the longitudinal direction and ${\textstyle \sigma _{2t}}$ and ${\textstyle \sigma _{2c}}$ are the tensile and compressive allowable stress in the transverse direction. ${\textstyle \tau _{12}}$ is the allowable shear stress.

see [4] and see [5] for a detailed description.

11.2 Tsai–Wu to calculate a minimum security factor

Another way of considering the problem is that the security factor ${\textstyle SF_{CS}}$ for the combined stress in any individual layer must be ${\textstyle SF_{CS}>SF_{csiapp}}$ where ${\textstyle SF_{csiapp}}$ can be calculated based on (16) as:

${\textstyle SF_{csiapp}F_{1}\sigma _{xx}+SF_{csiapp}F_{2}\sigma _{yy}+SF_{csiapp}^{2}F_{11}\sigma _{xx}^{2}+}$ ${\textstyle SF_{csiapp}^{2}F_{22}\sigma _{yy}^{2}+SF_{csiapp}^{2}F_{66}\tau _{xy}^{2}+}$ ${\textstyle SF_{csiapp}^{2}2F_{12}\sigma _{xx}\sigma _{yy}=1}$ (18)

then:

$$SF_{csiapp}={\frac {-b\pm {\sqrt {b^{2}-4ac}}}{2a}}$$ (19)

where:

${\textstyle a=F_{11}\sigma _{xx}^{2}+F_{22}\sigma _{yy}^{2}+2F_{12}\sigma _{xx}\sigma _{yy}+F_{66}\tau _{xy}^{2}}$

${\textstyle b=F_{1}\sigma _{xx}+F_{2}\sigma _{yy}}$

${\textstyle c=-1}$

The complete process is described at [6].

References

[1] Oñate E. (1991) Cálculo de estructuras por el método de los elementos finitos. ISBN 8487867006

[2] Shin M., Bommer A., Deaton, J.B. and Alemdar, B.N. (2009) Twisting Moments in Two-Way Slabs. Concrete international. Available: https://www.academia.edu/15396812

[3] OBrien E., O’Dwyer D. (2015). A review of the Wood and Armer equations for torsion in bridge slabs. 93. 40-43

[4] Tsai S.W. and Wu E.M. (1971) A general theory of strength for anisotropic materials. J Compos Mater, 5. pp. 58-80

[5] Hopkins P. (2005) Benchmarks for Membrane and Bending Analysis of Laminated Shells. ISBN 1874376085

[6] Bureau Veritas Marine & Offshore SA (2019) Engineering, production and life cycle management for the complete construction of large-length FIBRE-based SHIPs