Design‐oriented nonlinear‐elastic buckling analysis of reinforced concrete wall structures using convex optimization

The structural response of slender reinforced concrete (RC) structures may be highly nonlinear due to cracking, reinforcement yielding, and geometrically nonlinear effects. While advanced models can simulate the detailed response of such structures, they are generally ill‐suited for limit state verification in practical design scenarios due to a high computational and modeling effort. A design‐oriented method for evaluating the geometrically linear response of cracked RC wall structures was recently presented and demonstrated to allow the analysis of large‐scale models with more than 2400 finite shell elements within minutes on a standard personal computer. This paper proposes a design‐oriented numerical method for efficient instability analysis of slender RC wall structures, also enabling the inclusion of thermal effects. Based on a two‐step linearized buckling analysis, the method first determines the geometrically linear structural response by solving the complementary energy minimization problem, including, if relevant, thermal strains and reduced material stiffness. This solution is used to derive the sectional stiffness upon which a linearized buckling problem is formulated and subsequently solved as a linear eigenvalue problem. The model is validated using examples with exact solutions, and its applicability to large‐scale models is demonstrated through an example with a four‐story stairwell modeled using more than 3200 finite elements.

constitutive models for concrete, for example, based on fracture mechanics or damage-plasticity, are capable of taking into account these effects (see, e.g., Ref. [1]). However, the computational robustness and efficiency of these methods are generally not sufficiently high for use in practical design scenarios, and the high level of modeling details required by them often exceeds the information available in such design scenarios. Furthermore, as the loading history prior to the limit state considered is not known, taking into account the uncracked stiffness of parts of the structure may be unconservative and possibly less accurate than assuming fully cracked concrete.
Due to the lack of efficient numerical tools for design verification, a detailed analysis is, in practice, usually avoided by simplifying the wall into wall strips carrying primarily axial compression forces and without consideration of in-plane shear stresses. These wall strips can be analyzed by simplified methods for columns provided by design standards such as the Eurocode 2, 2 where, for example, the Euler load is adjusted using nominal stiffness parameters to account for cracking and material nonlinearity. However, such a simplification of the wall is far from ideal for several reasons: Firstly, in cases with substantial in-plane shear forces and/or the presence of openings, it may be difficult or even impossible in a sensible way to divide the wall into strips. Secondly, this type of simplification of the wall often leads to stability loads that are very conservative due to the disregard of possible bracing effects from, for example, transverse walls. On the other hand, it is difficult to ensure a safe load level since the stress distribution is influenced by the stiffness of the actual structure, and local buckling may occur.
To determine the collapse load as well as the loaddisplacement response of cracked RC wall structures, the authors have recently developed a design-oriented method based on convex optimization 3,4 exhibiting high computational robustness and efficiency. Using nonlinear elasticity to imitate elasto-plastic behavior for quasi-static loading, this method considers the principle of minimum complementary energy as a convex optimization problem (more specifically, a conic quadratic optimization problem) to solve the nonlinearelastic field problem with high efficiency and robustness. This method has been demonstrated to enable the analysis of FE models with more than 13,000 membrane elements or 2400 shell elements within minutes. In this paper, this approach is extended to include nonlinear buckling analysis, including the effects of high temperatures during fire.
The possibility of including the effects of high temperatures is highly relevant in the design of slender RC structures, as these are particularly prone to instability in fire scenarios due to the effects of heating. Modeling of fire scenarios consists of three levels of complexity: (1) the constitutive response to heating, (2) the heat development within the structural members, and (3) the fire development outside the structural members (see, e.g., Ref. [5] for an overview of issues in computational modeling of RC structures in fire). The fire development in buildings is in itself complex and depends on a range of unknowns, and it is difficult to translate to a realistic experimental setting. Consequently, although the temperature development in RC sections is strongly affected by the heating rate, duration and spatial variability, most experimental data is based on structural members subjected to standard fire conditions. 6 Based on such standard fire conditions, the development of the temperature profile within RC sections can be determined either through a computational thermal analysis (as done in the works by Refs. [7] and [8]) or using simplified closed-form models such as the model by Ref. [9] (as employed by Ref. [10]). Due to the layer of computational complexity added by the former option, closed-form models are much preferred in a design-oriented setting. The Eurocode draft, prEN 1992-1-2:2021, 11 provides closed-form expressions for the temperature profile and the thermal strain development in RC, as well as tabulated data for the strength and stiffness degradation as a function of the temperature. In addition to thermal strain and strength and stiffness degradation, the effect of fire-induced spalling might impact the load-carrying capacity of RC members, particularly in the case of highstrength concrete. Due to fire-induced spalling being complex and influenced by a range of interrelated parameters, prEN 1992-1-2:2021 suggests addressing spalling either by means of preventive measures or by taking into account the strength loss due to spalling by omitting the spalled layer of concrete from the geometry based on experimental evidence obtained from representative tests.
The method for instability analysis of slender RC wall structures proposed in this paper is based on the concept of a two-step linearized buckling analysis. In order to keep the modeling and computational complexity at a level suitable for application in design processes, the model assumes fully cracked concrete and nonlinearelastic behavior of the concrete and reinforcement. With these assumptions, the tangent sectional stiffness of the cracked structure is determined based on the solution to an initial geometrically linear analysis and subsequently used to estimate the critical buckling load. This procedure is repeated for a number of load steps of increasing magnitude until the predicted critical load corresponds to the applied load, at which point the linearization error vanishes. The effect of high temperatures is included through the constitutive model provided in the prEN 1992-1-2:2021. Although a reduced effective cross-section is trivial to implement, the effect of fire-induced spalling is not accounted for in the presented method as tabulated data on the extent of the spalled layers is not available.
The paper is structured as follows: First, a brief theoretical introduction is given to the variational methods used to formulate the proposed model, followed by a description of the analysis procedure. After this, the model basis is presented in terms of the constitutive model, the inclusion of fire-induced effects, and the finite elements (FEs) used for problem discretization. The expression for the tangent section stiffness is derived and used in a subsequent derivation of the linearized buckling problem. The model is verified using two verification examples addressing the modeling of thermal behavior for two simple stress states and the buckling behavior of a cruciform column, respectively. Finally, the model is applied to a four-storey stabilizing stairwell with door openings demonstrating the model's applicability to large-scale problems, and some aspects of reliability in design verification are discussed.

