Mindlin-Reissner Analytical Model with Curvature for Tunnel Ventilation Shafts Analysis

Abstract

The formulation and analytic solution of a new mathematical model with constitutive curvature for analysis of tunnel ventilation shaft wall is proposed. Based on the Mindlin–Reissner theory for thick shells, this model also takes into account the shell constitutive curvature and considers an expression of the shear correction factor variable $\left({\alpha }_n\right)$ in terms of the thickness $\left(h\right)$ and the radius of curvature $\left(R\right)$. The main advantage of the proposed model is that it has the possibility to analyze thin, medium and thick tunnel ventilation shafts. As a result, two comparisons were made: the first one, between the new model and the Mindlin–Reissner model without constitutive curvature with the shear correction factor ${\alpha }_n=5/$6 as a constant, and the other, between the new model and the tridimensional numerical models (solids and shells) obtained by finite element method for different slenderness ratios $(h/R)$. The limitation of the proposed model is that it is to be formulated for a general linear-elastic and axial-symmetrical state with continuous distribution of the mass.

1. Introduction

Tunnel ventilation shafts are vertical structural elements for accessing underground structures such as tunnels. The structural analysis of tunnel ventilation shafts is performed through a coaxial cylinder analysis scheme, in axial-symmetrical general state (load, geometry, material and boundary condition) (Figure 1). The importance of designing and building these structures lies in the increasing needs of hydrosanitary network and new means of communications. Thus, an optimum structural design is guaranteed if a reliable analysis method is used.

In the classical shell theory by Goldenveizer and Timoshenko [1,2,3], the three-dimensional continuous solid is modelled by the middle surface of the solid plus the thickness of the shell, Equation (1). This allows the use of the resultant internal forces per unit of arc longitude: bending moment ($M_Z,M_{\theta }$), shear forces ($Q_z$) and membrane axial forces $\left(N_{\theta },N_z\right)$ (Figure 1). In addition, Figure 1 also shows the set of variables to be considered.

$$\underset{\begin{array}{c}STRESS\\ DEFORMATIONS\\ DISPLACEMENTS\end{array}}{\underbrace{f\left({\alpha }_1,{\alpha }_2,{\alpha }_3\right)\in {ℝ}^3}}=\underset{\begin{array}{c}MIDDLESURFACE\\ OFTHESHELL\end{array}}{\underbrace{f\left({\alpha }_1,{\alpha }_2\right)\in {ℝ}^2}}+\underset{\begin{array}{c}FUNCTIONfBEHAVIOR\\ INSIDEOFTHE\\ SHELLTHICKNESS\end{array}}{\underbrace{h\left({\alpha }_3\right)\in ℝ}},$$

(1)

Resultant internal forces per unit of arc longitude and middle deformations of the shell $\left({\kappa }_z,{\kappa }_{\theta },{\gamma }_{zN},{\epsilon }_{\theta },{\epsilon }_z\right)$ are related by the integration process inside the thickness of the shell [2,3]. The integration process, defines the constitutive model that will be used in the construction of the mathematical model, in terms of displacement (direct operational model) or in resultant internal forces per unit of arc longitude (inverse operational model).

Equation (2) represents the classical linear-elastic constitutive model that relates the resultant internal forces per unit of arc longitude and the middle deformations of the shell inside the classical shell [1,2,3,4], which presents the following limitations [1,2,5]:

• The model involves a bending of a flat plate $\left(\frac{1}{R}=0\right)$ plus a membrane behavior without curvature (Figure 2);
• Nonexistence of the constitutive coupling between the membrane forces and the bending moments for any slenderness ratio $\left(\frac{h}{R}\right)$;
• Shear correction factor (${\alpha }_n$) is employed in the constant shear force, independently of the shell slenderness ratio.

$$\begin{array}{c}{N}_{\theta }=K\left({\epsilon }_{\theta }+\mu {\epsilon }_z\right),N_Z=K\left({\epsilon }_z+\mu {\epsilon }_{\theta }\right),M_z=D_f{\kappa }_z\\ M_{\theta }=\mu M_z,Q_Z={\alpha }_{n}K\frac{\left(1-\mu \right)}{2}{\gamma }_{ZN},D_f=\frac{1}{12}\frac{E{h}^3}{\left(1-{\mu }^2\right)},K=\frac{Eh}{1-{\mu }^2},{\alpha }_n=\frac{5}{6}\end{array}$$

(2)

In this paper, a new Mindlin–Reissner mathematical flexural model for tunnel ventilation shafts analysis is obtained; it takes into account the shell constitutive curvature effect [2,3]; along with an update in the shear correction factor (${\alpha }_n$) for the shear distribution. A new linear-elastic constitutive relationship between the resultant internal forces per unit of arc longitude and the middle deformations was obtained for the model construction. This constitutive relationship models the constitutive curvature effect in the flexo-compression behavior of the tunnel ventilation shafts. The new model is a generalization of the Love–Kirchhoff and Mindlin–Reissner mathematical models, without constitutive curvature for the analysis scheme of the coaxial cylinder (Figure 1).

2. Coupled Constitutive Model for General Bending of a Coaxial Cylinder

The term $d{s}^{*},$ represents the differential arc increased in the shell differential arc $\left(d\theta \right)$ of the middle surface (Figure 3). Figure 3 shows the stress state in a “differential point” that belong to the surface $S^{*}\left(d{s}^{*}dz\right)$.

From the stress–strain state definition for the “differential point”, which belong to the surface $S^{*}\left(d{s}^{*}dz\right)$, Equations (3) and (4), it is possible to integrate the surfaces ($d{A}_1\mathrm{and}d{A}_2)$ in order to obtain the resultant internal forces per unit of arc longitude: bending moments ($M_Z,M_{\theta }$), shear and membrane axial forces $\left(N_{\theta },N_z\right)$, Equation (5), Figure 3. The resultant internal forces per unit of arc longitude of the middle surface of the shell in terms of the deformations are obtained by integrating with respect to the axis n. The axis (n) is orthogonal to the middle surface of the shell (Figure 3), Equations (7) and (8). For the integration process, the term $(\frac{1}{1+\frac{n}{R}})$, Equation (6), was developed in n potential-function (Taylor series expansion) [3] affecting the integrand with the substitution of the stress-strain state for the “differential point”, Equations (3)–(5).

