A shell element for design‐oriented elasto‐plastic analysis of reinforced concrete wall structures using convex optimization

The iterative nature of the design processes for building structures requires computational models to be robust, efficient, and accessible while reflecting the actual structural behavior with sufficient accuracy. For limit state analysis of reinforced concrete structures, efficient linear‐elastic models are generally inaccurate, while the modeling and computational complexity of most high‐accuracy nonlinear models inhibit their use in design processes for large‐scale structures. A recently proposed framework for elasto‐plastic analysis of cracked reinforced concrete panels was demonstrated to be capable of analyzing models with more than 10,000 finite elements within minutes on a standard personal computer. This paper proposes an extension of this work in terms of a finite shell element for elasto‐plastic analysis of fully cracked reinforced concrete wall structures subjected to monotonic loading. Using nonlinear‐elastic constitutive models to imitate elasto‐plasticity, the proposed shell element couples the in‐plane section force variation to the nonlinear through‐thickness stress variation in a stress‐based formulation using a layer‐based submodel. The method is validated with exact solutions, and its applicability in design processes involving large‐scale structures is demonstrated using a finite element model of a four‐story stairwell with more than three million variables, which is solved within minutes on a personal computer.


| INTRODUCTION
The design of building structures takes place in the interface between a multitude of parties (including owners/ clients, architects, and engineers from different disciplines) where the design requirements are continuously Discussion on this paper must be submitted within two months of the print publication. The discussion will then be published in print, along with the authors' closure, if any, approximately nine months after the print publication. updated. This makes the process inherently iterative, requiring engineers to adopt tools for carrying out the various limit state analyses not only based on their accuracy but just as much their robustness, computational efficiency, and modeling complexity. As a result, traditional linear-elastic finite element (FE) models are still among the most widely used numerical tools by structural engineers. This is also the case for structures of reinforced concrete (RC) despite the behavior of cracked concrete being highly anisotropic and nonlinear.
State-of-the-art computational software for reinforced concrete based on fracture mechanics and elastoplasticity (see, e.g., reference [1] or [2]) can simulate the nonlinear behavior of cracked concrete rather accurately. However, the application of these types of models to large-scale structures is generally associated with a substantial modeling and computational complexity as well as the risk of convergence issues. In most practical applications, reinforced concrete structures are deliberately designed to obtain a ductile behavior with a stable crack development and failure governed by yielding of the reinforcement, for example, through minimum requirements to the reinforcement. 3 In these cases, the modeling and computational complexity associated with limit state verification can be significantly reduced, for example, by assuming quasi-static loading and fully cracked concrete.
For ultimate limit state (ULS) analysis of appropriately reinforced concrete structures, finite element limit analysis (FELA) has proven to be a reliable and computationally efficient method for a variety of types of RC structures. [4][5][6][7][8] By adopting rigid-plastic material models for the reinforcement and the concrete, the method determines the ultimate capacity and collapse mode of RC structures with great computational robustness and efficiency owing to the use of convex optimization algorithms. However, due to the rigid-plastic material models, FELA cannot provide information on finite deformations, ruling out its use in serviceability limit state (SLS) analysis or buckling analysis.
To resolve this issue, a framework for in-plane elastoplastic analysis of cracked RC panels, also based on convex optimization, was recently proposed by Vestergaard et al. 9 Being intended for application in limit state analysis of appropriately reinforced structures, the proposed method uses bilinear-elastic constitutive models for the reinforcement and the concrete to imitate elastoplasticity with isotropic hardening for quasi-static loading. In order to capture the softening effect of concrete in compression under transverse cracking, an iterative reduction of the concrete compression strength with the effectiveness factor was incorporated in the model in accordance with previous works on plasticity-based analysis. 10 With this approach, the proposed method applies the principle of minimum complementary energy to cast the stress-field problem as a convex optimization problem. This approach was demonstrated to allow the determination of finite deformation parameters, providing the basis for evaluating crack widths and structural displacements in SLS as well as assessing reinforcement ductility in ULS. With the use of convex solvers, the method was demonstrated to allow the analysis of a structure discretized into more than 13,000 finite elements and effectively more than one million variables on a personal computer within minutes.
In order to extend the framework to cover cases with plate bending, for example, walls subjected to transverse load or with adjoining wall corners, as well as enable a subsequent extension of the framework to cover buckling analysis in future works, this paper proposes a stressbased finite shell element for elasto-plastic analysis of RC wall structures. The proposed shell element couples the in-plane variation of the section forces and complementary energy to a submodel for the through-thickness stress variation in a discrete number of points. In the submodel, the through-thickness stress variation is modeled using a layer-based approach as proposed by several other authors (e.g., Larsen 5 and Jensen 11 in the context of FELA, and Backes 12 in the context of elastic-plastic stress fields). Assuming a plane state of stress in the layers, the same constitutive model as proposed in reference [9] is used, that is, bilinear-elasticity for reinforcement steel and fully cracked concrete, and with an iterative reduction of the concrete compressive strength to account for the effect of transverse cracking.
The paper is structured as follows: first, a brief description is given for the principles of variational methods and convex optimization employed in the paper. After this, the basis of the proposed model is described starting with the constitutive model, moving on to the layer submodel, and concluding with the equilibriumtype finite element. Finally, the model is verified, and its capabilities are demonstrated using four numerical examples of increasing complexity, dealing with a single crosssection, a beam, a slab, and finally, a four-story stairwell with openings.