| VARIATIONAL METHODS
Consider an elastic body with the domain Ω and the boundary Γ ¼ Γ t , Γ s f g, which is acted upon by the prescribed body forces f in Ω, the prescribed boundary tractions t on Γ t , and the prescribed boundary displacements u on Γ s . In the body, the stresses σ and tractions t must satisfy static admissibility, the strains ε and displacements u must satisfy kinematic admissibility, and the stresses and strains must satisfy the constitutive equations. While these conditions can be stated directly as a set of simple-looking governing equations, exact solutions to those equations are only available in simple cases. Consequently, most structural engineering problems of practical interest are solved using approximate methods, for example, FE methods.
A powerful basis for establishing such approximations is the principle of virtual work (PVW), of which two variations exist. The first (the most well-known variation) is the PVW by variation in shape, stating that a body for which the first variation in the work, produced by any imposed kinematically admissible virtual strain and displacement field δε, δu f g, vanishes, static admissibility of the stress field will be granted. The second one is the PVW by variation in stress stating that a body for which the first variation in the work, produced by any imposed statically admissible stress and traction field δσ, t f g, vanishes, the strain and displacement fields will satisfy kinematic admissibility. Although the difference between the two formulations is subtle, the fact that the kinematics enter (1) in terms of variations while they enter (2) as finite quantities means that (1) also holds in the case of geometric nonlinearity, while (2) does not.
As an alternative to considering (1) and (2) directly, they may be invoked implicitly by seeking the local extrema of the potential energy functional Π p ε, u ð Þ or the complementary energy functional Π c σ, t ð Þ, respectively. Utilizing that δu ¼ 0 must hold on Γ s to be kinematically admissible, the potential energy functional takes the form where is the potential elastic energy density. Similarly, since δt ¼ 0 must hold on Γ t , the complementary energy functional takes the form where is the complementary elastic energy density. Assuming a linear strain-displacement relation and an elastic constitutive model (i.e., a unique relation between σ and ε), the functionals (3) and (5) are convex, meaning they each have a unique minimizer. In this case, the functionals form the basis of the principle of minimum potential energy 12 and the principle of minimum complementary energy, 13 which can both be stated as convex minimization problems and solved efficiently by convex optimization algorithms. The form of the complementary energy minimization problem has been described in detail in Refs. [3 and 4].

