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

: 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 ( α n ) in terms of the thickness ( h ) and the radius of curvature ( R ) . 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 ﬁrst one, between the new model and the Mindlin–Reissner model without constitutive curvature with the shear correction factor α n = 5/6 as a constant, and the other, between the new model and the tridimensional numerical models (solids and shells) obtained by ﬁnite 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.


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. Nonexistence of the constitutive coupling between the membrane forces and the bending moments for any slenderness ratio h R ; 3.
Shear correction factor (α n ) is employed in the constant shear force, independently of the shell slenderness ratio.
Resultant internal forces per unit of arc longitude and middle deformations of the shell ( , , , , ) 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]  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 (α 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).

Coupled Constitutive Model for General Bending of a Coaxial Cylinder
The term ds * , represents the differential arc increased in the shell differential arc (dθ) of the middle surface ( Figure 3). Figure 3 shows the stress state in a "differential point" that belong to the surface S * (ds * dz). point" that belong to the surface * ( * ).
From the stress-strain state definition for the "differential point", which belong to the surface * ( * ) , Equations (3) and (4), it is possible to integrate the surfaces ( and ) in order to obtain the resultant internal forces per unit of arc longitude: bending moments ( , ), shear and membrane axial forces ( , ), Equation (5) (7) and (8). For the integration process, the term , Equation (6), was developed in n potentialfunction (Taylor series expansion) [3] affecting the integrand with the substitution of the stress-strain state for the "differential point", Equations (3)-(5).  From the stress-strain state definition for the "differential point", which belong to the surface S * (ds * dz), Equations (3) and (4), it is possible to integrate the surfaces (dA 1 and dA 2 ) in order to obtain the resultant internal forces per unit of arc longitude: bending moments (M Z , M θ ), shear and membrane axial forces (N θ , N z ), 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 ( 1 1+ n R ), Equation (6), was developed in n potentialfunction (Taylor series expansion) [3] affecting the integrand with the substitution of the stress-strain state for the "differential point", Equations (3)- (5).
γ zn σ θ * ndn, α n = 1/k n f : f orces unit; l : length unit 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 1 R = 0 is neglected in Equations (7) or (8), the classical shell model, Equation (2), is obtained. In addition, the linear-elastic coupled constitutive model, in  (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 [D C ], 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 The shear correction factor (α 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 α n = 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).
By inserting in Equation (11) the shell curvature geometrical relations (Figure 3), an expression for the shear correction factor variable (α n ) per unit of arc longitude in crosscylindrical and homogeneous sections is obtained, Equations (12) and (13), which depends on the thickness and the radius of curvature of the shell (Figure 3).
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).

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].
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 (U n , U z , ψ z ). 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.
The compatibility equations of deformations (Saint-Venant identities) and resultant internal forces functions per unit of arc longitude were obtained [12,13]  plying 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.
By substituting Equation (18) and the coupled constitutive model with curvature, Equation (7), in the compatibility equation of deformations, Equation (17), the new inverseoperational mathematical model is obtained, Equation (19).
The coefficients A 1 and B 1 from the Equation (19) characterize the influence of the shear and the membrane force facing the bending respectively.

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).
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 (weakformulation) 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 (z) and F SG (z), 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).
First order e f f ect F SG (z) If the shell constitutive curvature and the shear deformation [15] are disregarded (γ zn = 0) 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).
Since the dependency between the twist, ψ 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 (θ = 1) 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.

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, θc, due to the shear of each model (terms A and A 1 ). The general solution of the model is the sum of two solutions due to its linearity, Equation (24), i.e., the homogeneous solution, θc, and the particular solution, θ p [21,22].
By applying the general theory of differential equation [21], the homogeneous solution θc(z), 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 h R 1 .
Being : If the mathematical condition λ 1 = λ 2 = 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 (θ p 1 and θ p 2 ), 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 (z) and F SG (z). 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 (z), the horizontal effective stress, K 0 σ , and the overload action, q n 2 , are defined.' Once the general solution, Equations (27) and (30), of the resultant inte unit of arc longitude (θ) is obtained, it is necessary to evaluate the bound imposed on the mathematical model (Table 2) [23]. Six-boundary conditio to solve the flexo-compression problem with coupled constitutive curva Boundary conditions used in EN 1993-1-6 Eurocode [23] were adopted in th with the Mindlin-Reissner and Timoshenko bending theory particularity.
From a complete and fit definition of the function, , the internal for arc longitude ( , , , , ) are founded by applying Equation (18). Su using the coupled constitutive equations, Equation (10), the deformations surface of the shell ( , , , , ) were obtained, and to conclude, the i cess, the displacement field ( , , ) is determined by employing the ner kinematics equations, Equation (15). The stress-strain state of the coa finally obtained through Equations (3) and (4). 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) [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. Table 2. Comparison between the new shear correction factor for tunnel ventilation shafts, and the shear correction factor for plates for different slenderness ratio h R .

ID Simple Term
From a complete and fit definition of the function, θ, the internal forces per unit of arc longitude (M Z , M θ , Q z , N θ , 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 (κ z , κ θ , γ zN , ε θ , ε z ) were obtained, and to conclude, the integration process, the displacement field (U n , U z , ψ z ) 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 MATH-LAB [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).

R PEER REVIEW 10 of 18
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).

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) of a typical tunnel ventilation shaft built in the soft soil of Mexico City from Espejel [25] were employed for numerical comparison.
, mechanical and geometrical properties of the tunnel ventilation shaft N°1 of the project-"Río La Comjel. 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] (Figures 6 and 7) and considering different formulations (Table 4).

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) of a typical tunnel ventilation shaft built in the soft soil of Mexico City from Espejel [25] were employed for numerical comparison. 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
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] (Figures 6 and 7) and considering different formulations (Table 4).  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 ( ) 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 (Figures 8-13; Tables 5-8).   It is known, [7,27,28] that the standard 8-nodes-solid element with complete integ tion (C3D8), presents a locking-shear problem when used for shells modeling. This pro lem is bigger when the slenderness increases, and its bending behavior prevails. This n merical locking indicates the incapacity of the interpolation functions (and its gradien to adapt to solid behavior that often invalidates the obtained solutions [28]. The eleme of reduced integration of first and second order (C3D8R and C3D20R), present more p cision in flexo-compression problems with bending predominance than the ones cor sponding to complete integration-elements, in addition to a reduction in the computi time [26]. The second-order solid elements allow modeling of the curved surfaces w fewer elements, and a greater approximation in areas of high stresses concentration. T shell element S4R is valid for thin or thick-shell problems [26], but it has the disadvanta that it does not obtain the shear per unit of thickness ( ) as output, which allows t subsequent calculation of the transversal stress. The S8R element is valid for thick sh analysis.
The numerical models used in this paper considered a different slenderness ratio the shell, and the employment of the shell elements (S4R and S8R) and solid eleme (C3D8R and C3D20R) ( Table 4 and Figure 6).
The analytic solution with constitutive curvature shows numerical superiority co pared with shell-elements (S8S and S4R) for all results in the entire high of the tunn ventilation shaft (Figures 8-13; Tables 5-8).  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 (Figures 8-13; Tables 5-8).

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 ( , , ) 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 axialsymmetric 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 , resulting in the following significant conclusions: 1. When the isotropic shell is thin, ≤ 0.05 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 ≥ 0.12 it is recommendable to insert the constitutive curvature;  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 5 for legend 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) (Tables 7  and 8). Although, the elements are provide by the Hourglass control [26] there is a stiffness of the shell in that section (Figures 10 and 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 h R ≤ 0.05 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 (Figures 10 and 11; Tables 7 and 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 h R > 0.05 , 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.

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 (U n , U z , ψ z ) 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 axialsymmetric 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 h R , resulting in the following significant conclusions:

1.
When the isotropic shell is thin, h R ≤ 0.05 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 h R ≥ 0.12 it is recommendable to insert the constitutive curvature; 2.
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. Data Availability Statement: Not report any data.