| VARIATIONAL METHODS AND CONVEX OPTIMIZATION
For an elastic body with the domain Ω and the boundary Γ = {Γ s , Γ t } which is subjected to the prescribed displacements u on Γ s , the prescribed traction t on Γ t , and the prescribed body forces f in Ω, the governing equations are where σ and ε are the stress and strain vectors, u is the displacement vector, T is the boundary stress-to-traction transformation matrix, and ∂ is the differential operator. In some simple cases, an exact solution to this set of equations may be found analytically, while most other cases demand the use of approximate solution methods, for example, FE methods. A common approach for approximating the governing equations in FE analysis is based on the use of certain variational principles where some of the Equations (1) to (5) are enforced implicitly by requiring stationarity of a functional; the most widely used principle being the principle of minimum potential energy, 13 which approximates the equations of equilibrium, (1) and (5), by requiring stationarity of the potential energy functional, Π p (ε, u). This approach yields the equations for the displacement-based FE analysis. A less used principle is the principle of minimum complementary energy 14 which approximates the strain-displacement relation and the kinematic boundary conditions, (3) and (4), by requiring stationarity of the complementary energy functional, Π c (σ, t), leading to the equations for stress-based FE analysis. While the two principles, in the limit, converge towards the same exact solution, they do so from different directions meaning that displacement-based FE solutions generally exhibit decreased stiffness with refined discretization, whereas the stiffness of stress-based FE solutions generally increases with refined discretization. This property, combined with the compatibility between the principle of minimum complementary energy and the lower bound theorem of plasticity (i.e., that a safe and statically admissible stress field cannot cause collapse), means that stress-based FE analysis, in many applications, can be seen as leading to conservative results for coarse mesh discretizations. The complementary energy functional is defined as where is the complementary elastic energy density which is popularly interpreted as the area to the left of the stressstrain curve. With the definition (6), the principle of minimum complementary energy can be stated as the optimization problem: If C 0 is convex, this problem is a convex optimization problem (see, e.g., reference [15] for a general treatment of convex optimization), meaning that it can be solved using efficient and robust commercial algorithms. Furthermore, the kinematics of the problem can be interpreted directly from the dual solution, which is found simultaneously by so-called primal-dual interior-point algorithms. 16 In this paper, the optimization problem (8) is reformulated as a second-order cone program and solved using the commercial conic solver MOSEK. 17