| ANALYSIS PROCEDURE
Analogous to the modification of Euler's critical buckling load suggested by Ref. [14], in which the elastic modulus is substituted with the tangent modulus to account for material nonlinearity, the method proposed in this paper makes use of the tangent sectional stiffness for posing the linearized buckling problem as a linear eigenvalue problem. The concept, which is similar to what was proposed by Ref. [15], is illustrated in Figure 1 for the onedimensional case. As illustrated, the method estimates the critical buckling load P cr based on the response to the applied load P 0 . Note that the linearization is exact at the point of linearization (and a good estimate in the vicinity), that is, when the critical load corresponds to the applied load.
Since the tangent section stiffness depends on the structural response to the applied load, the method consists of a sequential procedure composed of two steps for each load step, namely (1) determination of the tangent section stiffness based on the response to the current load, and (2) estimation of the critical buckling load based on a linearized buckling analysis. Note that since this procedure considers decreasing material stiffness with increasing strain, the proposed method captures failure governed by stiffness, strength, and combinations of the two.
As mentioned in Section 1, the framework proposed by Ref. [4] has proven extremely efficient and robust for stressbased analysis of cracked RC wall structures with nonlinear material behavior due to the use of state-of-the-art commercial convex optimization solvers. Being based upon the principle of minimum complementary energy, this framework cannot be used directly as the basis for geometrically nonlinear analysis. However, the framework is suitable for computing the response to the current load used to determine the tangent section stiffness.
In order to account for the adverse effects of initial imperfections, a scaled version of the initial buckling mode prediction is imposed upon the original geometry. Furthermore, to determine the tangent section stiffness as reliably as possible, the equilibrium equations are updated iteratively in each load step to ensure equilibrium in the deformed configuration.

| Material model
The method proposed in this paper models the behavior of RC using separate material models for concrete and reinforcement, respectively. The material models are defined in terms of uniaxial stress-strain curves for both materials. Since the reinforcement stress state is assumed uniaxial along each reinforcement direction, the uniaxial model relates the axial strains ε xx , ε yy È É directly to the reinforcement stresses σ sx ,σ sy È É . For the concrete, which is assumed fully cracked and with a plane stress state, the uniaxial model relates the principal strains ε I ,ε II f g, to the corresponding principal concrete stresses σ cI , σ cII f g. As Poisson's ratio for cracked concrete is taken as zero, these principal stresses and strains are assumed to be coaligned. Note that for a plane stress state, the principal stress components σ I , σ II f gare related to the plane stress components σ xx , σ yy ,τ xy È É as where is the radius of Mohr's circle, and σ d ¼ σ xx À σ yy À Á =2 is the difference between σ m and the axial stress components. This also applies to a plane strain state.
In order to suit a conic quadratic formulation of the complementary energy minimization problem, the uniaxial stress-strain curves are defined as piecewise-linear, which allows approximation of any nonlinear stressstrain curve. Figure 2 illustrates this concept for a segment of a stress-strain curve with positive stresses. As indicated, the piecewise-linear approximation of the strain function can, in this case, be stated as where the stress components σ i satisfy and f iÀ1 and f i are the bounding stress levels for σ i . Here, the constant term ε 0 is relevant, for example, when thermal strains are present. Note that in the general case where the stress is not necessarily positive, the F I G U R E 2 Piecewise-linear approximation of a convex strain function.
inequalities in (9) should be reversed for the negative stress components. The corresponding approximation of the complementary elastic energy can be stated as the sum of two components related to the positive and negative stress components: Denoting with the combined subscript "AE" the positive/ negative components, each component of C 0 can be stated as where 1 ¼ 1,…,1 ½ T , σ AE contains the positive/negative stress components ordered by absolute magnitudes, and is the uniaxial compliance matrix corresponding to σ AE . The function defined in (11) is convex if and only if E i ≥ E iþ1 ≥ 0, that is, for stress-strain curves with a nonnegative and non-increasing stiffness.
With the slack variable α, the complementary elastic energy components can be expressed equivalently as by including the convex constraints where M 1 2 AE can be chosen as any matrix satisfying the Cholesky factorization:

| Section model
Since buckling takes place in slender structures, the section model considered in this paper is based on Kirchhoff shell theory for thin plates, that is, only the plane stress and strain components, σ ¼ σ xx ,σ yy , τ xy Â Ã T and ε ¼ ε xx , ε yy ,2ε xy Â Ã T , are considered. Considering a section of a generic plane RC shell (Figure 3), the inplane section forces n ¼ n xx , n yy , n xy Â Ã T and the section moments m ¼ m xx , m yy , m xy Â Ã T are related directly to the concrete stresses σ c ¼ σ cx , σ cy ,τ cxy Â Ã T and the reinforce- where J is the componentwise multiplication operator, and a s,i ¼ A sx,i ,A sy,i Â Ã T contains the cross-section areas per unit width of the i'th reinforcement layer with the coordinate z s,i . The in-plane strain components are assumed to vary linearly over the shell thickness, allowing them to be stated as where ε 0 ¼ ε 0 xx ,ε 0 yy ,2ε 0 xy h i T contains the center-plane strain components, and κ ¼ κ xx , κ yy ,2κ xy Â Ã T contains the curvature components.
Composition of a section in a generic plane RC shell.
Assuming moderate displacements, the centerplane strain is related to the displacement field u ¼ u x , u y ,u z Â Ã T by the Green strain measure, while the curvature κ is described with sufficient accuracy using the small-strain measure. With comma derivative notation, that is, u ,x ¼ ∂u=∂x, these relations can be stated as where u ε ¼ u x ,u y Â Ã T and u κ ¼ u z ½ T are the in-plane and out-of-plane displacement fields, and are the in-plane and out-of-plane differential operators. Note that (19) can be stated compactly as where e ε ¼ ε 0 , κ ½ T is introduced as the generalized strain, and the definition of the combined differential operator ∂ is implied.

| FIRE-INDUCED EFFECTS
The risk of buckling in RC structures is generally increased when exposed to fire due to the effects caused by heating. While the behavior of RC under fire exposure is acknowledged as a complicated matter, it is generally accepted in structural analysis to base simplified models on the material temperature ϑ evaluated based on standard fire conditions and to model the total strain of a material fiber subject to high temperatures using the decomposition: Here, the first two components are the mechanical and thermal strain, while the third component is the Load-Induced Thermal Strain (LITS) covering creep and socalled transient strain caused by heating of the concrete. Expressions for the LITS are proposed in the literature; however, the model suggested in the Eurocode draft for fire design of concrete structures (prEN 1992-1-2:2021) 11 accounts for LITS implicitly in terms of temperaturedependent material curves. With this approach, the fireinduced effects corresponding to a given temperature ϑ can essentially be taken into account using material stress-strain curves of the type (8). This paper considers fire-induced effects using the approach suggested in prEN 1992-1-2:2021.

| Material model
The material curves suggested by prEN 1992-1-2:2021 are specified for a set of discrete temperatures ranging from 20 C to 1200 C which may be interpolated linearly, while the thermal strain variation is specified as a function of the temperature. For illustration, the material curves suggested by prEN 1992-1-2:2021 for cold-worked reinforcement and calcareous-aggregate normal-strength concrete are shown for different temperatures in Figures 4 and 5, including thermal strain. Note that because the material curves are specified in the model for the principal stress components, the thermal strain only contributes to volumetric expansion. The piecewise-linear approximations used in this paper are made using 2 Â 4 line segments for the reinforcement curves (2 Â 3 with nonzero slope) and four line segments for the concrete curves (two with nonzero slope) such that the maximum error is limited to 5% of f Y and f c , respectively. The approximations are indicated in Figures 4 and 5 as well.

| Section model
Using the expression for the section heat profile ϑ z ð Þ provided by prEN 1992-1-2:2021, the material curve parameters can be evaluated at any point z Àt=2, t=2 ½ . In the layer model used in the stress-based analysis, the stress limitations are enforced in the layer interfaces, while the complementary elastic energy is based on a linear interpolation of the compliance, that is, the elastic modulus inverse.
In a section, the complementary elastic energy per unit area due to the thermal strain is Given the layer model's piecewise-linear stress variation, this can be restated using integration by parts as is the thermal strain antiderivative, and is the average value of η z ð Þ in the i'th layer.