$${\sigma }_z{}^{*}=\frac{E}{1-{\mu }^2}\left({\epsilon }_z{}^{*}+\mu {\epsilon }_{\theta }{}^{*}\right),{\sigma }_{\theta }{}^{*}=\frac{E}{1-{\mu }^2}\left({\epsilon }_{\theta }{}^{*}+\mu {\epsilon }_z{}^{*}\right),{\tau }_{zn}{}^{*}=\frac{{\alpha }_{n}E}{2\left(1+\mu \right)}{\gamma }_{zn}{}^{*}$$

(3)

$$\begin{array}{c}{\epsilon }_z{}^{*}={\epsilon }_z+n{\kappa }_z,{\epsilon }_{\theta }{}^{*}=\frac{{\epsilon }_{\theta }}{1+\frac{n}{R}},{\gamma }_{zn}{}^{*}={\gamma }_{zn}\\ n\in ℝ:\left[-\frac{h}{2},\frac{h}{2}\right]\end{array}$$

(4)

$$\begin{array}{c}{N}_z\left(\frac{f}{l}\right)={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_z{}^{*}\left(1+\frac{n}{R}\right)dn,N_{\theta }\left(\frac{f}{l}\right)={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_{\theta }{}^{*}dn\\ Q_z\left(\frac{f}{l}\right)={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\tau }_{zn}{}^{*}\left(1+\frac{n}{R}\right)dn,M_z\left(\frac{fl}{l}\right)={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_z{}^{*}n\left(1+\frac{n}{R}\right)dn\\ M_{\theta }\left(\frac{fl}{l}\right)={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_{\theta }{}^{*}ndn,{\alpha }_n=1/k_n\\ f:forcesunit;l:lengthunit\end{array}$$

(5)

$$\frac{1}{1+\frac{n}{R}}=1+{\displaystyle \sum }_{i=1}^{\infty }{\left(-\frac{n}{R}\right)}^i$$

(6)

Equations (7) and (8) show a linear-elastic coupled constitutive model between the membrane forces and the moments for any slenderness ratio of the shell. If the shell curvature $\left(\frac{1}{R}=0\right)$ is neglected in Equations (7) or (8), the classical shell model, Equation (2), is obtained. In addition, the linear-elastic coupled constitutive model, in Equations (7) and (8), is consistent with the coupled models defined by Galerkin, Novozhilov, Finkel’shtein and Lur’e, and other authors [2,3].

The deformation expressions in terms of the forces and moments, Equation (9), are obtained by inverting the coupled constitutive matrix $\left[D_C\right]$, Equation (8). It shall be noted that the symmetry on the constitutive matrices, Equations (8) and (9) guarantees the fulfillment of the shell theory general theorems [2,3,6].

Equation (10) represents the coupled constitutive matrix, Equation (8), for an unstiffened orthotropic shell. It is noteworthy that Equation (10) generalizes the constitutive model employed by authors Rotter and Sadowski [4], opening the possibility for future research on this field

$$\begin{array}{c}{N}_z=K\left({\epsilon }_z+\mu {\epsilon }_{\theta }\right)+\frac{D_f}{R}{\kappa }_z,N_{\theta }=K\left({\epsilon }_{\theta }{\alpha }_1+\mu {\epsilon }_z\right)\\ M_z=D_f\left({\kappa }_z+\frac{{\epsilon }_z}{R}\right),M_{\theta }=D_f\left(\mu {\kappa }_z-\frac{{\epsilon }_{\theta }}{R}\right),Q_Z={\alpha }_{n}K\frac{\left(1-\mu \right)}{2}{\gamma }_{ZN}\\ D_f=\frac{1}{12}\frac{E{h}^3}{\left(1-{\mu }^2\right)},K=\frac{Eh}{1-{\mu }^2},{\alpha }_1=\left(1+\frac{h^2}{12R^2}\right),{\alpha }_n\end{array},$$

(7)

$$\begin{array}{cc}\left[\begin{array}{c}\begin{array}{c}{N}_z\\ N_{\theta }\\ M_z\end{array}\\ Q_z\end{array}\right]=\left[\begin{array}{cccc}K& \mu K& \frac{D_f}{R}& 0\\ \mu K& K{\alpha }_1& 0& 0\\ \frac{D_f}{R}& 0& D_f& 0\\ 0& 0& 0& \frac{Eh{\alpha }_n}{2\left(1+\mu \right)}\end{array}\right]\left[\begin{array}{c}{\epsilon }_z\\ {\epsilon }_{\theta }\\ \begin{array}{c}{\kappa }_z\\ {\gamma }_{zn}\end{array}\end{array}\right]& \underset{Coupledconstitutivematrix}{\underbrace{D_C=\left[\begin{array}{cccc}K& \mu K& \frac{D_f}{R}& 0\\ \mu K& K{\alpha }_1& 0& 0\\ \frac{D_f}{R}& 0& D_f& 0\\ 0& 0& 0& \frac{Eh{\alpha }_n}{2\left(1+\mu \right)}\end{array}\right]}}\end{array}$$

(8)

$$\begin{array}{c}\left[\begin{array}{c}{\epsilon }_z\\ {\epsilon }_{\theta }\\ \begin{array}{c}{\kappa }_z\\ {\gamma }_{zn}\end{array}\end{array}\right]=\left[\begin{array}{cccc}{H}_{11}& H_{12}& H_{13}& 0\\ H_{21}& H_{22}& H_{23}& 0\\ H_{31}& H_{32}& H_{33}& 0\\ 0& 0& 0& \frac{\left(1+\mu \right)}{2Eh{\alpha }_n}\end{array}\right]\left[\begin{array}{c}{N}_z\\ N_{\theta }\\ \begin{array}{c}{M}_z\\ Q_z\end{array}\end{array}\right]\\ H_{11}=-\frac{R^2{\alpha }_1}{F},H_{22}=\frac{D_f}{KF}-\frac{R^2}{F},H_{33}=-\frac{12R^2\left({\alpha }_1-{\mu }^2\right)}{KF},H_{12}=H_{21}=\frac{\mu R^2}{F}\\ H_{13}=H_{31}\frac{R{\alpha }_1}{F},H_{23}=H_{32}\begin{array}{c}-\frac{\mu R}{F}\end{array},F=\frac{D_f{h}^2}{12R^2}-R^{2}Eh\end{array}$$

(9)

$$\begin{array}{cc}\left[\begin{array}{c}\begin{array}{c}{N}_z\\ N_{\theta }\\ M_z\end{array}\\ Q_z\end{array}\right]=\left[\begin{array}{cccc}{K}_1& {\mu }_{z\theta }{K}_2& \frac{D_f{}_1}{R}& 0\\ {\mu }_{z\theta }{K}_2& K_2{\alpha }_1& 0& 0\\ \frac{D_f{}_1}{R}& 0& D_f{}_1& 0\\ 0& 0& 0& \frac{E_{z}h{\alpha }_n}{2\left(1+{\mu }_{z\theta }\right)}\end{array}\right]\left[\begin{array}{c}{\epsilon }_z\\ {\epsilon }_{\theta }\\ \begin{array}{c}{\kappa }_z\\ {\gamma }_{zn}\end{array}\end{array}\right];& \begin{array}{c}{K}_1=\frac{E_{z}h}{1-{\mu }_{z\theta }{\mu }_{\theta z}}\\ \begin{array}{c}{K}_2=\frac{E_{\theta }h}{1-{\mu }_{z\theta }{\mu }_{\theta z}}\\ D_f{}_1=\frac{1}{12}\frac{E_z{h}^3}{\left(1-{\mu }_{z\theta }{\mu }_{\theta z}\right)}\end{array}\end{array}\end{array}$$

(10)

The shear correction factor (${\alpha }_n$) of constitutive equations, Equations (1), (5), (7)–(9), for parabolic and uniform distributions [7] are obtained from the tangential stress of the energetic equilibrium. If a rectangular and homogeneous section (plate) is considered, the shear correction factor per unit of length is ${\alpha }_n=\frac{5}{6}$ Equation (11) [7], which is the default value of the FEM software [8,9], and corresponds to the classic constitutive model for plates in the Mindlin–Reissner bending state, Equation (1).

$$\begin{array}{c}{\alpha }_n=\frac{I^2}{A_1}{\left[{{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\left[g_1\left(n\right)\right]}^{2}dn\right]}^{-1}\\ g_1\left(n\right)={{\displaystyle \int }}_{-\frac{h}{2}}^{z}ndn,A_1={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}dn,I={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{n}^{2}dn\end{array}\}{\alpha }_n=\frac{5}{6}$$

(11)

By inserting in Equation (11) the shell curvature geometrical relations (Figure 3), an expression for the shear correction factor variable $\left({\alpha }_n\right)$ per unit of arc longitude in cross-cylindrical and homogeneous sections is obtained, Equations (12) and (13), which depends on the thickness and the radius of curvature of the shell (Figure 3).

$$\begin{array}{c}{\alpha }_n=\frac{I^2}{A_1}{\left[{{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\left[g_1\left(n\right)\right]}^2\left(1+\frac{n}{R}\right)dn\right]}^{-1}\\ g_1\left(n\right)={{\displaystyle \int }}_{-\frac{h}{2}}^z\left(1+\frac{n}{R}\right)ndn;A_1={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}\left(1+\frac{n}{R}\right)dn;I={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}\left(1+\frac{n}{R}\right)n^{2}dn\end{array}$$

(12)

$${\alpha }_n\left(h,R\right)=\frac{140R^3}{168R^3-140R^{2}h+34R{h}^2+7h^3}$$

(13)

The advantage of the new shear correction factor, Equation (13), is that it takes into account the slenderness ratio of the tunnel ventilation shaft in a geometrical optimization analysis (

Table 1. Comparison between the new shear correction factor for tunnel ventilation shafts, and the shear correction factor for plates for different slenderness ratios $\left(\frac{h}{R}\right)$ .

$\left(\frac{\mathit{h}}{\mathit{R}}\right)$	${\mathit{\alpha }}_{\mathit{n}}=\frac{5}{6}\left(\mathbf{Plate}\right)$	${\mathit{\alpha }}_{\mathit{n}}\left(\mathbf{Shafts}\right)$	Relative Error (%)
$0.05$	$0.83\overline{3}$	$0.869$	$4.142$
$0.10$	$0.83\overline{3}$	$0.907$	$8.159$
$0.12$	$0.83\overline{3}$	$0.923$	$9.751$
$0.15$	$0.83\overline{3}$	$0.947$	$12.038$
$0.20$	$0.83\overline{3}$	$0.990$	$15.859$
$0.22$	$0.83\overline{3}$	$1.008$	$17.361$
$0.25$	$0.83\overline{3}$	$1.035$	$19.517$

$\left({\alpha }_n\right)$ Shafts was assumed as a pattern.

).

3. New Mathematical Operational Model

For this study, the internal equilibrium equations generate an equation set of one hyperstatic degree (HD), Equation (14) [1]. The kinematics equations, Equation (15), (Mindlin–Reissner) in cylindrical coordinate associated to Equation (14) can be obtained through the principle of the virtual works (PVW) but also with the transposition method [1,10,11].

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i=1}^n}{F}_{\overrightarrow{e_z}}=0\Rightarrow \frac{d{N}_z}{dz}+q_{z_{PP}}=0,{\displaystyle {\displaystyle \sum }_{i=1}^n}{F}_{\overrightarrow{e_n}}=0\Rightarrow \frac{N_{\theta }}{R}-\frac{d{Q}_z}{dz}+q_n=0\\ {\displaystyle {\displaystyle \sum }_{i=1}^n}{M}_{\overrightarrow{e_{\theta }}}=0\Rightarrow \frac{d{M}_z}{dz}-Q_z=0,GH=1\left(N_z,N_{\theta },M_z,Q_z\right)\end{array}$$

(14)

$$\begin{array}{c}{\epsilon }_{\theta }=-\frac{U_n}{R},{\epsilon }_z=\frac{d{U}_z}{dz},{\kappa }_z=-\frac{d{\psi }_z}{dz},{\gamma }_{zn}=-\frac{d{U}_n}{dz}-{\psi }_z\\ DF=3\left(U_n,U_{z}y{\psi }_z\right)\end{array}$$

(15)

By substituting the kinematics equations, Equation (15), and the coupled constitutive equations, Equation (7), into the equilibrium equations, Equation (14), a new mathematical operational model, Equation (16), is obtained, whose unknowns are the middle surface displacement of the shell $\left(U_n,U_z,{\psi }_z\right)$. Equation (16) generalizes the mathematical model for the isotropic shell presented by Rotter and Sadowski [4].

The analytical solution of the differential set equations, Equation (16), for complex boundary conditions and the inclusion of sections (loads, stiffness, etc.) is inapproachable. However, if the new mathematical model is formulated in terms of the resultant internal forces per unit of arc longitude (inverse operational model), it is possible to obtain an analytical solution of the model.

$$\begin{array}{c}K\frac{d^2{U}_z}{d{z}^2}-\frac{\mu }{R}\frac{d{U}_n}{dz}-\frac{D_f}{R}\frac{d{\psi }_z}{dz}=-q_z\\ \frac{2\mu }{R{\alpha }_n\left(1-\mu \right)}\frac{d{U}_z}{dz}+\frac{d^2{U}_n}{d{z}^2}-\frac{2{\alpha }_1}{{\alpha }_n{R}^2\left(1-\mu \right)}{U}_n+\frac{d{\psi }_z}{dz}=-\frac{2q_n}{{\alpha }_{n}K\left(1-\mu \right)}\\ \frac{D_f}{R}\frac{d^2{U}_z}{d{z}^2}+{\alpha }_{n}K\frac{\left(1-\mu \right)}{2}\frac{d{U}_n}{dz}-D_f\frac{d^2{\psi }_z}{d{z}^2}+{\alpha }_{n}K\frac{\left(1-\mu \right)}{2}{\psi }_z=0\end{array}$$

(16)

The compatibility equations of deformations (Saint-Venant identities) and resultant internal forces functions per unit of arc longitude were obtained [12,13] for mathematical modeling in terms of the resultant internal forces per unit of arc longitude.

The new functions for the resultant internal forces per unit of arc longitude, Equation (18), and the new compatibility equation of deformations, Equation (17), were obtained by applying the indeterminate coefficient method (ICM) [14]. By substituting Equation (18) in the internal equilibrium equations of the shell, Equation (14), the equation satisfies itself, which is the general solution of the internal equilibrium of the shell.

$$R\frac{d^2{\epsilon }_{\theta }}{d{z}^2}+{\kappa }_Z-\frac{d{\gamma }_{ZN}}{dz}=0$$

(17)

$$\begin{array}{c}{N}_{\theta }\left(z\right)=-R\frac{d^2\theta }{d{z}^2}-R{q}_n\left(z\right),M_z\left(z\right)=-\theta ,Q_Z=-\frac{d\theta }{dz}\\ N_z={\gamma }_m\left(Hz-\frac{z^2}{2}-\frac{H^2}{2}\right)\end{array}$$

(18)

By substituting Equation (18) and the coupled constitutive model with curvature, Equation (7), in the compatibility equation of deformations, Equation (17), the new inverse-operational mathematical model is obtained, Equation (19).

$$\begin{array}{c}\underset{BENDING CONTRIBUTION}{\underbrace{\frac{d^4\theta }{d{z}^4}}}+\underset{SHEARFORCECONTRIBUTION}{\underbrace{A_1\frac{d^2\theta }{d{z}^2}}}+\underset{\begin{array}{c}MEMBRANALFORCE\\ CONTRIBUTION\end{array}}{\underbrace{B_1\theta }}=F_{CG}\left(z\right)\\ F_{CG}\left(z\right)=\frac{H_{21}}{R{H}_{22}}\frac{d^2{N}_z}{d{z}^2}-\frac{d^2{q}_N}{d{z}^2}+\frac{H_{31}}{R^2{H}_{22}}{N}_z-\frac{H_{32}}{R{H}_{22}}{q}_N\\ A_1=-\frac{1}{R{H}_{22}}\left(H_{23}+H_{32}-\frac{1}{R}\frac{\left(1+\mu \right)}{2Eh{\alpha }_n}\right),B_1=\frac{H_{33}}{R^2{H}_{22}}\end{array}$$

(19)

The coefficients A₁ and B₁ from the Equation (19) characterize the influence of the shear and the membrane force facing the bending respectively.

4. Particular Cases

If the constitutive curvature of the shell is disregarded in the mathematical model, Equation (19), the following mathematical model is obtained, Equation (20).

$$\begin{array}{c}\frac{d^4\theta }{d{z}^4}+A\frac{d^2\theta }{d{z}^2}+B\theta =F_{SG}\left(z\right);F_{SG}\left(z\right)=-\frac{\mu }{R}\frac{d^2{N}_z}{d{z}^2}-\frac{d^2{q}_N}{d{z}^2}\\ A=-\frac{12}{5}\frac{\left(1+\mu \right)}{R^2}yB=\frac{4}{C^4};C^4=\frac{h^2{R}^2}{3\left(1-{\mu }^2\right)}\end{array}$$

(20)

The mathematical model expressed in Equation (20), is a strong-operational formulation in terms of the resultant internal forces per unit of arc longitude of the Mindlin–Reissner flexural model for a coaxial cylinder; and its numerical resolution by FEM (weak-formulation) represents a reference for the analysis of thick shells.

Furthermore, Equation (21) shows a comparative analysis between the second member of the mathematical models, $F_{CG}\left(z\right)\mathrm{and}{F}_{SG}\left(z\right),$ with and without constitutive curvature respectively, Equations (19) and (20). Equation (21) shows the additional effects originated by the axial forces and the orthogonal load model that are not considered in Equation (20).

$$\begin{array}{cc}{F}_{CG}\left(z\right)=\underset{\begin{array}{c}Firstordereffect\\ F_{SG}\left(z\right)\end{array}}{\underbrace{\underset{\cong -\frac{\mu }{R}}{\underbrace{\frac{H_{21}}{R{H}_{22}}}}\frac{d^2{N}_z}{d{z}^2}-\frac{d^2{q}_N}{d{z}^2}}}+& \underset{ConstitutiveSecondordereffect}{\underbrace{\underset{Oftheaxialforce}{\underbrace{\frac{H_{31}}{R^2{H}_{22}}{N}_z}}-\underset{Ofthecircumferentialforce}{\underbrace{\frac{H_{32}}{R{H}_{22}}{q}_N}}}}\end{array}$$

(21)

If the shell constitutive curvature and the shear deformation [15] are disregarded $\left({\gamma }_{zn}=0\right)$ in the general mathematical model, Equation (19), the bending model in operational formulation of Love–Kirchhoff in term of the resultant internal forces per unit of arc longitude is obtained, Equation (22).

$${\gamma }_{ZN}=0\Rightarrow {\psi }_z=\frac{d{U}_n}{dz}\Rightarrow \frac{d^4\theta }{d{z}^4}+\frac{4}{C^4}\theta =-\frac{\mu }{R}\frac{d^2{N}_z}{d{z}^2}-\frac{d^2{q}_N}{d{z}^2}$$

(22)

Since the dependency between the twist, ${\psi }_z,$ and the fundamental displacement, $U_n,$ are conditioned, the Love–Kirchhoff model for the bending of a coaxial cylinder generates an equality between the hyperstatic degree $\left(\theta =1\right)$ and the degree of freedom, $U_n,$ of the model. This means that the formulation in terms of the resultant internal forces per unit of arc longitude, Equation (22), and the formulation in terms of fundamental displacement, Equation (23), endure to only one differential equation.

Equation (23) is the most widely used mathematical model for analytical analysis and structural flexo-compression design, whose scheme of analysis is reduced to a coaxial cylinder (Figure 1) [1,2,16,17,18,19,20].

$$\begin{array}{c}\boxed{\frac{d^4{U}_n}{d{z}^4}+\frac{4}{C^4}{U}_n=\frac{q_N}{D_f}+\frac{\mu }{R{D}_f}{N}_z}{C}^4=\frac{h^2{R}^2}{3\left(1-{\mu }^2\right)};\\ S.F.S=\left\{e^{\left(\frac{z}{c}\right)}cos\left(\frac{z}{c}\right)z,e^{\left(\frac{z}{c}\right)}sin\left(\frac{z}{c}\right)z+e^{\left(-\frac{z}{c}\right)}cos\left(\frac{z}{c}\right)z,e^{\left(-\frac{z}{c}\right)}sin\left(\frac{z}{c}\right)z\right\}\end{array}$$

(23)

5. Analytical Resolution of the New Mathematical Model

The mathematical model, with and without constitutive curvature in terms of the resultant internal forces per unit of arc longitude, Equations (19) and (20), is the nonhomogeneous fourth-order-linear differential equations with constant coefficients. The main difficulty is to obtain its homogeneous solution, $\theta c,$ due to the shear of each model (terms A and A₁). The general solution of the model is the sum of two solutions due to its linearity, Equation (24), i.e., the homogeneous solution, $\theta c,$ and the particular solution, $\theta p$ [21,22].

By applying the general theory of differential equation [21], the homogeneous solution $\theta c\left(z\right)$, Equations (27) and (28), is obtained, and the characteristic equation, Equation (25), and the fundamental solution system (F.S.S), Equation (26), are obtained. Out of the two F.S.S, the most employed in engineering is the case $D<0$ due to the fact that the thickness is smaller in comparison with the principal radius of curvature $\left(\frac{h}{R}\ll 1\right)$.

$$\theta \left(z\right)=\theta c\left(z\right)+\theta p\left(z\right),$$

(24)

$$\begin{array}{c}{m}^4-A{m}^2+B=0\Rightarrow m_1=\sqrt{\alpha },m_2=\sqrt{\beta },m_3=-\sqrt{\alpha },m_4=-\sqrt{\beta }\\ \mathrm{Being}:D=A^2-4B,\alpha =\frac{A}{2}+\frac{\sqrt{D}}{2}y\beta =\frac{A}{2}-\frac{\sqrt{D}}{2}\end{array}$$

(25)

$$\begin{array}{c}IfD>0\Rightarrow \left(\mathrm{real}\mathrm{solutions}\right)\Rightarrow S.F.S=\left\{e^{\sqrt{\alpha }z},e^{-\sqrt{\alpha }z},e^{\sqrt{\beta }z},e^{-\sqrt{\beta }z}\right\}\\ IfD<0\Rightarrow \left(\mathrm{complex}\mathrm{solutions}\right)\Rightarrow S.F.S=\begin{array}{c}\{e^{{\lambda }_{1}z}cos{\lambda }_{2}z,e^{{\lambda }_{1}z}sin{\lambda }_{2}z+\\ e^{-{\lambda }_{1}z}cos{\lambda }_{2}z,e^{-{\lambda }_{1}z}sin{\lambda }_{2}z\}\end{array}\end{array}$$

(26)

$$\theta c\left(z\right)=A_1{e}^{{\lambda }_{1}z}\cos {\lambda }_{2}z+A_2{e}^{{\lambda }_{1}z}\sin {\lambda }_{2}z+A_3{e}^{-{\lambda }_{1}z}\cos {\lambda }_{2}z+A_4{e}^{-{\lambda }_{1}z}\sin {\lambda }_{2}z$$

(27)

$$\mathrm{Being}:\begin{array}{c}{\lambda }_1=\sqrt{r}cos\frac{\theta }{2}\\ {\lambda }_2=\sqrt{r}sin\frac{\theta }{2}\end{array}\}\begin{array}{c}r=\frac{1}{2}\sqrt{A^2+\left|D\right|}\\ \theta =Ta{n}^{-1}\left(\frac{\sqrt{\left|D\right|}}{A}\right)\end{array}$$

(28)

If the mathematical condition ${\lambda }_1={\lambda }_2=\frac{1}{c}$ is evaluated in the new F.S.S, Equation (26), then the general solution of the Love–Kirchhoff model with constitutive isotropy is obtained, Equation (23).

The particular or nonhomogeneous solutions $\left(\theta p_1\mathrm{and}\theta p_2\right)$, Equation (30), are obtained by applying the Lagrange’s parameters variation method [21], and depend on the function that operates on the left-hand side term of the differential equations, $F_{CG}\left(z\right)\mathrm{and}{F}_{SG}\left(z\right)$.

Figure 4 and Equation (29) shows the results when a load model characterized by the own weight of the tunnel ventilation shaft wall, $q_{z_{pp}}\left(z\right),$ the horizontal effective stress, $K_0{\sigma }^{\prime },$ and the overload action, $q_{n_2},$ are defined.’

$$\begin{array}{c}{q}_{z_{pp}}\left(z\right)=-{\gamma }_{m}h\\ \begin{array}{c}{q}_n\left(z\right)=q_{m}z+q_{n_1};q_m=\frac{q_{n_2}-q_{n_1}}{H};q_{n_1}=K_{0}H\left({\gamma }_s-{\gamma }_a\right)+{\gamma }_{a}H+q_{n_2}\end{array}\end{array}$$

(29)

$$\begin{array}{c}\theta p_1=f_{1}z+f_2\Rightarrow ForEquation\left(18\right)\\ \theta p_2=0\Rightarrow ForEquation\left(19\right)\\ Being:f_1=\frac{-H_{32}{q}_m}{B_{1}R{H}_{22}}+\frac{H_{31}{\gamma }_{m}h}{B_1{R}^2{H}_{22}};f_2=-\frac{H_{32}{q}_{n_1}}{B_{1}R{H}_{22}}-\frac{H_{31}{\gamma }_{m}hH}{B_1{R}^2{H}_{22}}\end{array}$$

(30)

Once the general solution, Equations (27) and (30), of the resultant internal forces per unit of arc longitude (θ) is obtained, it is necessary to evaluate the boundary conditions imposed on the mathematical model (

Table 2. Comparison between the new shear correction factor for tunnel ventilation shafts, and the shear correction factor for plates for different slenderness ratio $\left(\frac{h}{R}\right)$ .

ID	Simple Term	$\mathbf{Radial}\mathbf{Displacement}\left({\mathit{U}}_{\mathit{n}}\right)$ $\mathbf{Twist}\mathbf{Angle}\left({\mathit{\psi }}_{\mathit{z}}\right)$ $\mathbf{Shear}\left({\mathit{Q}}_{\mathit{z}}\right)$ $\mathbf{Bending}\mathbf{Moment}\left({\mathit{M}}_{\mathit{z}}\right)$	$\mathbf{Axial}\mathbf{Displacement}\left({\mathit{U}}_{\mathit{z}}\right)$
BC1 r	Clamped	$U_n=0$ and ${\psi }_z=0$	$U_z=0$
BC1 f		$U_n=0$ and $M_z=0$	$U_z=0$
BC2 r	Pinned	$U_n=0$ and ${\psi }_z=0$	$U_z\ne 0$
BC2 f		$U_n=0$ and $M_z=0$	$U_z\ne 0$
BC3	Free edge	$M_z=0$ and $Q_z=0$	$U_z\ne 0$

) [23]. Six-boundary conditions are required to solve the flexo-compression problem with coupled constitutive curvature (Table 2). Boundary conditions used in EN 1993-1-6 Eurocode [23] were adopted in this paper, along with the Mindlin–Reissner and Timoshenko bending theory particularity.

From a complete and fit definition of the function, $\theta ,$ the internal forces per unit of arc longitude $(M_Z,M_{\theta },Q_z,N_{\theta },N_z$) are founded by applying Equation (18). Subsequently, by using the coupled constitutive equations, Equation (10), the deformations on the middle surface of the shell $\left({\kappa }_z,{\kappa }_{\theta },{\gamma }_{zN},{\epsilon }_{\theta },{\epsilon }_z\right)$ were obtained, and to conclude, the integration process, the displacement field $\left(U_n,U_z,{\psi }_z\right)$ is determined by employing the Mindlin–Reissner kinematics equations, Equation (15). The stress–strain state of the coaxial cylinder is finally obtained through Equations (3) and (4).

The analytical solutions were implemented in the mathematical assistant MATHLAB [24], developing the user graphic interfaces with GUIDE tool. A computer program was created for the stress–strain state analysis of tunnel ventilation shafts with and without coupled constitutive curvature (Figure 5).

6. Numerical Results and Discussions

In this section, fundamental concepts in the shells-analytical and numerical analyses (FEM) in general bending state are illustrated. The physical–mechanical parameters (

Table 3. Physical, mechanical and geometrical properties of the tunnel ventilation shaft N°1 of the project-“Río La Compañia” from Espejel.

Physical and Mechanical Properties

$E=2.378\times {10}^7\frac{kN}{m^2}$${\gamma }_a=10\frac{kN}{m^3}$${\gamma }_m=20.46\frac{kN}{m^3}$${\gamma }_s=10.23\frac{kN}{m^3}$$\begin{array}{c}\mu =0.2\\ K_0=0.8\end{array}$$q_{n_1}=206.735\frac{kN}{m^2}$$q_{n_2}=0$

Geometrical Properties

$h=0.70m$$H=20.30m$$R=6m$$\mathrm{Ratio}\left(\frac{h}{R}\right)=0.12$

) of a typical tunnel ventilation shaft built in the soft soil of Mexico City from Espejel [25] were employed for numerical comparison.

Bearing in mind that there is a continuity between the base and the tunnel ventilation shaft wall, a clamped condition on the cylinder base (BC1r in z = 0) was assumed, along with the free edge on the top (BC1r in z = H). The analytical results (with and without constitutive curvature) were compared with FEM results, using the ABAQUS software [26] (Figure 6 and Figure 7) and considering different formulations (

Table 4. Finite element types and statistical meshing.

Finite Element Type Description		$\frac{\mathit{h}}{\mathit{R}}=0.05$		$\frac{\mathit{h}}{\mathit{R}}=0.12$		$\frac{\mathit{h}}{\mathit{R}}=0.25$
		Nodes	Elements	Nodes	Elements	Nodes	Elements
Solid elements 3D	C3D20R: A 20-node quadratic brick, reduced integration	59,644	8432	289,708	57,528	326,636	71,928
	C3D8R: An 8-node linear brick, reduced integration, hourglass control	17,112	8432	77,456	57,528	84,952	71,928
Shell elements	S8R: An 8-node doubly curved thick shell, reduced integration	8694	8568	8694	8568	8694	8568
	S4R: A 4-node doubly curved thin or thick shell, reduced integration, hourglass control, finite membrane strains	25,956	8568	25,956	8568	25,956	8568

).

It is known, [7,27,28] that the standard 8-nodes-solid element with complete integration (C3D8), presents a locking-shear problem when used for shells modeling. This problem is bigger when the slenderness increases, and its bending behavior prevails. This numerical locking indicates the incapacity of the interpolation functions (and its gradients) to adapt to solid behavior that often invalidates the obtained solutions [28]. The elements of reduced integration of first and second order (C3D8R and C3D20R), present more precision in flexo-compression problems with bending predominance than the ones corresponding to complete integration-elements, in addition to a reduction in the computing time [26]. The second-order solid elements allow modeling of the curved surfaces with fewer elements, and a greater approximation in areas of high stresses concentration. The shell element S4R is valid for thin or thick-shell problems [26], but it has the disadvantage that it does not obtain the shear per unit of thickness ($Q_z$) as output, which allows the subsequent calculation of the transversal stress. The S8R element is valid for thick shell analysis.

The numerical models used in this paper considered a different slenderness ratio in the shell, and the employment of the shell elements (S4R and S8R) and solid elements (C3D8R and C3D20R) (Table 4 and Figure 6).

The analytic solution with constitutive curvature shows numerical superiority compared with shell-elements (S8S and S4R) for all results in the entire high of the tunnel ventilation shaft (Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13;

Table 5. Maximum circumferential stress in extrados for different slenderness ratio.

$\left(\frac{\mathit{h}}{\mathit{R}}\right)$	Length (m)	$\mathbf{Maximum}\mathbf{Circumferential}\mathbf{Stress}\mathbf{in}\mathbf{Extrados}\left({\mathit{\sigma }}_{\mathit{\theta }}\right)\left(\raisebox{1ex}{$\mathit{k}\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{m}}^{\mathbf{2}}$}\right.\right)$
		Analytical Solution			FEM Shells				FEM 3D
		With CC	Without CC	% EA	S4R	% ES4R	S8R	% ES8R	C3D8R	% EC3D8R	C3D20R	% EA-C3D20R
0.05	2.57	3833.002	3923.135	2.35	3596.92	6.16	3773.11	1.59	3730.3	0.514	3749.56	2.225
0.12	3.73	1523.094	1601.1821	5.13	1447.13	4.99	1447.63	5.21	1640.83	0.787	1628.02	6.445
0.25	5.1	644.19	706.634	9.69	586.868	8.90	587.301	9.69	752.557	0.557	748.391	13.923
	Where:	CC: Constitutive curvature; EA: relative error between the analytical results (pattern: analytical solution with CC); ES4R: relative error between the analytical result with CC and shell element S4R (pattern: analytical solution with CC); ES8R: relative error between the analytical result with CC and shell element S8R (pattern: analytical solution with CC); EC3D8R: relative error between the elements C3D8R and C3D20R (pattern: C3D20R); EA-C3D20R: relative error between the analytical result with CC and the element C3D20R (pattern: C3D20R).

,

Table 6. Maximum circumferential stress in intrados for different slenderness.

$\left(\frac{\mathit{h}}{\mathit{R}}\right)$	Length (m)	$\mathbf{Maximum}\mathbf{Circumferential}\mathbf{Stress}\mathbf{in}\mathbf{Intrados}\left({\mathit{\sigma }}_{\mathit{\theta }}\right)\left(\raisebox{1ex}{$\mathit{k}\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{m}}^{\mathbf{2}}$}\right.\right)$
		Analytical Solution			FEM Shells				FEM 3D
		With CC	Without CC	% EA	S4R	% ES4R	S8R	% ES8R	C3D8R	% EC3D8R	C3D20R	% EA-C3D20R
0.05	3.06	3770.237	3673.337	2.57	3834.26	1.70	3599.85	4.52	3619.5	4.517	3790.72	0.543
0.12	4.38	1560.739	1471.875	5.69	1335.59	14.43	1336.09	14.39	1665.84	0.471	1673.72	7.239
0.25	5.85	712.314	634.328	10.95	649.409	8.83	649.813	8.77	828.22	1.459	840.48	17.993
See Table 5 for legend

,

Table 7. Maximum longitudinal stress in extrados for different slenderness ratio.

$\left(\frac{\mathit{h}}{\mathit{R}}\right)$	$\mathbf{Maximum}\mathbf{Longitudinal}\mathbf{Stress}\mathbf{in}\mathbf{Extrados}\left({\mathit{\sigma }}_{\mathit{z}}\right)\left(\raisebox{1ex}{$\mathit{k}\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{m}}^2$}\right.\right)\mathbf{z}=0$
	Analytical Solution			FEM Shells				FEM 3D
	With CC	Without CC	% EA	S4R	% ES4R	S8R	% ES8R	C3D8R	% EC3D8R	C3D20R	% EA-C3D20R
0.05	6500.02	6556.31	0.87	4429.12	31.860	6006.91	7.586	412.305	93.334	6184.77	5.097
0.12	2452.9	2506.165	2.17	1623.46	33.815	1917.51	21.827	1265.24	53.086	2696.94	9.049
0.25	861.4145	910.11	5.65	432.582	49.782	537.565	37.595	535.274	57.637	1263.55	31.826
See Table 5 for legend

and

Table 8. Maximum longitudinal stress in intrados for different slenderness ratio.

$\left(\frac{\mathit{h}}{\mathit{R}}\right)$	$\mathbf{Maximum}\mathbf{Longitudinal}\mathbf{Stress}\mathbf{in}\mathbf{Intrados}\left({\mathit{\sigma }}_{\mathit{z}}\right)\left(\raisebox{1ex}{$\mathit{k}\mathit{N}$}\!\left/ \!\raisebox{-1ex}{${\mathit{m}}^{\mathbf{2}}$}\right.\right)\mathbf{z}=\mathbf{0}$
	Analytical Solution			FEM Shells				FEM 3D
	With CC	Without CC	% EA	S4R	% ES4R	S8R	% ES8R	C3D8R	% EC3D8R	C3D20R	% EA-C3D20R
0.05	7444.063	7386.986	0.77	5253.69	29.424	6837.6	8.147	412.305	94.214	7126.37	4.458
0.12	3390.943	3336.8407	1.60	2450.06	27.747	2748.18	18.955	2179.69	40.199	3644.93	6.968
0.25	1790.302	1740.786	2.77	1258.13	29.725	1368.24	23.575	1444.26	33.889	2184.61	18.049
See Table 5 for legend

).