| CONSTITUTIVE MODEL
The framework proposed in this paper is oriented towards use in design processes involving orthogonally reinforced concrete panels whose thickness is much smaller than their planar dimensions. Due to the constraints on modeling and computational complexity that are implicitly present in these iterative processes, the uniaxial behavior of the reinforcement and the concrete are taken as piecewise-linear elastic, as illustrated in Figure 1, in order to imitate elasto-plastic behavior for quasi-static loading. As illustrated, zero stiffness is assumed for the concrete in tension (corresponding to fully cracked concrete), while the compressive stiffness is E c for stresses σ c not exceeding the plastic concrete Similarly, the reinforcement stiffness is E s for absolute values of the stress σ s smaller than the yield stress f Y and the reduced value H s for stresses beyond this limit. This model corresponds to what was used in the framework presented in reference [9] for RC panels with in-plane forces. Note that the effect of transverse cracking is included in the model through a reduction of f cp based on the concrete effectiveness factor. For a bilinear-elastic material, a uniaxial stress state may be decomposed as where σ ℓ and σ h are the uniaxial stress components related to the linear and hardening response, respectively. With this decomposition, the corresponding complementary elastic energy can be expressed as The reinforcement stress state is assumed uniaxial in each direction, while the concrete stress state is assumed biaxial and in-plane at all points within the structure, that is, transverse shear stresses are neglected. This assumption greatly reduces the number of variables as well as the complexity of the concrete failure criterion, and since the transverse shear stress only contributes marginally to the complementary energy for thin structures (through transverse shear deformation), sufficient accuracy is maintained. Furthermore, Poisson's ratio for cracked concrete is taken as zero. With these assumptions, the uniaxial stressstrain curves for reinforcement and concrete can be applied directly to the reinforcement stresses {σ sx , σ sy } and the concrete principal stresses {σ cI , σ cII }, respectively.
By introducing the auxiliary variables {σ cm , σ cd , σ cr } (representing the mean axial stress, the difference between σ cm and the axial stress components, and the radius of Mohr's circle) using the constraints the concrete principal stresses can be introduced in the model as Thus, by decomposing the concrete principal stresses and the reinforcement stresses, applying the appropriate linear-elastic restraints to the linear stress components, as well as enforcing σ cI,h ≤ 0 since the concrete is only assumed to be ductile in compression, the formulation (11) can be applied directly to the concrete and reinforcement models.
F I G U R E 2 Complementary elastic energy of a material point with proportionality limit f

| Concrete effectiveness factor
As argued in reference [9], the plastic concrete compressive strength f cp appearing in (16) should be based on a reduction of the uniaxial compressive strength f c with the effectiveness factor ν, where ν accounts, inter alia, for the strength reduction due to cracking. In reference [9], this was taken into account by linking f cp = ν(ε I )f c to the stress components via the transverse strain ε I = ε xx + ε yy À ε II in a linearized expression, requiring a sequence of optimization problems to be solved in an iterative approach. The same approach is used in this paper.

| LAYER SUBMODEL
The method proposed in this paper uses a submodel approach where the section stress variation is modeled in a discrete set of points, and the resulting integral values (i.e., section forces and complementary elastic energy) are interpolated between these points. The concept is illustrated for a generic plane reinforced concrete shell with the thickness t located in the xy-plane with two layers of reinforcement oriented along the x-and y-axes in Figure 3.
Since the stress variation in a reinforced concrete section can be highly nonlinear, the discretization approach used in this paper is based on a piecewiselinear interpolation of the concrete stress σ c = [σ cx , σ cy , τ cxy ] T over a set of fixed layers: where I ℝ 3Â3 is the identity matrix, and z [z c,i , z c,i+1 ]. This concept is illustrated in Figure 4 for a section discretized into five layers. The reinforcement stress σ s = [σ sx , σ sy ] T is taken as constant in the reinforcement layers.

| Equilibrium
The in-plane section forces n = [n xx , n yy , n xy ] T and section moments m = [m xx , m yy , m xy ] T are related to the stress components by the relations: With a discretization of the section into n c layers, this can be stated in terms of the concrete and reinforcement stress components as m ¼ F I G U R E 3 Composition of a section in a generic plane reinforced concrete shell F I G U R E 4 Stress interpolation in a section using five layers where t c,i = z c,i+1 À z c,i is the layer thickness, n s is the number of reinforcement layers, a s,i = [A sx,i , A sy,i ] T contains the reinforcement areas per unit width, and J is the component-wise multiplication operator. Note that the reinforcement does not contribute to the shear forces and torsional moments.

| Complementary energy
The complementary elastic energy per unit area can be stated in terms of the concrete and reinforcement components, where the subscripts cI and cII refer to the concrete principal components, and the subscripts sx and sy refer to the x-and y-directional reinforcement components. With the discrete layer model, the complementary elastic energy is a piecewise-quadratic polynomial that can be integrated over each of the layers using simple quadrature: Here, z ip and w p are coordinates and weight for the integration point p in layer i according to Gaussian quadrature, and C 0 (z ip ) is given by (10) with the decomposed concrete stress variation represented in terms of the discrete components; see Equation (18).

| Kinematics
Despite the complementary energy minimization problem being formulated in terms of statics, it is possible to retrieve the kinematics from the dual solution (which is related to the primal solution through the Lagrangian function; see, e.g., reference [15] for the theoretical details). For a single section, this can be illustrated by considering the complementary energy minimization problem for a section with the in-plane section forces n and section moments m: The corresponding Lagrangian function is where {μ n , μ m } are the dual variables, giving the stationarity conditions: Given the definition r σ C 0 (σ) = ε(σ), it is seen that {μ n , μ m } can be interpreted as the center-plane strain Alternatively, given at least two reinforcement layers, the axial strain components {ε xx , ε yy } can be computed directly from {σ sx , σ sy } in the reinforcement layers and interpolated over the height, while the lowest principal strain ε II can be computed in each layer interface from σ cII . This approach is advantageous when a direct coupling between f cp = ν(ε I )f c and ε I = ε xx + ε yy À ε II is desired.

| EQUILIBRIUM FINITE ELEMENTS
An equilibrium-type finite element using n and m as the primary degrees-of-freedom (dofs) is used to model the section force variation. To ensure compatibility between the in-plane and out-of-plane section forces of noncoplanar elements, a linear and quadratic variation are assumed for n and m, respectively. The element contains 10 submodel points in which n and m are coupled to the stress variations in the discrete layer model. The element is illustrated in Figure 5. As described in the following, the element ensures continuity of in-plane forces, shear forces, and moments across element boundaries, although some researchers, notably in the context of limit analysis, have suggested to replace the torsional moment and shear force continuity conditions by equilibrium between Kirchhoff shear forces and corner forces. 18,19

| Equilibrium
Equilibrium between the generalized in-plane section forces and moments {n, m} and the domain load

be established based on Figures 6 and 7 as
where ∂ 1 and ∂ 2 are defined as ,

5: ð29Þ
Equilibrium across an element boundary is ensured as are the resulting boundary bending and torsional moment components lying in the element plane with the transformation given in terms of the boundary unit nor- Note that the element cannot represent moments about the element surface normal, meaning that torsional moments can only be transferred between two elements when these are coplanar. By considering the Figures 8 and 9, the relations between the boundary force and moment components and the section forces and moments are found as the following: Equations (32) to (34) are evaluated in the nodes containing the dofs for each adjoining boundary, again with v being eliminated using (29).

| Complementary energy
By interpolating the complementary elastic energy per unit area, C(x, y), between the values in a set of discrete points, C = [C 1 , …, C n ] T , its integral over the element area A can be approximated using Gaussian quadrature as where N(ζ) is an interpolation matrix given in terms of the area coordinates ζ = {ζ 1 , ζ 2 , ζ 3 }, and w p and ζ p are the Gauss integration weights and coordinates.
In the case of an isotropic and linear-elastic material, the polynomial variation of the section forces will carry over to the complementary energy, and (35) can be represented exactly using polynomial interpolation. In the case of an anisotropic, bilinear-elastic material as the ones considered in this paper, polynomial interpolation is still seen as the most appropriate choice when the material is not stressed beyond the linear-elastic limit. When stressed beyond the linear-elastic limit, a piecewise-linear interpolation between the submodel points may be more appropriate due to the local nature of plastic material behavior. For simplicity, however, polynomial interpolation is employed for this paper.
The choice of the assumed polynomial order of C(x, y) within a finite element is not only based on the usual consideration of accuracy versus computational efficiency but also the fact that the components of the coefficient vector resulting from (35) must be strictly positive in order to guarantee the conic constraint in (11) to be satisfied as an equality at optimality. This requirement is fulfilled for a cubic polynomial, while neither a quadratic nor a fourth-order polynomial fulfill this condition. Based on this, C(x, y) is assumed to be a cubic polynomial requiring 10 submodel points. It is worth mentioning that this number of points has been demonstrated to be reasonable for imposing yield criteria in the context of limit analysis. 11

| Kinematics
The interpretation of the dual variables for the optimization problem (24), {μ n , μ m }, as the generalized section strains, {ε 0 , κ}, can be extended to the complementary energy minimization problem for a plane domain rather straightforwardly: If the planar variation of C(x, y) is interpolated from a set of discrete points, its integral corresponds simply to a weighted sum of the discrete values of C, that is, in which case the dual variables should be normalized with the corresponding weight w i . Similar to how the dual variables for the section equilibrium equations can be interpreted as the generalized section strains, the dual variables for the boundary equilibrium equations can be interpreted as the boundary displacements. It was demonstrated by Vestergaard et al. 9 that a linear in-plane boundary displacement field can be computed as F I G U R E 8 Boundary force equilibrium F I G U R E 9 Boundary moment equilibrium u 1 where L is the boundary length, {u 1 , u 2 } are the boundary end-point displacements along one of the in-plane axes, and {μ q1 , μ q2 } are the dual variables for the in-plane boundary equilibrium equations along the corresponding axes. The transverse components of the boundary displacement fields can be computed using a similar approach by viewing each boundary as a onedimensional member acted upon by a linear transverse load (shear force) and a quadratic bending moment (torsional moment). The situation is illustrated in Figure 10. At the bottom, the discrete boundary equilibrium parameters performing work through the transverse displacement field (i.e., the transverse shear forces q nzi and torsional moments m nti ) are shown along with the corresponding dual variables, μ qnzi and μ mnti . At the top of the figure, the continuous polynomial variations, q nz t ð Þ, m nt t ð Þ È É , interpolated from the discrete boundary equilibrium parameters are shown together with a displacement field, w(t), interpolated from the values of the displacements and rotations in the equilibrium nodes, {w i , θ nti }. Since the dual variables, {μ qnzi , μ mnti }, are the work conjugates of the discrete boundary equilibrium parameters in the optimization problem, q nzi , m nti È É , these two situations are equivalent in terms of performed work, and the dual variables can be transformed to the discrete nodal values of the corresponding continuous polynomial displacement and rotation fields.
Assuming a polynomial displacement field, and given that θ(t) = dw(t)/dt, the five dual variables can be related uniquely to a displacement field of fourth-order at most. However, since a fourth-order polynomial with zero value and gradient at the end-points will also have a zero gradient at the midpoint, no shape function exists for the midpoint rotation. Consequently, the dual variable for the midpoint torsional moment equilibrium equation is disregarded, and a cubic transverse displacement variation is assumed. Based on the equations for work equivalence, the transverse displacement field is given in terms of the dual variables as where {w 1 , w 2 } and {θ nt1 , θ nt2 } are the transverse displacements and rotations, respectively, of the boundary endpoints as indicated in Figure 10.

| Spurious modes
The issue of spurious kinematic modes caused by linear dependencies in the in-plane boundary equilibrium equations was addressed in reference [9]. As described, these linear dependencies appeared in the four different types of mesh configurations illustrated in Figure 11, and could be resolved by introducing an additional kinematic restraining equation in the dual system for each case. As seen from (32) and (34), the boundary moment equilibrium equations have a structure similar to the inplane boundary equilibrium equations, and consequently, linear dependencies arise for the same types of mesh configurations. Thus, a similar approach is used in this paper to remove these linear dependencies, where additional kinematic restraints are imposed on the endpoint rotations for each case via the dual system. The boundary end-point rotations around the boundary normal (in the element plane) are given by (40) Here, the product of [μ mnn1 , μ mnn2 , μ mnn3 ] T and m nn1 , m nn2 , m nn3 ½ T is assumed to be equal to the work performed by the θ nn (t) and m nn t ð Þ. The following constraints are applied for the mesh cases (a)-(d): a. The change in angle between the two boundaries, around one boundary's normal, is zero. b. The change in angle between the two parallel boundaries is zero. c, d. The changes in angle between the boundaries in each pair are equal.

| Section with constant moment and variable axial load
The effect of the layer discretization on the solution accuracy was investigated for stress-based finite beam elements in reference [20], where as few as five layers were found to give reasonably accurate results. This example investigates the effect of the number of layers on the solution accuracy for the shell element. For this, consider a plane section subjected to a constant bending moment m and a variable axial compressive force p; see Figure 12. With the main challenge for the layer model being to accurately represent the nonsmooth through-thickness stress variations occurring, for example, at the neutral axis, this example is well-suited for assessing the effect of the number of layers on the solution accuracy.
The center-plane (z = 0) strain ε 0 and the curvature κ of the section are computed using five and 50 layers, respectively, and are compared to the analytical solution in Figure 13. In the figure, the axial load p has been normalized with respect to p cr , and ε 0 and κ have been normalized with respect to the quantities ε 0,Y = ε 0 (p = 0) and κ Y = κ(p = 0).
The analytical curves demonstrate how the curvature approaches a constant value when the axial compression is increased (corresponding to the result for an uncracked section) and how the center-plane strain goes from tension to compression accompanied by a change in axial stiffness. For both ε 0 and κ, the solution computed using 50 layers follows the analytical solution closely, with the maximum relative error (with respect to ε 0,Y and κ Y ) being in the order of 0.01%. For the solution computed using five layers, the maximum relative error of 4% is found for p = 0, which is reduced to 1% for p = 0.4p cr and further to 0.1% for p = 0.6p cr .
The variation in error magnitudes is explained by the concrete stress profiles for p = 0, p = p cr /2, and p = p cr in Figure 14. The figure clearly shows how the stress profile is generally well-approximated by 50 layers, while the approximation is less accurate, particularly for p = 0, when using only five layers. It should be kept in mind that even in this case, the resulting relative error on the F I G U R E 1 2 Section subjected to a constant bending moment and a variable axial load section level is less than 5%, and the accuracy can presumably be increased significantly by simply dividing each of the outermost layers into two, that is, effectively using just seven layers, or by merely using layers of varying sizes.

| Beam with clamped supports
In this example, a beam with clamped supports subjected to a uniformly distributed transverse load is considered; see Figure 15. As illustrated, the vertical support reactions on the beam are modeled as concentrated reactions and uniformly distributed reactions over a given length a, respectively. The beam has the length L = 5 m, the width b = 100 mm, and the height t = 200 mm, and it is modeled as a shell with two layers of Ø12 mm longitudinal reinforcement per 50 mm positioned 30 mm from each of the faces. The material strength and stiffness parameters are taken as f Y = 500 MPa, E s = 210 GPa, and H s = E s /250 for the reinforcement (this corresponds to a reinforcement stress magnitude at the ultimate strain ε u = 0.05, which is 8% higher than f Y corresponding to Class B steel cf. Eurocode 2 3 ), and f c = 35 MPa, E c = 34 GPa, and H c = E c /2500 for the concrete. The beam is modeled using symmetry conditions and a regular mesh of 24 square four-element patches over the first 2.4 m from the line of symmetry and six rectangular fourelement patches over the last 0.1 m. The cross-section variation is modeled using five concrete layers. The load-displacement curves for the beam midpoint are generated for the case with a concentrated support reaction and two cases with distributed reactions over a = t/4 and a = t/2, respectively. These curves are shown in Figure 16 in which the load p has been normalized with the limit load: where m u = 174 kNm/m is the ultimate moment (with 5% reinforcement strain, occurring at the support), and m Y = 158 kNm/m is the yield moment (occurring at the beam midpoint). the ultimate strain ε u . It is seen that in the case of a concentrated support reaction, the reinforcement reaches the ultimate strain at p/p u ' 0.9, while this happens at p ≥ p u for both cases with distributed support reactions. This demonstrates the potential impact of the model idealization on the results; here, in terms of high localized reinforcement strain values at a concentrated support that, if taken at face value, suggests that the limit load cannot be attained due to reinforcement rupture. In this example, the effect is decreased by modeling the supports using distributed reactions, but it is important to note that several other idealizations such as full strain compatibility and negligible transverse shear also have an effect. In practical design scenarios, a simple approach for handling idealization effects might be based on local strain averaging. This approach is also illustrated in Figure 16 for the case with concentrated support reactions, where the load levels for which ε max s ¼ ε u are indicated assuming strain averaging over d = t/4 and d = t/2, respectively. Here, it is seen that ε max s < ε u at p = p u even for d = t/4.

| EXAMPLE: CLAMPED SQUARE SLAB
Consider a clamped square slab of reinforced concrete loaded by the uniformly distributed, transverse load p as illustrated in Figure 17. This is a classic example from the field of limit analysis which Fox 21 demonstrated to have the exact load-carrying capacity p u = 42.851m u /L 2 , where m u is the ultimate moment capacity of the slab, as well as a rather complicated yield line pattern.
Being based upon the theory of rigid plasticity, the solution only grants information on the behavior in the ultimate state and not the path leading to this state. Still, it is of interest to compare this solution to the results predicted by the method presented in this paper. In order to allow a direct comparison with the exact solution, a constant effectiveness factor of ν = 1 is assumed throughout this example.
In the current example, the slab has the side lengths L = 5 m and the uniform thickness t = 200 mm, and it is reinforced with two layers of orthogonal mesh reinforcement of Ø20 mm bars per 150 mm positioned 30 mm from each face. As in the previous example, the materials are taken as being bilinear-elastic with E s = 210 GPa, H s = E s /250, and f Y = 500 MPa for the reinforcement, and E c = 34 GPa, H c = E c /2500, and f c = 35 MPa for the concrete. A quarter of the slab is modeled using symmetry conditions, a regular mesh with 8 Â 8 square patches of four elements, and five concrete layers over the slab thickness. Figures 18 and 19 show the principal moment magnitudes and orientations computed using the proposed method for three different load levels, namely p/ p u = 0.30, p/p u = 0.60, and p/p u = 0.95. In these figures, the transition between the elastic and plastic regimes is quite clear, with redistribution of the moments taking place as certain areas reach the yield capacity. Figure 20 shows a set of load-displacement curves for the slab midpoint. The four different curves are generated based on cases with and without in-plane axial supports, and with and without material hardening. For each curve, the load levels corresponding to reinforcement rupture (i.e., a maximum reinforcement strain corresponding to the ultimate strain) are indicated by a change in the line type. In the two cases with in-plane axial supports, the response is significantly more stiff than the remaining cases, and the analytical ultimate load is surpassed with only a small decrease of the initial stiffness. While these results may not be surprising, it does show the ease with which cases not covered by existing analytical solutions can be analyzed with the proposed shell element.
In the cases without in-plane axial supports, the stiffness reduces significantly at higher load levels, and the case with hardening reaches the ultimate load with a midpoint displacement of 100 mm (in the case without hardening, the ultimate load is a true upper bound and cannot be attained). It is also seen that the ultimate strain is reached for the reinforcement before 60% of the limit load. As in the previous example, this is a consequence of the idealized model, and simple strain averaging over d = t/4 and d = t/2 from the supports increases this to 85% and 92% of the ultimate load in the case with hardening, and more than 95% in the case without hardening. These levels are only indicated in Figure 20 for the former case since the displacements in the latter case exceed the axis limits.
In order to evaluate the method's ability to represent the collapse mode at high load levels, the clamped square slab without in-plane axial supports is now modeled using a regular mesh with 20 Â 20 square patches of four elements (again using symmetry conditions). As a visualization of the collapse mode, the approximated principal F I G U R E 1 7 Square clamped slab curvature fields for p/p u = 0.97 in the case without hardening are shown in Figure 21 where the color axis has been cut off at the ultimate curvature ± κ u (corresponding to ε s = ε u ) to allow a more explicit visual representation of the plastic zones. In the figure, the zones from the exact solution 21 are also indicated using enumerated points and connecting dotted lines. These zones consist of the ruled region AEC, the anticlastic region AED, and the undeformed region DBE. Although the zone locations do not correspond fully with the numerical solution, they are seen to be qualitatively similar, with a ruled region dominated by one principal component, an anticlastic region containing both principal components, and an undeformed region. Note that also in the exact solution, the curvatures in the region ACD are small relative to CDE.

| EXAMPLE: FOUR-STORY STAIRWELL WITH OPENINGS
In this example, the proposed method's applicability for modeling large structures is demonstrated for a four-story In the analysis performed in reference [22], the limit load was found using a fine mesh and design values for the material strengths as p ' 90 kN/m, with the reinforcement being the limiting factor. In the following, this value adjusted with the reinforcement safety factor, that is, p lim = 1.2p = 108 kN/m, is taken as the reference limit load.
In the current example, the stairwell is modeled using a mesh of 2439 finite elements and five concrete layers over the wall thickness. This discretization brings the number of optimization problem variables up to 3.43 million when material hardening is included and 1.77 million when material hardening is not included. The problem is solved on a personal computer with an Intel Xeon Silver 4110 CPU 2.10 GHz (4 processors) and 32 GB RAM with a solution time of 190 s per optimization problem when material hardening is included and 105 s when material hardening is not included. Note that inclusion of the concrete effectiveness factor requires load stepping with 1-7 iterations in each step depending on the effect of the change in the plastic concrete strength in the given step, and that material hardening is essentially a requirement for relating the concrete effectiveness factor to the stresses as described for membrane elements in reference [9]. However, due to the stepwise and iterative nature of that approach, the majority of the hardening variables and stress bounds can be omitted for any given iteration.
The load-displacement curves are generated for four different cases, namely assuming shell and membrane (i.e., in-plane action only) behavior, respectively, and with and without material hardening. In the cases without material hardening, ν = 1 is assumed. These curves are shown in Figure 23. In the case with in-plane action only and no hardening, the response is bounded above by p lim , while the material hardening is seen to increase the post-yielding stiffness with similar amounts in both cases. Regardless of the hardening, the bending capacity affects the stiffness significantly, and the ultimate load (not considering strain limitations) is increased by more than 50%.
The principal moment fields in the deformed stairwell at p/p lim = 1.5 are shown in Figure 24 for the case with hardening where the moment levels are given in terms of the yielding moment m Y = 26.4 kNm/m. From the figure, the bending action is quite evident in terms of the deformed geometry and the magnitude of the moments, which are concentrated around the openings and the wall edges.
F I G U R E 2 2 Four-story stairwell with mesh F I G U R E 2 3 Load-displacement curves for the stairwell modeled using shell and membrane elements, respectively The maximum of the principal strain fields is taken over the wall thickness for p/p lim = 1.5 and shown in Figure 25, and the development of these extreme principal strains is shown for the section of the stairwell front with coordinates (x/b, z/h) = (0.1, 0.7) in Figure 26. In both figures, the strain levels are given in terms of ε u = 0.05 and ε cu = 0.0035 for the largest and smallest principal strain component, respectively. Here, it is seen that significant strains are present, especially near the openings and the wall edges. This illustrates the price accompanying the increased stiffness provided by the bending capacity, namely that it must be ensured that sufficient deformation capacity is present in the parts subjected to moment action.

| CONCLUSIONS
A stress-based shell finite element framework for efficient elasto-plastic analysis of cracked reinforced concrete wall structures has been presented. With a scope limited to limit state analysis of thin, appropriately reinforced structures ensuring a stable crack growth, nonlinear elastic constitutive models for reinforcement steel and cracked concrete was used in combination with the principle of minimum complementary energy to cast the general stress-field problem as a second-order cone program which was solved efficiently using commercial convex optimization algorithms. The proposed element couples the variation of the section forces and the complementary elastic energy to a layer submodel for the throughthickness stress variation. The model was verified by a sequence of examples with exact solutions of increasing complexity, and the computational efficiency of the model in the context of design scenarios involving large structures was demonstrated on a model of a four-story stairwell with more than three million variables which was analyzed within minutes on a personal computer.

ACKNOWLEDGMENTS
The work presented in this paper is partly funded by the Innovation Fund Denmark under File No. 9065-00170B, and the Ramboll Foundation under File No. 2018-103.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.
ORCID Daniel Vestergaard https://orcid.org/0000-0002-9645-6904 F I G U R E 2 4 First (left) and second (right) principal moment fields in stairwell at p/p lim = 1.5 F I G U R E 2 5 Maximum first principal strain field (left) and minimum second principal strain field (right) in stairwell at p/ p lim = 1.5 F I G U R E 2 6 Extreme principal strain development in the section of the stairwell front with coordinates (x/b, z/h) = (0.1, 0.7)