| Relaxed stress-based element
Reference [4] proposed using a planar equilibrium-type shell element for geometrically linear stress-based analysis of cracked RC wall structures. Assuming linear and quadratic variations of n and m, respectively, the element enforces stress equilibrium throughout the element domain as where f ¼ ½p x , p y ,p z T is the body force vector, and e σ ¼ n, m ½ T is introduced as the generalized stress. Furthermore, the element ensures traction equilibrium across element boundaries for the in-plane forces components t n ¼ n nx , n ny Â Ã T , the out-of-plane force component where is the out-of-plane shear force vector, n ¼ n x ,n y Â Ã T is the in-plane boundary unit normal vector, and is the in-plane stress-to-traction transformation matrix. As pointed out by the authors, because the element does not contain in-plane moment components (i.e., acting about the element's normal axis), torsional traction can only be transferred between two elements of this type if the two elements are coplanar. Consequently, any geometrical changes (e.g., imposed initial imperfections or displacements) will disrupt element coplanarity (as well as cause problems with the elimination of spurious F I G U R E 4 prEN 1992-1-2:2021 curves ( ) and piecewiselinear approximation ( ) for cold-worked reinforcement subjected to high temperatures.
F I G U R E 5 prEN 1992-1-2:2021 curves ( ) and piecewiselinear approximation ( ) for calcareous-aggregate, normalstrength concrete subjected to high temperatures. kinematic modes; see Ref. [4] for details), making its use in the framework proposed in this paper problematic. In order to resolve this issue, a relaxed equilibrium-type element is introduced, which preserves equilibrium between bending traction components, m nn ¼ n T t m while replacing the rigorous equilibrium conditions for t n , t v , and m nt with equilibrium requirements for a set of threedimensional concentrated force quantities which satisfy traction equilibrium in an integral sense. In essence, the approach corresponds to a Kirchhoff relaxation of m nt followed by a distribution of the in-plane and resulting out-of-plane traction components to the concentrated nodal quantities.
The boundary continuity components are illustrated for the rigorous and relaxed equilibrium FEs in Figure 6, where the circumflex accent "ð Þ" is used to signify discrete nodal variables.
Integral-sense equilibrium across boundaries is ensured by requiring the work performed by the integral traction quantities b N nx , b N ny , b V nz n o to equal the work performed by the continuous traction variation n nx t ð Þ, n ny t ð Þ,v nz t ð Þ È É and torsional moment traction variation m nt t ð Þ: Here, L is the boundary length, θ nt t ð Þ ¼ du z t ð Þ=dt is the boundary rotation around the boundary normal, and b u x , b u y , b u z È É contain the discrete values of the continuous boundary displacement components u x t ð Þ, u y t ð Þ, u z t ð Þ È É in the locations of the concentrated force quantities. Assuming a quadratic displacement field (which is accurate for the in-plane displacement field in the case of a linear-elastic material), these requirements translate to the following relations between the rigorous and the relaxed element continuity components:

| Displacement-based element
In displacement-based FE analysis, the well-known Constant Strain Triangle (CST) element (see, e.g., Ref. [16]) and Specht's plate bending element 17 have proven simple and robust for modeling in-plane and out-of-plane (bending) behavior. Based on this, these elements are chosen for the displacement-based, linearized buckling analysis in the proposed method. The elements and their dofs are illustrated in Figure 7. Although the CST and Specht elements can only represent constant and nearly linear states of in-plane strain and curvature, respectively, they are computationally inexpensive, allowing them to be used in fine meshes. In order to utilize this property without requiring an equivalent, unnecessarily small mesh size for the more accurate stress-based elements, the mesh for the displacement-based analysis is generated by subdividing each FE used in the stress-based analysis into 9 FEs of equal size. Collecting the CST and Specht's element dofs in b u ε and b u κ , respectively, the three-dimensional element displacement field can be represented as where N ε ξ ð Þ and N κ ξ ð Þ are the CST and Specht displacement interpolation matrices. In the following, the dependence on the spatial coordinates ξ is made implicit to simplify notation.
The linear and nonlinear strain interpolation matrices B and G xx , G yy , G xy È É are defined in terms of N as F I G U R E 6 Discrete equilibrium components for a boundary of a rigorous (top) and relaxed (bottom) equilibrium finite element.
and by introducing the indicator vector I i ¼ δ 1i ,…, δ 6i ½ T (in terms of the Kronecker delta δ ij ) and the following short-hand function notation: the generalized strain-displacement relation (19) can be expressed as where e ε L and e ε NL are the linear and nonlinear generalized strain functions, respectively. For later reference, note that e ε NL b u ð Þ can be linearized as Any virtual displacement field δu is represented analogously to the physical displacement field, that is, as δu ¼ Nδb u, giving rise to a virtual generalized strain field of the form where δe ε L and δe ε NL are the linear and nonlinear generalized virtual strain functions, respectively.

| Transformation of rotational degrees-of-freedom
In the relaxed stress-based element, the displacements and rotations are extracted from the part of the dual solution corresponding to the equilibrium equations (see Ref. [4]). Since bending moment equilibrium is enforced for each boundary, the bending rotational dofs around the boundary axes define a quadratic variation along each boundary. In order to use the rotations from the stressbased analysis in the displacement-based buckling analysis, these rotations should be converted from boundary rotations to the Specht element nodal rotations. By considering the end-point rotational dofs b θ n 1 n 1 , b θ n 2 n 2 n o from each of the adjoining boundaries with the unit normal vectors n 1 , n 2 f g (Figure 8), the corresponding Specht where n 1 Â n 2 is used as a short-hand notation for the out-of-plane component of the cross-product of n 1 and n 2 in three-dimensional space.