The C3D8R element shows a bad convergence in the obtaining of the maximum longitudinal stress at the base of the cylinder (clamped) (bending concentration) (Table 7 and Table 8). Although, the elements are provide by the Hourglass control [26] there is a stiffness of the shell in that section (Figure 10 and Figure 11). Due to the linearity of the interpolation functions, it is necessary to generate more than one C3D8R element inside the thickness of the shell, which will be able to model the bending concentration in that cross section. The C3D20R element employs trigonometric interpolation functions, which allow modeling with one element into the small thickness of the shell $\left(\frac{h}{R}\le 0.05\right)$ and the stress distribution (circumferential and longitudinal) in the alteration zones [1,3].

For the longitudinal stress (longitudinal flexor-compression), in all shell domains, the analytical solution with constitutive curvature shows numerical superiority than the one obtained through the C3D8R element with more than one element inside the thickness of the shell (Figure 10 and Figure 11; Table 7 and Table 8).

An important aspect for structural design in reinforced concrete is the change in the circumferential flexural moment configuration with the increment of the slenderness ratio of the shell (Figure 13) for analytical solution with constitutive curvature. For a slenderness ratio of 0.25, the studied shell did not have an inversion in the circumferential flexural moment (Figure 13), exalting the physical phenomena that occur: a decreasing of the circumferential negative moment when increasing the slenderness ratio. That is, when the shell is not thin anymore $\left(\frac{h}{R}>0.05\right)$, the bending moment analysis in the circumferential direction is insufficient if the classical constitutive relationship is applied, Equation (31), [1,4,5] on a revolution shell (isotropic, orthotropic and anisotropic) under a general distribution of axial-symmetric pressures.