| Tangent section stiffness
In the neighborhood of a given state of generalized strain and stress, e ε 0 ,e σ 0 f g, the generalized stress-strain relation is assumed to be well-approximated by the first-order Taylor expansion: where is the tangent constitutive matrix. As seen, D T is fully determined by the concrete stress-strain gradient ∂σ c =∂ε and the force-strain gradients of the reinforcement layers, where E sT ε ð Þ is the reinforcement tangent stiffness and ε i ¼ ε xx,i , ε yy,i ,2ε xy,i Â Ã T is the strain state in the reinforcement layer i. Since the model defines the concrete tangent stiffness E cT ε ð Þ in the principal directions (which are co-aligned for both stresses and strains), the concrete stress-strain gradient should be expressed using the chain rule. Viewing σ c as a function of the principle stress magnitudes σ cp ¼ σ cI , σ cII ½ T and orientation φ, the chain rule takes the form where cos 2φ has been chosen as the principal orientation angle variable for convenience. By introducing xy q , the double-angle sine and cosine can be expressed as and the derivatives of ε I ¼ ε m þ ε r and ε II ¼ ε m À ε r can be found directly as Based on (46) and (48), the remaining components of (47) can be established as and ∂σ c ∂ cos 2φ where it is utilized that d sin2φ ð Þ=d cos 2φ ð Þ¼ Àcot 2φ ¼ Àε d =ε xy . In effect, (47) yields a positive semidefinite matrix for ε r > 0. Since ε r ¼ 0 corresponds to an unstrained material, the concrete stiffness can be taken as uncracked and linear-elastic in this special case. Due to the complexity of the expression for ∂σ c =∂ε, the integrals in (43) are evaluated numerically.

| STABILITY ANALYSIS BY THE PRINCIPLE OF VIRTUAL WORK
The displacement field within an element is decomposed into a component related to the element nodal forces, u k b q ð Þ, and a load-independent perturbation, ϵu ⊥ where ϵ ( 1: Note that any constant displacement term (e.g., due to thermal strains) appears implicitly in u k b q ð Þ. In the neighborhood of a given state u 0 , b q 0 f g, the relation between the nodal load vector b q ¼ λb q 0 and the load-dependent nodal displacement vector b u k is assumed to be wellapproximated by the first-order Taylor expansion: Assuming no load-independent displacements in the current point, the strain field can be approximated as F I G U R E 8 Rotational dofs at a vertex of the stress-based element (black) and Specht's plate bending element (white). where is the tangent strain interpolation matrix. The virtual strains are assumed to be locally independent of λ, that is, With the linearized generalized stress-strain relation (42), the virtual internal work in a FE can be approximated as where B T b u k,0 Àe ε NL b u k,0 À Á ¼ e ε 0 by definition, and e σ T ¼ D T B T b u T . By expanding the expression and neglecting terms containing ϵ 2 , the virtual internal work can be approximated as To allow evaluation of the integrals containing D T , its components are assumed to vary linearly between the submodel points. By expanding (59) to the system level, and stating the external virtual work in terms of the current load vector b R 0 , the principle of virtual works leads to the following equation: It is seen that the first two terms require equilibrium of the current state, while the last term is a linear eigenvalue problem which must be satisfied for any nonzero perturbation ϵ. Note that this simplifies to the wellknown linear buckling problem when K T ¼ K and K gT ¼ K g .

| Thermal behavior for simple stress states
To validate the implementation of the stress-strain curve approximation, consider a square RC plate with the side length L ¼ 5 m and thickness t ¼ 200 mm, which is reinforced with two layers of orthogonal mesh reinforcement in the form of Ø12 mm per 150 mm with the distance z s ¼ AE70 mm to the center plane. The plate is modeled using four FEs as illustrated in Figure 9, and its in-plane response to a state of uniaxial tension/compression and pure shear, respectively, is investigated for a range of constant temperatures. The material parameters are based on the stress-strain curves defined by prEN 1992-1-2:2021 for cold-worked reinforcement and calcareous-aggregate C30 concrete. The load-strain curves are shown in Figures 10 and 11, demonstrating F I G U R E 9 A square plate subjected to uniaxial tension/ compression (left) and pure shear (right).
the numerical model to provide good approximations in both cases.

| Buckling of a cruciform column
Consider a column of length L ¼ 30 m with a cruciform cross-section loaded by the compressive load P; see Figure 12. The column is formed by two intersecting RC walls with the width b ¼ 3:4 m and thickness t ¼ 100 mm, and reinforced with two layers of orthogonal mesh reinforcement in the form of Ø8 mm per 100 mm positioned with the distance z s ¼ 20 mm to the wall center-plane. The column is fixed at x ¼ 0 and free at x ¼ L.
The material parameters are based on the stress-strain curves defined by prEN 1992-1-2:2021 for cold-worked reinforcement and calcareous-aggregate C30 concrete at 20 C. The analytical flexural buckling load is where EI is the effective bending stiffness of the cross-section, while the analytical torsional buckling load is where GA is the effective torsional stiffness of the cross-section. If the section is assumed uncracked with a constant stiffness corresponding to the initial moduli of elasticity E c ¼ 2G c ¼ 16:1 GPa and E s ¼ 200 GPa, the critical buckling load for these two modes are Note that these loads correspond to average stresses of 22.6 and 27.8 MPa, respectively, which are lower than (but comparable to) the concrete crushing strength. Each column wall is modeled using a fine structured mesh with 8 Â 60 tiles of four FEs. The modes corresponding to the five lowest buckling loads are shown in Figure 13, with the first two modes being flexural (the direction of mode 2 is normal to the viewing plane) and the remaining three modes being torsional. It is seen that the buckling loads for the flexural modes are 2.7% higher than the analytical result, while the buckling loads for the torsional modes deviate from the analytical value by 0:3%, 0:6%, and 2:0%.
If the mesh size is doubled, the buckling load deviation is increased by another 1:0% for the flexural modes, while the deviations for the torsional modes remain within 0:2% of the fine mesh deviations. Thus, despite predicting buckling loads with a decent accuracy for both mesh sizes, it is seen that the buckling load has converged further towards the exact solution for the torsional modes than the flexural modes. This illustrates the relatively poor performance of the CST element in representing in-plane due to its linear interpolation of displacements.
To determine the nonlinear critical buckling load of the cruciform column, a linearized buckling analysis is performed for an imperfect geometry generated by imposing upon the original geometry scaled versions of the initial modes. The analysis is carried out on four different geometries with imperfections based on the four cases: Mode 1 with max. eccentricity e ¼ 100 mm (M1-100), mode 3 with max. eccentricity e ¼ 100 mm (M3-100), mode 1 with max. eccentricity e ¼ 30 mm (M1-30), and mode 3 with max. eccentricity e ¼ 30 mm (M3-30). The analysis employs a structured mesh with 4 Â 30 tiles of four FEs. Figure 14 shows the development of the predicted critical buckling load P cr .
As seen from Figure 14, the initial critical buckling load prediction is P cr ¼ 16:2MN in all four cases, corresponding to a flexural buckling mode. Note that this is only slightly higher than the critical buckling mode for the perfect geometry. Around P ¼ 5MN, the critical buckling load for case M3-100 begins to deviate as the critical buckling mode transitions to a torsional mode, while the  ½ . Consequently, in contrast to the prediction based on uncracked concrete, the torsional mode is more critical than the flexural mode when the initial imperfection corresponds to case M3-100. Note that the calculation is terminated in the first iteration for which P cr < P, that is, the curve points beyond this limit have not necessarily converged.
In the cases failing in bending, that is, M1-100, M1-30, and M3-30, the ranges for the critical buckling load can be compared with analytical estimates based on a cross-sectional analysis and an assumed linear variation of curvature along x. This approach estimates the critical buckling load as P est cr ¼ 12:1MN for case M1-100 (which is slightly lower than the above result), and P est cr ¼ 14:7MN for cases M1-30 and M3-30 (which corresponds well with the above results).

| EXAMPLE: STAIRWELL WITH OPENINGS
Consider the four-storey stabilizing stairwell of RC with door openings illustrated in Figure 15  The structure is supported against translation in all three directions along the bottom edge of the faces normal to the y-axis, while the bottom edge of the faces normal to the x-axis is free. It is loaded in the y-direction by the scalable, uniform distributed load p y ¼ 112:5 kN=m and in the z-direction by a uniform distributed load consisting of a constant component p z0 ¼ 300 kN=m and a scalable component p z ¼ 100 kN=m.
The walls are t ¼ 100 mm thick and contain two layers of Ø10 mm orthogonal mesh reinforcement per 150 mm, with the rebar centers positioned 30 mm from the wall faces. The material parameters are based on the stressstrain curves defined by prEN 1992-1-2:2021 for coldworked reinforcement and calcareous-aggregate C30 concrete at 20 C.
The structure is modeled using an unstructured mesh with 3268 stress-based FEs (i.e., 29,412 displacement-based FEs) and five concrete layers over the wall thickness, bringing the number of optimization problem variables and eigenvalue problem variables to 7,604,978 and 90,318, respectively. The geometric imperfection is based on the critical mode for the perfect geometry with an amplitude of e ¼ 50 mm. The entire analysis required 19 iterations distributed between 9 load steps, leading to a total computation time of 115 min.
The predicted critical load factor λ cr is shown as a function of the applied load factor λ in Figure 16, and the initial and final critical buckling modes are shown in Figure 17.
It is seen that the curve for the predicted critical load factor crosses the line λ ¼ λ cr for a critical load factor slightly lower than λ cr ¼ 0:84, which is approx. 12% lower than the initial prediction. Furthermore, the buckling mode is seen to localize in the final load step due to the concentration of the compressive forces near the bottom corners of the stairwell front.
The stairwell is now subjected to a fire scenario where the external faces are exposed to 30 min of standard fire cf. prEN 1992-1-2:2021, resulting in a through-thickness variation of the concrete compressive strength f c , stiffness E c , and thermal strain ε th . The thermal effects are illustrated in Figure 18 in terms of the deformations of the unloaded stairwell.
The load parameters are taken as p y ¼ 67:5 kN=m, p z0 ¼ 90 kN=m, and p z ¼ 60 kN=m, while the mesh and imperfection amplitude remain unchanged. The structure is analyzed using the procedure described above, leading to 21 iterations in total distributed between 10 load steps and a combined computation time of 160 min in this scenario. Note that the computation time per iteration in this scenario is 24% higher than the previous one due to the increased complexity of the reinforcement stressstrain curve approximations.
F I G U R E 1 9 Development of the predicted buckling load for the stairwell with openings after 30 min of fire exposure.
F I G U R E 2 0 Buckling mode and minor in-plane principal force field for the stairwell after 30 min of fire exposure.
The predicted critical load factor λ cr is shown as a function of the applied load factor λ in Figure 19, while the initial and final critical buckling modes are shown in Figure 20.
In this scenario, the predicted critical load factor increases from the initial value of 0:74 to λ cr ¼ 0:96 when the line λ ¼ λ cr is crossed. This behavior contrasts the previous examples where the predicted critical buckling load generally decreased with the applied load and is presumably a consequence of the thermal deformations which do not align with the buckling mode or the imposed imperfection. In this case, the difference between the critical buckling modes in the initial and final steps is relatively small, as the concentration of forces in the final step is less pronounced than in the case without thermal loading. It is seen that the deformation in both modes is concentrated around the first-and secondfloor door openings, which is presumably a result of the effect of the thermal deformations shown in Figure 18.

| RELIABILITY ASPECTS IN DESIGN VERIFICATION
The reliability of design verifications of RC structures is most frequently ensured by assuming design values for the modeled loads and material properties, for example, based on the partial safety factor format. When using nonlinear methods, extra attention should be paid to choosing these design values, as the solution might be sensitive to the model input. This includes the combination and positioning of loads, structural imperfections, and the influence of the probabilistic distributions of material properties. Following the Eurocode 2 suggestions, whether a nonlinear model used for design verification is sensitive to the material properties may be investigated by checking whether the predicted loadcarrying capacity based on design values for either the reinforcement or concrete properties and mean values for the other is lower than the load-carrying capacity based on design values for both materials. As the presented model does not make use of the tensile strength or fracture energy of concrete, it is expected to be less sensitive to the probabilistic distributions of material properties than traditional nonlinear models. However, as highlighted in the verification example concerning buckling of a cruciform column, the buckling load predicted by the proposed method may depend considerably on the imposed imperfection. Consequently, actual design verifications carried out by the proposed method should make sure to include the most critical initial imperfection modes.

| CONCLUSIONS
A finite-element framework for design-oriented nonlinearelastic buckling analysis of RC wall structures has been presented. Using an iterative two-step procedure, where the geometrically linear problem is posed and solved as a convex optimization problem in the first step, and the linearized buckling problem is posed as a linear eigenvalue problem based on the tangent section stiffness in the second step, the method enables robust and efficient stability analysis of RC wall structures. The effect of thermal strains and stiffness degradation resulting from thermal loading was included in the model utilizing piecewise-linear approximations of the stress-strain curves suggested by the Eurocode, making the model suitable for limit state analysis of RC buildings in fire scenarios. The model was verified using two simple examples with analytical solutions, and its applicability in practical design scenarios was demonstrated for a FE model of a four-storey stabilizing stairwell with openings using more than 3200 FEs, which was analyzed on a personal computer in 2-2.5 h.