$$M_{\theta }=\pm \mu M_z,$$

(31)

7. Conclusions

The formulation and analytic solution presented in this paper displays the generality of combining the Mindlin–Reissner theory, the shear correction factor properly for cylindrical shell and the general coupled constitutive model (isotropic and orthotropic) without disregarding the shell curvature from the constitutive.

The proposed mathematical model presents three independent degrees of freedom $\left(U_n,U_z,{\psi }_z\right)$ and it is formulated in terms of the resultant internal forces per unit of arc longitude. This has the advantage of allowing the analysis of cylindrical shell by analytical formulation with thin, medium and thick thickness under a general distribution of axial-symmetric pressures with the very best numerical reliability in the entire shell domain.

Potentially engineering applications include tunnel ventilation shafts, tubular piles under earth pressures, tanks under hydrostatic pressures, reinforced concrete silos under granular solid pressures, gas-filled cisterns, simple alteration effect analysis in thin and medium shells and chimneys. A practical example has been illustrated in a tunnel ventilation shaft analysis for different slenderness ratio $\left(\frac{h}{R}\right)$, resulting in the following significant conclusions:

• When the isotropic shell is thin, $\left(\frac{h}{R}\le 0.05\right)$ its predominant resistant mechanism is the circumferential membrane force with an inversion of bending moments in the both main directions. From an increase in the slenderness ratio, the flexural contribution in the two main directions dominated the internal equilibrium of the shell and the inserting of the constitutive curvature acquires the biggest importance in the structural response of the shell. For analysis and tunnel ventilation shafts design with slenderness ratio of $\frac{h}{R}\ge 0.12$ it is recommendable to insert the constitutive curvature;
• The equations and the general methodology displayed in this paper might be usefully employed in the analysis and design of the cylindrical shell (isotropic and orthotropic) under general distribution of axial-symmetric pressures. The mathematical model formulated in terms of the internal forces per unit of arc longitude allows to solve differential equation systems of multiple degrees of freedom and to model the complex boundary conditions by the Saint-Venant simplification. For this study case, the equations can be applied by means of the basic spreadsheets as tools to assist in the design.