Tuan-Anh Nguyen

Quang-Huy Nguyen
email: quang-huy.nguyen@insa-rennes.fr

Hugues Somja

An enhanced finite element model for reinforced concrete members under torsion with consistent material parameters

Keywords: Reinforced concrete, torsion, warping, torque-twist, large rotation kut nagiba tlačnih štapova; posmik

This paper deals with the development of a non-linear finite element model for reinforced concrete members under torsion. Using multi-fiber approach and displacement-based formulation, an enhanced multi-fiber 3D beam is proposed for predicting the behavior of reinforced concrete elements under torsion. The sectional analysis under elastic torsion is considered following Saint-Venant torsional theory for beam in order to take into account the effect of warping phenomenon. Elastic behavior under large displacements is also investigated in this model using a second-order approximation of the Green-Lagrange strains. In the inelastic domain, after cracking, following the space truss theory, the whole element is assumed to act as a tube, meaning the applied torsional moment is resisted only by the shear flow in the wall of the tube. Then, the effective wall thickness of the member after cracking is determined by an empirical formulation developed by the authors. Moreover, the section is discretized into different regions following its material response.

In each region, depending on its characteristics, an appropriate constitutive material model is applied. For the concrete, the proposed behavior models are based on the Modified of Compression Field Theory and its extension. In order to correctly predict the torsional response, the authors proposed some modifications in the tensile behavior of concrete, based on experimental tests in torsion. In the elastic domain, the model is validated by comparing to the analytical solution and some results from other researches. In the inelastic field, a good agreement is obtained between a series of experimental tests and numerical results.

Introduction

Although considered as a major factor in the design of many kinds of civil engineering structures (such as bridges, electric pylons...), torsion in reinforced concrete (RC) elements did not draw much attention by design engineers and researchers before the 1960s [6]. Prior to this period, the basic approach for calculating the torsional strength of RC members was the space truss concept proposed firstly by Rausch in 1929, in which the torsional effect was included within a one-dimensional frame element based on membrane models and the inclination of diagonal concrete cracks was assumed to be 45˚, so the model was called fixed-angle truss model [START_REF] Raush | Design of reinforced concrete in torsion[END_REF]. Then, with the introduction of the first building code design rules for torsion in the 1971 ACI Building Code [48], the space truss theory was developed widely by many researchers: the variable-angle truss model which considered the combination of torsion and bending and/or shear, developed by Lampert & Thurlimann [START_REF] Lampert | Ultimate Strength and Design of Reinforced Concrete Beams in Torsion and Bending[END_REF], Elfgren [START_REF] Elfgren | Reinforced concrete beams loaded in combined torsion, bending and shear : a study of the ultimate loadcarrying capacity[END_REF], Mitchell & Collins [START_REF] Mitchell | Diagonal compression field theory -A rational model for structural concrete in pure torsion[END_REF], [START_REF] Collins | Shear and torsion design of prestressed and non-prestressed concrete beam[END_REF] and Vecchio & Collins [START_REF] Vecchio | The Modified Compression-Field Theory for Reinforced Concrete Elements Subjected to Shear[END_REF]; or the softened truss model developed by Hsu & Mo [START_REF] Hsu | Softening of concrete in torsional members-theory and tests (Part I)[END_REF] which focused on the torsional behavior.

In the space truss models, the peculiar and featured deformation of cross-section under the torsional effect, called the warping phenomenon, is taken into account by the classical Saint-Venant torsional constant, whose expression is explicit and can be found in the literature [START_REF] Timoshenko | Theory of Elasticity[END_REF]. Even though the warping phenomenon was included successfully in the space truss model by this torsional constant, its influence in the shear strain and stress state distribution can not be represented completely. For this, a sectional analysis model with equilibrium kinematic conditions and appropriate constitutive material laws is considered suitable to obtain the distribution of stresses and strains at sectional level. In this approach, the warping phenomenon is included by a warping function, which is not explicit and must be mathematically determined by an infinite series. For the elements with thin-walled or open cross-sections, as the warping is strongly restrained at supports, some extended parameters are required to consider this effect called non-uniform torsion [13], [14].

Whereas, for solid cross-sections, this warping restrained effect can be neglected and the phenomenon is called free warping.

For the sectional analysis approach, since more than thirty years, the fiber (or multi-fiber) models, based on the subdivision of cross section into a system of point elements (or fibers), have been developed and applied successfully in the analysis of normal forces in RC elements [17], [START_REF] Spacone | Fiber beam-column model for non-linear analysis of R/C frames: Part 1[END_REF]. Recently, the tangential forces including torsional effect are investigated in the fiber models, using different approaches for the constitutive material model. Using damage mechanics for the non-linear constitutive model, Mazars et al. proposed a displacement-based multi-fiber beam element in which the warping kinematic for torsion is investigated initially in elasticity and then extended in a locally non-linear way [START_REF] Mazars | Using multifiber beams to account for shear and torsion: Applications to concrete structural elements[END_REF]. The strain field due to pure torsion is obtained for each fiber by using a warping-conduction analogy method to solve the problem of warping for sections of various materials. Le Corvec presented a mixed-formulation of 3D multi-fiber beam that is able to capture the local effects due to constrained warping of the section, so it can represent accurately 2 A c c e p t e d M a n u s c r i p t the torsional response of beam under warping constraints conditions [START_REF] Corvec | Nonlinear 3d frame element with multi-axial coupling under consideration of local effects[END_REF]. For this, the additional warping degrees of freedom are introduced at each section and then interpolated by Lagrange polynomials to ensure equilibrium on the entire element. In this model, the effects of boundaries on the warping distribution under torsion can be represented for arbitrary cross-sections, but the constitutive material model is limited only to steel. Based on this idea of warping interpolation using additional degrees of freedom, Addessi & Di Re extended Le Corvec's model to plain concrete member under torsion using an isotropic damage model for brittle-like materials [START_REF] Addessi | A 3D mixed frame element with multi-axial coupling for thin-walled structures with damage[END_REF], then Di Re used Hermite polynomials instead of Lagrange polynomials for the interpolation functions [START_REF] Di Re | 3D beam-column finite element under non-uniform shear stress distribution due to shear and torsion[END_REF]. Another model for torsional responses of RC element using damage approach for constitutive law was developed by Capdevielle, in which the warping was included by a governing equation following Saint-Venant's theory and solved by finite element formulations [22].

Besides damage mechanics, the smeared-crack approach was also employed by many researchers to take into account the torsional effect. Bairan & Mari proposed a general sectional model in which the torsional moment is included by solving the equilibrium between normal and shear stresses at the element and crosssection level [START_REF] Bairán | Multiaxial-coupled analysis of RC cross-sections subjected to combined forces[END_REF]. Warping effect is defined internally and included in an enhanced displacement field in addition to the classical displacement field following the plane-section theory. A non linear hypoelastic constitutive model was used and post-cracked concrete was modeled by a rotating-smeared crack approach.

The same approach of smeared-crack was also proposed by Navarro Gregori et al., with a discretization of the section into different regions following its stress state and an appropriate constitutive model based on the Modified of Compression Field Theory (MCFT) [START_REF] Navarro-Gregori | A 3D numerical model for reinforced and prestressed concrete elements subjected to combined axial, bending, shear and torsion loading[END_REF]. However this model neglect the effect of warping for torsion. Another frame model using forced-based formulation and 3D orthotropic smeared crack material model was also introduced by Kagermanov & Ceresa [START_REF] Kagermanov | 3D Fiber-Based Frame Element with Multiaxial Stress Interaction for RC Structures[END_REF].

In this present work, an enhanced multi-fiber 3D beam model is proposed for RC members under pure torsion. This model takes into account all the aspects of torsional effects, including the warping phenomenon, the behavior before and after cracking, the behavior under large rotation conditions, as well as the contribution of stirrups. It is worth to note that this proposed model is the improved version of the multi-fiber frame model proposed by Navarro Gregori [START_REF] Navarro-Gregori | A 3D numerical model for reinforced and prestressed concrete elements subjected to combined axial, bending, shear and torsion loading[END_REF], taking advantage of its sectional discretization which is very suitable for the specific behavior of RC elements under pure torsion. For the constitutive model of concrete, the MCFT is employed, due to its suitability for multi-fiber section discretization, as well as its simplicity and its wide use in engineering applications. Knowing that Navarro Gregori had not defined a rule to determine the thickness of the regions in his model, a parametric formulation for calculating this effective thickness was developed by the authors, based on a calibration study of numerous specimens of torsion tests. In addition, from this experimental calibration, a new constitutive model for concrete under torsional effect was proposed by the authors, based on the fact that numerical cracking torsional moments are reduced about half of the experimental values when using the original MCFT as concrete constitutive model [START_REF] Jeng | Softened membrane model for torsion in reinforced concrete members[END_REF]. The numerical examples are then compared to the analytical and experimental results for the 3 A c c e p t e d M a n u s c r i p t validation of the proposed model.

Theory of torsion

General

The problem of torsion in a straight member was first investigated in 1784 by Coulomb [START_REF] Coulomb | Recherches théoriques et expérimentales sur la force de torsion et sur l'élasticite des fils de metal[END_REF], during his study of the property of electric charges. When carrying out the torsional oscillation of a weight suspended on a wire, Coulomb found that torsional moment is proportional to the angle of twist. This discovery was then introduced in a theoretical equation for torsion by Navier in 1826 [START_REF] Navier | Résumé des lecons données à l'école des ponts et chaussées sur l'application mécanique à l'établissement des constructions et des machines[END_REF], in which the proportionality between torque and twist angle is represented by a constant called the torsional rigidity. This constant is equal to the multiplication of the the shear modulus G and the polar moment of inertia I p : T = GI p θ. However, this theory was limited to members with a circular section: torsional tests made by Duleau in 1820 noted a 20% difference between the moduli of rigidity for a circular and a square section [START_REF] Duleau | Essai théorique et expérimentale sur la résistance du fer forgé[END_REF]. This difference was then explained by a phenomenon called warping phenomenon, meaning that the cross sections, under torsional effect, are deformed and come out of their original plan. As a consequence, the sectional rigidity is changed explaining the difference of 20% in Duleau's observation. This warping phenomenon was first discovered by Saint-Venant in 1855. He solved the problem of torsion of a straight element with an arbitrary cross section by substituting the polar moment I p by a new constant C called Saint-Venant's torsional constant [START_REF] Saint-Venant | Mémoire sur la torsion des prismes[END_REF]. Due to the work of Saint-Venant, torsion was then considered a classical problem by Timoshenko in his famous theory of elasticity [START_REF] Timoshenko | Theory of Elasticity[END_REF].

Theory of torsion for Reinforced concrete members

The torsional response of a RC element can be divided into 2 different phases: before and after cracking of concrete. Before cracking, the response of the section is considered elastic and the behavior is very similar to that of a plain concrete member, which can be predicted by the Saint-Venant torsional theory. After cracking, the material is not continuous anymore, the theory of elasticity becomes useless and a new mechanism is required to interpret the torsional response in this phase. The most basic and common mechanism for a beam element under pure torsion is the space truss analogy.

The concept of truss model was first proposed by Ritter and Morsch in the early of 20th century, in order to simulate a RC beam subjected to shear. Under shear effect, diagonal cracks occur and separate the beam into a series of concrete struts. An analogy of truss model was therefore assumed, where the longitudinal reinforcement steels act like the top and bottom chords of the truss, while the transversal stirrups and concrete struts hold the role of web members (Figure 1a). In reality, the inclination of the concrete strut may vary, but in this initial model, it was assumed to be 45˚. This 45˚truss model, also called fixed-angle truss model, was then developed to carry out the modeling of the post-cracking behavior of RC members Since the early days of reinforced concrete the socalled classical truss analogy developed by Ritter and Mörsch (Ritter 1899, Mörsch 1908) [START_REF] Coulomb | Recherches théoriques et expérimentales sur la force de torsion et sur l'élasticite des fils de metal[END_REF][START_REF] Navier | Résumé des lecons données à l'école des ponts et chaussées sur l'application mécanique à l'établissement des constructions et des machines[END_REF] was proposed for shear design of reinforced concrete members (Fig. 1). The truss analogy is based on a truss model with parallel chords and web members connected by means of pin joints, where the concrete compressive struts are inclined at 45° with respect to the longitudinal axis of the beam while the shear reinforcement represents the tensile web members. According to Zilch and Zehetmaier [START_REF] Navier | Résumé des lecons données à l'école des ponts et chaussées sur l'application mécanique à l'établissement des constructions et des machines[END_REF], when the shear reinforcement (stirrups) is placed closely to each other the simple truss becomes a statically indeterminate truss (Fig. 2b). Generally, the truss model may be considered as a statically determinate simple truss composed of resultant forces from parallel tension and compression stress fields with pinned joints (Fig. 2c). subjected to bending-shear-torsion, with an assumption that the concrete struts take only compressive stress after cracking [START_REF] Elfgren | Reinforced concrete beams loaded in combined torsion, bending and shear : a study of the ultimate loadcarrying capacity[END_REF]. The theory is therefore named the Compressive Stress Field Theory.

Instead of using the theory of plasticity in the truss model, Mitchell & Collins focused on the strain compatibility of the truss model, and then derived a compatibility equation to determine the angle of the compression stress field. This theory permits the strain to be calculated by Mohr's circle, and was called the Compression Field Theory [START_REF] Mitchell | Diagonal compression field theory -A rational model for structural concrete in pure torsion[END_REF]. Then, based on this theory, Vecchio & Collins developed a constitutive model to evaluate the stress-strain relation of material points in RC members subjected to shear and axial stresses. This theory, called the Modified Compression Field Theory, was calibrated on various experimental tests of RC panels under plane-stress loading and was widely used to evaluate the constitutive models for concrete. Focusing on the torsional behavior, Hsu & Mo noted that in the variable-angle truss model, when using the stress-strain relation obtained from the concrete compression cylinder test, the predicted torsional strength becomes very unconservative [START_REF] Hsu | Softening of concrete in torsional members-theory and tests (Part I)[END_REF]. Hence, a softened-truss model was proposed by Hsu which is capable of correctly predicting the torsional strength and strains throughout loading history [START_REF] Hsu | Softened Truss Model Theory for Shear and Torsion[END_REF].

5

A c c e p t e d M a n u s c r i p t

Among the available models and theories above, the MCFT is very suitable to apply to the multi-fiber finite element model, because in this theory cracked concrete is simulated continuously, and cracking is handled as a distributed effect following its direction. Hence in an integration point of section there is enough information of the concrete strut's direction. In general, this approach is referred as to the smeared crack approach, regardless particularly suitable for the sectional analysis of members under torsion as well as combined loading [START_REF] Ceresa | A fibre flexure-shear model for seismic analysis of RC-framed structures[END_REF], [START_REF] Kagermanov | 3D Fiber-Based Frame Element with Multiaxial Stress Interaction for RC Structures[END_REF]. Otherwise, many constitutive models have been proposed and can be applied to RC members using multi-fiber approach, such as plasticity model or damage mechanics ( [START_REF] Mazars | Using multifiber beams to account for shear and torsion: Applications to concrete structural elements[END_REF], [22], [START_REF] Di Re | A mixed 3D corotational beam with cross-section warping for the analysis of damaging structures under large displacements[END_REF])...

Warping phenomenon

Under the torsional effect, cross-sections warp and come out of their own plane. In Figure 2, the signal + means that the cross-section come out following the positive direction of x axis, while the signal -indicates that the warping follows the negative direction of x axis. The iso-curves in solid represent the points having the same warping magnitude in the positive direction, while the iso-curves in dashed line represent the points having the same warping magnitude in the negative direction. This peculiar deformation reduces the sectional rigidity, generates the additional normal stresses which decrease the tangential stresses and so strongly influence the twist deformation. Named warping phenomenon, it was firstly studied by Saint-Venant in 1855 with the help of mathematical tools like Fourier series. By inventing a semi-inverse method, Saint-Venant represented the warping phenomenon by the so-called Saint-Venant warping function, which is restricted to linear behavior with two assumptions:

1. Cross-section shape remains unchanged after twisting.

2. Warping of the cross-section is identical throughout the length of the member. This assumption recalls to the problem of uniform torsion.

+ + + + - - - - + + + + + - - -- - z y

6

A c c e p t e d M a n u s c r i p t

According to Saint-Venant, the warping function is described as a solution of the Neumann problem:

∂ 2 ψ ∂y 2 + ∂ 2 ψ ∂z 2 = 0 in Ω (Laplace equation) (1a) ∂ψ ∂y n y + ∂ψ ∂z n z = zn y -yn z on G (Neumann boundary conditions) (1b)
where Ω represent the domain of cross section and G its boundary, n y and n z are the vector normal to the boundary. Exact solution can be expressed for any arbitrary kind of cross-section, for example the rectangular cross-section [16]:

ψ(y, z) = 4 π 3 ∞ n=1 (-1) n-1 (2n -1) 3 b 2 sinh 2n-1 b πy cos 2n-1 b π h 2 sin 2n -1 b πz -h 2 sinh 2n-1 h πz cos 2n-1 h π b 2 sin 2n -1 h πy ( 2
)
where b is the section width and h is the section height, with b ≤ h. A simplified exponential expression was proposed by Xu et al. [START_REF] Xu | Saint-Venant torsion of orthotropic bars with inhomogeneous rectangular cross section[END_REF]:

ψ(y, z) = A n exp β 1n z - h 2 + B n exp -β 2n z + h 2 - λC n hα 2 ξ 2 n -yz (3)
where h is the section height; α is the shear modulus ratio; λ is the gradient factor; and β 1n , β 2n , C n and ξ n are defined as follows:

β 1n = λ 2 4h 2 + α 2 ξ 2 n - λ 2h ; β 2n = λ 2 4h 2 + α 2 ξ 2 n + λ 2h ; ξ n = nπ b ; (n = 1, 3, 5, ...) C n = 8b π 2 (-1) n-1 2 n 2 (n = 1, 3, 5, ...)
with b is the section width and b ≤ h. If the material is isotropic, λ = 0, in this case:

A n = [1 -exp (-hβ 2n )] C n β 1n {1 -exp [-h (β 1n + β 2n )]} ; B n = [1 -exp (-hβ 1n )] C n β 2n {exp [-h (β 1n + β 2n )] -1} ; (n = 1, 3, 5, ...)
Otherwise, there are other approaches to interpolate the warping function to solid or composite crosssection, such as using the finite element technique [22] or using Lagrange and Hermite polynomials [START_REF] Corvec | Nonlinear 3d frame element with multi-axial coupling under consideration of local effects[END_REF], [START_REF] Di Re | A mixed 3D corotational beam with cross-section warping for the analysis of damaging structures under large displacements[END_REF]. For example, Capdevielle et al. presented a governing equation of warping function for a homogeneous material first (following Saint-Venant warping function), then expressed 2 conditions of Neumann (equation ( 1)) and finally applied finite element formulations to solve the warping function. Le Corvec, when developing a model for steel members with torsion, described the warping displacement as the product of two interpolation functions using Lagrange polynomials, and considered warping displacement as independent degrees of freedom of element. Di Re then developed this approach using Hermite polynomials for RC members. Solving the warping problem is indispensable to figure out the torsional response of a straight element with an arbitrary cross section in the elastic behavior.

7

A c c e p t e d M a n u s c r i p t

Enhanced 3D beam FE formulation

Based on the principle of discrete finite element method, the idea of multi-fiber finite element is to divide the structure into several longitudinal fibers and some control sections along the element. At the intersection of longitudinal fibers and control sections, a system of integration points is obtained (Figure 3). Each point, called fiber and considered as a material point, has its own coordinates, surface and an appropriate material law in order to determine the strain and stress from the element's displacements. Therefore, the proposed model is described as a frame element with a set of cross-sections along its longitudinal axis. The local frame system of coordinates (x, y, z) is defined in Figure 3. In every single point of the element axis subjected to pure torsion, the displacement field consists in only one rotation θ x (x) about the longitudinal axis. In a displacement-based formulation, this generalized displacement is determined from the nodal displacements by shape functions that are chosen linear in the present work [START_REF] Friedman | An Improved Two-Node Timoshenko Beam Finite Element[END_REF]:

θ x (x) = 1 -x L θ I x + x L θ J x = 1 -x L x L   θ I x θ J x   = N s (x)q e ( 4
)
where θ I x and θ J x are the twist angle at node I and J respectively; x is the coordinate of the section and L the element length; q e is the vector of nodal twist angle. The twist κ x is then equal to:

κ x = ∂θ x (x) ∂x = - 1 L 1 L   θ I x θ J x   = B s q e (5)

Sections kinematics and variational formulation under small displacement

The displacement of a material point has 3 components U, V, W representing the displacement field in

x, y and z directions of the adopted coordinate system, respectively, which are collected in a single vector d f (x, y, z). Under the hypothesis of small displacements, the following kinematic relation is established according to Saint-Venant torsional theory:

U (x, y, z) = ψ(y, z)α(x) V (x, y, z) = -zθ x (x) W (x, y, z) = yθ x (x) (6) 8
A c c e p t e d M a n u s c r i p t

Where the term ψ(y, z) is the Saint-Venant warping function which describes the profile of warping over the cross-section, while α(x) represents the distribution of warping along the element length, called the warping distribution. Remembering that without warping, there is no axial displacement for material points of element subjected to pure torsion. The axial displacement due to warping (called the warping displacement)

is assumed by Saint-Venant to be proportional to his warping function and to the distribution of warping over the element length. According to Saint-Venant torsional theory, the profile of warping is considered constant along the element, the warping displacement is thus independent of x, hence the warping distribution becomes constant and equal to the derivative of the twist angle:

α(x) = ∂θ x ∂x = κ x , and ∂α(x) ∂x = ∂κ x ∂x = 0.
It's worth to notice that, in the case of a thin-walled cross-section, as the role of warping becomes important, the warping distribution cannot be considered constant, but depends on the cross-section's position. The warping phenomenon is then called restrained warping. Thus, Vlasov proposed a new theory of torsion for thin-walled cross-section with α(x) = ∂θ x ∂x and ∂α(x) ∂x = 0 . Consequently, an additional degree of freedom needs to be added to consider the contribution of the derivative of parameter α(x). Otherwise, Benscoter proposed another torsional theory in which the warping distribution parameter α is independent of the torsional angle α(x) = ∂θ x ∂x . These two theories of torsion for thin-walled cross-sections, can be compared as the analogy of two classical bending theories of Navier-Bernoulli and Timoshenko. Vlasov's assumption of neglecting the shear warping deformation is compatible with neglecting the shear bending strain in the Navier-Bernoulli's theory, while Benscoter's assumption of incorporating shear warping deformation can be regarded similar to the Timoshenko's assumption of taking into account the shear bending strain. However, in the case of a solid cross-section, the warping effect is limited, but cannot be neglected, warping is considered free and its distribution is constant over the element length. The displacement field of free warping becomes:

U (x, y, z) = ψ(y, z)κ x V (x, y, z) = -zθ x (x) W (x, y, z) = yθ x (x) (7)
The strains of any material point of the cross-section are then evaluated with only 3 components considered in the sectional analysis as follows: one normal strain and two transverse strains collected in a single strain vector e f (y, z):

ε xx = ∂U ∂x = 0 γ xy = ∂U ∂y + ∂V ∂x = -z + ∂ψ ∂y κ x γ xz = ∂U ∂z + ∂W ∂x = y + ∂ψ ∂z κ x                    ⇒ e f (y, z) = a f (y, z)κ x (8) 9
A c c e p t e d M a n u s c r i p t

With the compatibility matrix:

a f (y, z) = 0 -z + ∂ψ ∂y y + ∂ψ ∂z T (9)
As can be seen, in the case of solid cross-sections, the use of Saint-Venant warping function to represent the warping has a disadvantage when it makes appear a null value of the axial strain. This drawback can be solved using different approaches for warping, such as Lagrange interpolation polynomials, which will be introduced in another paper by the authors.

Once the strain vector is obtained at each material point, an appropriate material law is applied to determine the material stresses, which are collected in a single stress vector

s f = σ xx τ xy τ xz T .
The behavior law will be introduced in the next section and is represented by a material constitutive matrix k f :

s f = k f e f (10)
Next, the element equilibrium is considered between internal and external potential energy. Let the element be subjected to a virtual displacement δd, then the principle of virtual work gives an equation between internal and external energy:

δW i = δW e ⇒ V δe T f s f dV = δq T e M x,e + L δθ x (x)T u dx (11)
where the virtual internal work is represented by the left-hand side, while the virtual external work is expressed by the right-hand, with the contribution of the external nodal torsional moment M x,e and the external uniform torque T u . Considering firstly the left-hand side, with the aid from equation (8), we obtain:

W i = V δe T f s f dV = V δκ x a T f s f dV (12)
As the term a f (y, z) T s f is a function of the fiber coordinates y, z, the virtual internal work can be decomposed into an integral over the element length and another over the cross-section area:

W i = L δκ x A a T f s f dA dx (13)
As a consequence, at sectional level, the torsional moment can be determined as an integral over the crosssection area of the stress field in the section:

M x,s (x) = A a T f s f dA = A -z + ∂ψ ∂y τ xy + y + ∂ψ ∂z τ xz dA (14)
With the aid of the constitutive relation in equation [START_REF] Collins | Shear and torsion design of prestressed and non-prestressed concrete beam[END_REF], the sectional stiffness matrix can be obtained in a compact form:

M x,s (x) = A a T f s f dA = A a T f k f e f dA = K s κ x (15) 10
A c c e p t e d M a n u s c r i p t so:

M x,s = K s κ x (16)
with:

K s = A a T f k f a f dA (17)
K s is defined as the sectional stiffness matrix.

Next, the right-hand side is utilized with the aid of equation ( 4):

W e = δq e T M x,e + L δθ T x T u dx = δq e T M x,e + L δq e T N T s T u dx = δq e T M x,e + L N T s T u dx (18)
Hence, the virtual work equilibrium from equation ( 11) has now become:

L δκ x M x,s dx = δq e T M x,e + L N T s T u dx (19)
With the relation of equation ( 5) and ( 16), now we have:

δq T e L B T s K s B s dx q e = δq T e M x,e + L N T s T u dx ⇒ L B T s K s B s dx q e = M x,e + L N T s T u dx (20)
The element stiffness matrix is thus defined as:

K e = L B T s K s B s dx (21)
And the element equilibrium becomes:

K e q e = M x,e + L N T s T u dx (22)

Section kinematics and Finite element formulation under large displacements

Under small displacements hypothesis, the element response is assumed in the linear geometrical conditions. In normal requirement of use, this assumption is enough for the kinematic and equilibrium conditions of structures. However, when subjected to extreme loads (such as earthquake or natural hazards...), structures may undergo large displacements, and the geometrical analysis of the structure becomes non-linear.

According to Battini & Pacoste, for members under torsional effects the geometrical nonlinearity is generated by a term included in the local strain definition, which is defined by a second-order approximation of the Green Lagrange strains [START_REF] Battini | Co-rotational beam elements with warping effects in stability problems[END_REF].

Under pure torsion, in a solid symmetric section, the second-order approximation of displacement vector of an arbitrary material point has 3 components:

U = ψκ x (23a) V = -zθ x - 1 2 zθ 2 x (23b) W = yθ x - 1 2 yθ 2 x (23c)

11

A c c e p t e d M a n u s c r i p t

The Green-Lagrange strain components are given by:

ε xx = ∂U ∂x + 1 2 ∂U ∂x 2 + ∂V ∂x 2 + ∂W ∂x 2 = 1 2 (y 2 + z 2 )κ 2 x ( 24a
)
γ xy = ∂U ∂y + ∂V ∂x + ∂U ∂x ∂U ∂y + ∂V ∂x ∂V ∂y + ∂W ∂x ∂W ∂y = -z + ψ y κ x (24b)
γ xz = ∂U ∂z + ∂V ∂x + ∂U ∂x ∂U ∂z + ∂V ∂x ∂V ∂z + ∂W ∂x ∂W ∂z = y + ψ z κ x (24c)
Unlike the material strain in linear geometry condition, the axial strain under large displacements is not zero and is called Wagner term, and it causes a non linearity in the response in pure torsion. Because of this term, the local strain can not be related to the generalized twist κ x in a compact form as in the above section. Instead, the nodal torsional moments and element stiffness matrix in a finite element framework will be derived from the strain energy function.

The strain energy is expressed as a function of the local strains:

Φ = L 0 Φ A dx = L 0 1 2 A Eε 2 xx dA + 1 2 A G γ 2 xy + γ 2 xz dA dx = 1 2 L 0 EI rr κ 4 x + GJκ 2 x dx (25)
With:

EI rr = A E(y, z) 1 4 (y 2 + z 2 )dA GJ = A G(y, z) ∂ψ ∂y -z 2 + ∂ψ ∂z + y 2 dA
Using equation ( 5), equation ( 25) becomes:

Φ = 1 2 L 0 EI rr κ 4 x + GJκ 2 x dx = 1 2 EI rr L 0 κ 4 x dx + GJ L 0 κ 2 x dx = 1 2 EI rr L 0 - θ I x L + θ J x L 4 dx + GJ L 0 - θ I x L + θ J x L 2 dx = 1 2
EI rr L 3 (θ J x -θ I x ) 4 + GJ L (θ J x -θ I x ) 2 (26)
The nodal torsional moment in each element is then evaluated by:

M x,e = ∂Φ ∂q e =      ∂Φ ∂θ I x ∂Φ ∂θ J x      =      -2 EI rr L 3 (θ J x -θ I x ) 3 - GJ L (θ J x -θ I x ) 2 EI rr L 3 (θ J x -θ I x ) 3 + GJ L (θ J x -θ I x )      =      M I x M J x      (27) 12
A c c e p t e d M a n u s c r i p t

And the element stiffness matrix:

K e = ∂M x,e ∂q e =      ∂M I x ∂θ I x ∂M I x ∂θ J x ∂M J x ∂θ I x ∂M J x ∂θ J x      =      6 EI rr L 3 (θ J x -θ I x ) 2 + GJ L -6 EI rr L 3 (θ J x -θ I x ) 2 - GJ L -6 EI rr L 3 (θ J x -θ I x ) 2 - GJ L 6 EI rr L 3 (θ J x -θ I x ) 2 + GJ L      (28)

Mechanical model for RC members

Multi-fiber model has proven its power and advantage in the modeling of frame members with solid or thin-walled cross-section under many conditions of load. Knowing that the element is discretized in the longitudinal direction, taking into account the contribution of transversal reinforcement in a multi-fiber model is not an obvious work. A force-based formulation finite element was proposed by Saritas & Filippou to include the transverse strain of shear effects in RC members, by satisfying the internal equilibrium between concrete and stirrups in the constitutive model [START_REF] Saritas | A beam finite element for shear critical RC beams[END_REF]. Navarro-Gregori et al. proposed a displacementbased finite element in which he discretized the section in some separated regions following the direction of transversal reinforcement and the contribution of stress state to the sectional response [START_REF] Navarro-Gregori | A 3D numerical model for reinforced and prestressed concrete elements subjected to combined axial, bending, shear and torsion loading[END_REF]. More recently, Khoder proposed a multi-fiber approach considering the lateral confinement of concrete due to stirrups, in which the transversal stresses of stirrups are interpolated into stress state using Lagrange polynomials [23].

The numerical results give good correlations, but limited in the two-dimensional field under shear-bending.

In this work, the authors borrow the idea proposed by Navarro-Gregori, knowing that this approach is accorded to the sectional analysis, so that it can be applied to the multi-fiber FE frame model. Moreover, the constitutive model is formulated general enough for being used for any shape of cross-section in further studies. The cross-section is discretized in 3 zones according to the stress state that they are subjected to (Figure 4). In each zone, an appropriate constitutive model is established to determine the regional stiffness matrix

K 1D s , K 2D s or K 3D s .
These stiffness matrices are computed by a secant-stiffness-based formulation, based 13

A c c e p t e d M a n u s c r i p t on the constitutive relations of MCFT [START_REF] Vecchio | Toward Compression-Field Analysis of Reinforced Concrete Solids[END_REF]. In this theory, the concrete relationship is formulated in the principal direction of crack, with the following primary assumptions:

1. Principal directions of strains and stresses are coincident.

2. The effect of concrete softening in compression is included.

3. The tension-stiffening effect is also included.

Then, for the entire section, the sectional stiffness matrix in equation ( 17) is obtained by summarizing all these regional stiffness matrices:

K s = K 1D s + K 2D s + K 3D s (29)

1D-Zone

This region takes into account only the contribution of the longitudinal reinforcing steel (rebar), and the only stress accounted for is the axial component σ xx , which can be easily computed from the axial strain using an uniaxial behavior law of steel. In this paper, the authors used a simple elasto-plastic steel model presented in Figure 5. The stress vector is therefore expressed by:

s 1D f =      σ 1D xx 0 0      =      E s ε m xx 0 0      =      E s 0 0 0 0 0 0 0 0           ε 1D xx γ xy γ xz      = k 1D f e 1D f ( 30
)
Where E s is the secant modulus of steel. The 1D-regional stiffness matrix K 1D s can then be determined using the kinematic condition in equation ( 17).

sy f s E sh E σ = s s E ε s ± εy ± ε s σ s 0 s E Figure 5: Steel's behavior

2D-Zone

This region corresponds to the portion in which the transverse steel crosses in one direction and may also have the contribution of longitudinal steel bars. As mentioned above in the theory of torsion, after

14

A c c e p t e d M a n u s c r i p t cracking of concrete, the whole member is assumed to act like a tube, and the wall thickness of this tube corresponds to the width of the 2D-zone. The 2D-zone therefore has a strong impact on the torsional behavior of RC members, and the determination of the effective wall thickness t e becomes an important and essential problem that will be addressed later.

Corresponding to the sectional analysis of a frame element, the stress state of this zone has only two non-zero components instead of three: one axial and one transversal which correspond to the direction of stirrups. This results in a stress vector

s 2D f -h = σ xx τ xy 0 T (zone in which stirrups are disposed in horizontal direction or horizontal stirrups -called 2D-horizontal zone) or s 2D f -v = σ xx 0 τ xz T (zone in
which stirrups are disposed in vertical direction or vertical stirrups -called 2D-vertical zone) (Figure 4). It's worth to note that, for a RC element, the contribution of transversal reinforcement is taken into account by considering a third stress component: σ yy or σ zz , depending on the direction of stirrups. This component must be considered firstly, and then will be imposed to zero in order to satisfy the internal transversal equilibrium of RC members. The stress vector to be considered in the element state therefore becomes

s 2D f -h = σ xx τ xy 0 σ yy 0 0 T or s 2D f -v = σ xx 0 τ xz 0 σ zz 0 T .
For the sake of generality, the full stress vector for the following will take the expression of s

2D f = σ xx τ xy τ xz σ yy σ zz τ yz T ,
where some components are zero depending on the stirrup's direction. Corresponding to this 2D-general stress state, the 2D-strain vector must also have 6 components e

2D f = ε xx γ xy γ xz ε yy ε zz γ yz T ,
but unlike in the stress state, ε yy (2D-horizontal zone) and ε zz (2D-vertical zone) are not imposed to zero, and must be determined in the process of satisfying the transversal equilibrium. However, in the sectional analysis thereafter, these strain components ε yy and/or ε zz could not be put in the strain vector e 2D f , which contains only 3 components: ε xx , γ xy and/or γ xz . In the other hand, ε yy and/or ε zz could be included in the strain vector e 2D f if the distortion phenomenon is taken into account, but this phenomenon is not considered by the Saint-Venant theory of torsion.

In this 2D-zone, the original MCFT is applied to determine the regional constitutive matrix K 2D s . Knowing that the original MCFT is designed for membrane members, it is practical to explore the constitutive formulation of this theory in an in-plane frame system of coordinate axes (l, t) following the longitudinal and transversal direction (Figure 6).

For this, a change of reference is applied to transform the strain vector from the local frame system (x, y, z) to the in-plane system (l, t):

e 2D ip = T ip e 2D f (31)
Where T ip is the transformation matrix, here as sequence, the subscript ip denotes the index of parameters expressed in the in-plane system:

T ip =               1 0 0 0 0 0 0 0 0 cos 2 α sin 2 α sin α cos α 0 cos α sin α 0 0 0 0 0 0 sin 2 α cos 2 α -sin α cos α 0 -sin α cos α 0 0 0 0 0 0 -sin 2α sin 2α cos 2α               (32)
α is the rotation angle between the two frames whose value is either 0 for the 2D-horizontal zone or π 2 for the 2D-vertical one. As a consequence, the in-plane strain vector becomes:

e 2D ip-h = ε xx ε yy γ xy ε zz γ xz γ yz T 2D-horizontal (33a) e 2D ip-v = ε xx ε zz γ xz ε yy -γ xy -γ yz T 2D-vertical (33b)
Only the first three components are considered in the in-plane coordinate system, the others are given a null value, so:

e 2D ip-h =      ε l ε t γ lt      =      ε xx ε yy γ xy      (2D-horizontal) or e 2D ip-v =      ε xx ε zz γ xz      (2D-vertical) (34)
Among these 3 in-plane strain components, ε l and γ lt can be obtained from the kinematic condition in equation ( 8) (small displacement condition) or equation [START_REF] Navarro-Gregori | A 3D numerical model for reinforced and prestressed concrete elements subjected to combined axial, bending, shear and torsion loading[END_REF] (large displacement condition), while ε t must be handled separately by satisfying the transversal equilibrium conditions. The determination follows an iterative process and will be described later. Corresponding to this strain vector, the in-plane stress vector has 3 components:

s 2D ip = σ l σ t τ lt T
. The stress and strain vectors are related by a material stiffness matrix D 2D ip :

s 2D ip = D 2D ip e 2D ip ( 35
)

16

A c c e p t e d M a n u s c r i p t

The contribution of concrete and reinforcement are added separately to the material stiffness matrix and the stress vector of the in-plane frame (Figure 7):

s 2D ip = s 2D ip,c + s 2D ip,s (36a)
D 2D ip = D 2D ip,c + D 2D ip,s (36b)
Where D 2D ip,c is the concrete stiffness matrix and D 2D ip,s is the stiffness matrix from reinforcement. While the steel stiffness matrix D 2D ip,s can be set up easily in the in-plane system, it's more convenient to express the concrete relation in principal directions (Figure 7). A change of reference from principal direction to in-plane axes systems is therefore required and equation (36b) becomes:

D 2D ip = T 2D c T D 2D prin,c T 2D c + D 2D ip,s (37)
Where T 2D c is the transformation matrix composed of the direction cosines, which define the direction of [START_REF] Vecchio | The Modified Compression-Field Theory for Reinforced Concrete Elements Subjected to Shear[END_REF] the principal concrete strain:

T 2D c =      cos 2 θ sin 2 θ sin θ cos θ sin 2 θ cos 2 θ -sin θ cos θ -sin 2θ sin 2θ cos 2θ      (38)
The direction of principal strains can be determined from the in-plane strain vector using Mohr's circle:

θ = 1 2 arctan γ lt ε l -ε t (39) D 2D
prin,c is the concrete secant-stiffness matrix in the principal directions:

D 2D prin,c =      E 1 0 0 0 E 2 0 0 0 E 1 E 2 E 1 + E 2      (40) 17
A c c e p t e d M a n u s c r i p t E 1 and E 2 are the secant moduli, defined as follows:

E 1 = σ 1 ε 1 ; E 2 = σ 2 ε 2 ; ( 41
)
Where ε 1 and ε 2 are the concrete strains in principal directions, obtained from the in-plane concrete strain vector and the direction of principal strains:

e 2D prin,c =      ε 1 ε 2 γ 12      = T 2D c e 2D ip ( 42
)
The concrete stresses in principal directions σ 1 and σ 2 are then conducted from the stress-strain relations expressed below in Section 4.4:

s 2D prin,c = D 2D prin,c e 2D prin,c
The hypothesis of same principal directions between strains and stresses gives:

s 2D ip,c = (T 2D c ) T s 2D prin,c (43)
For the reinforcement's contribution, the steel stiffness matrix is expressed in the in-plane coordinate system:

D 2D ip,s =      ρ sl E sl 0 0 0 ρ st E st 0 0 0 0      (44)
Where E sl is the secant modulus of longitudinal reinforcement and E st the secant modulus of transverse steel which can be determined from the steel behavior law in figure 5; ρ sl and ρ st are the reinforcement ratio in the in-plane directions l, t respectively, which have been determined as follows:

ρ sl = A sl A 2D c (45a) ρ st = A st P st A 2D c s (45b)
where A sl is the total area of longitudinal reinforcing steel situated in 2D-zone,

A 2D c is the area of 2D-zone,
A st is the area of one leg of a transverse steel bar, P st is the perimeter of the stirrup centerline in 2D-zone and s denotes the average spacing of stirrups. It should be noted that the reinforcement percentages are evaluated with respect to the discretized cross-section area. Now all the terms of the in-plane constitutive relation in equation [START_REF] Vecchio | Toward Compression-Field Analysis of Reinforced Concrete Solids[END_REF] are expressed. The process to determine the in-plane transverse strains ε t is however still to be defined. This is done remembering that to satisfy the condition of transverse equilibrium between concrete and stirrup at each material point, the

18

A c c e p t e d M a n u s c r i p t in-plane transverse stress σ t must be zero. From equation ( 35) we have:

s 2D ip = D 2D ip e 2D ip ⇔      σ l σ t τ lt      =      D 2D ip,11 D 2D ip,12 D 2D ip,13 D 2D ip,21 D 2D ip,22 D 2D ip,23 D 2D ip,31 D 2D ip,32 D 2D ip,33           ε l ε t γ lt      (46)
By imposing σ t = 0, one obtains:

ε t = - D 2D ip,21 ε l + D 2D ip,23 γ lt D 2D ip,22 (47)
This equation must be solved by an iterative process, that is presented in Figure 8. After determining the transverse strains and stresses in the in-plane system following the equations above, another change of reference needs to be applied to transform the stress, strain vector and the stiffness matrix to the local frame system:

, l lt   t   2 , D ip c D 2 , s D ip D 2 D ip D 0 t   YES NO   2 D ip l t lt     e

Transverse equilibrium condition

D 2D f = T T ip   D 2D ip 0 3 0 3 0 3   T ip (48a) s 2D f = D 2D f e 2D f ( 48b
)
Where 0 3 is a (3 × 3) zero matrix, resulting from the exclusion of unnecessary strain components in equation

19

A c c e p t e d M a n u s c r i p t [START_REF] Battini | Co-rotational beam elements with warping effects in stability problems[END_REF]. The 2D-stiffness matrix is therefore:

D 2D f -h =               D 2D ip,11 D 2D ip,13 0 D 2D ip,12 0 0 D 2D ip,31 D 2D ip,33 0 D 2D ip,32 0 0 0 0 0 0 0 0 D 2D ip,21 D 2D ip,23 0 D 2D ip,22 0 0 0 0 0 0 0 0 0 0 0 0 0 0               or D 2D f -v =               D 2D ip,11 0 D 2D ip,13 0 D 2D ip,12 0 0 0 0 0 0 0 D 2D ip,31 0 D 2D ip,33 0 D 2D ip,32 0 0 0 0 0 0 0 D 2D ip,21 0 D 2D ip,23 0 D 2D ip,22 0 0 0 0 0 0 0               (49)
After the iteration process of satisfying the internal equilibrium above (Figure 8), the stresses σ yy or σ zz become zero and as a consequence:

     σ yy σ zz τ yz      2D f -h =      D 2D ip,21 D 2D ip,23 0 D 2D ip,22 0 0 0 0 0 0 0 0 0 0 0 0 0 0                    ε xx γ xy γ xz ε yy ε zz γ yz               =      D 2D ip,21 ε xx + D 2D ip,23 γ xy + D 2D ip,22 ε yy 0 0      =      0 0 0      (50a)      σ yy σ zz τ yz      2D f -v =      0 0 0 0 0 0 D 2D ip,21 0 D 2D ip,23 0 D 2D ip,22 0 0 0 0 0 0 0                    ε xx γ xy γ xz ε yy ε zz γ yz               =      0 D 2D ip,21 ε xx + D 2D ip,23 γ xz + D 2D ip,22 ε zz 0      =      0 0 0      (50b)
From this condensation, the strain vector ε yy ε zz γ yz T is expressed in the local frame as:

     ε yy ε zz γ yz      2D f -h = - 1 D 2D ip,22      D 2D ip,21 D 2D ip,23 0 0 0 0 0 0 0           ε xx γ xy γ xz      (51a)      ε yy ε zz γ yz      2D f -v = - 1 D 2D ip,22      0 0 0 D 2D ip,21 0 D 2D ip,23 0 0 0           ε xx γ xy γ xz      (51b)
The value of ε yy or ε zz correspond to the value of ε t determined in equation (47). The stress vector s 2D

s 2D f =      σ xx τ xy τ xz      =      D 2D ip,11 D 2D ip,13 0 D 2D ip,31 D 2D ip,33 0 0 0 0           ε xx γ xy γ xz      +      D 2D ip,12 0 0 D 2D ip,32 0 0 0 0 0           ε yy ε zz γ yz      (52a) s 2D f =      σ xx τ xy τ xz      =      D 2D ip,11 0 D 2D ip,13 0 0 0 D 2D ip,31 0 D 2D ip,33           ε xx γ xy γ xz      +      0 D 2D ip,12 0 0 0 0 0 D 2D ip,32 0           ε yy ε zz γ yz      (52b)
Using the expression of equation ( 51), we get the material constitutive relation in 2D-zone:

s 2D f = k 2D f e 2D f ( 53
)
with the expression of material stiffness matrix:

k 2D f =      D 2D ip,11 D 2D ip,13 0 D 2D ip,31 D 2D ip,33 0 0 0 0      - 1 D 2D ip,22      D 2D ip,12 D 2D ip,21 D 2D ip,12 D 2D ip,23 0 D 2D ip,21 D 2D ip,32 D 2D ip,23 D 2D ip,32 0 0 0 0      (54a) k 2D f =      D 2D ip,11 0 D 2D ip,13 0 0 0 D 2D ip,31 0 D 2D ip,33      - 1 D 2D ip,22      D 2D ip,12 D 2D ip,21 0 D 2D ip,12 D 2D ip,23 0 0 0 D 2D ip,21 D 2D ip,32 0 D 2D ip,23 D 2D ip,32      (54b)
Using equation ( 17), the 2D-regional stiffness matrix K 2D s can then be determined.

3D-Zone

This zone corresponds to the regions of concrete in which transverse steels come across in two directions (the four corners of section) and the regions of concrete in the core of section without any reinforcement (Figure 4). The stress state of this zone has three components: are also related by a material stiffness matrix D 3D f :

s 3D f = D 3D f e 3D f (55)

21

A c c e p t e d M a n u s c r i p t

This stiffness matrix D 3D f is determined by a secant-stiffness-based formulation, which is extended from the constitutive relations of the MCFT for membrane members in section 4.2. After cracking the contribution of concrete and reinforcement can always be added separately to the stiffness matrix:

D 3D f = D 3D f,c + D 3D f,s (56)
Where D 3D f,c is the concrete stiffness and D 3D f,s is the stiffness from reinforcement. As before, while the steel stiffness can be set up easily in the local frame system, for concrete stiffness matrix, it is required a change of reference from principal directions to local axes. Equation ( 56) becomes:

D 3D f = T T c D 3D prin,c T c + D 3D f,s (57)
Where D 3D prin,c is the concrete stiffness matrix evaluated in the principal direction:

D 3D prin,c =                 E 1 0 0 0 0 0 0 E 2 0 0 0 0 0 0 E 3 0 0 0 0 0 0 E 1 E 2 E 1 + E 2 0 0 0 0 0 0 E 2 E 3 E 2 + E 3 0 0 0 0 0 0 E 1 E 3 E 1 + E 3                 (58)
E 1 , E 2 and E 3 are the secant moduli, defined as follows:

E 1 = σ 1 ε 1 ; E 2 = σ 2 ε 2 ; E 3 = σ 3 ε 3 ; ( 59
)
T c is the transformation matrix composed of the direction cosines which define the direction of the principal concrete strains (Figure 9):

T c =               l 2 1 l 1 m 1 m 1 n 1 m 2 1 n 2 1 n 1 l 1 l 2 2 l 2 m 2 m 2 n 2 m 2 2 n 2 2 n 2 l 2 l 2 3 l 3 m 2 m 3 n 2 m 2 3 n 2 3 n 3 l 3 2l 1 l 2 l 1 m 2 + l 2 m1 m 1 n 2 + m 2 n 1 2m 1 m 2 2n 1 n 2 n 1 l 2 + n 2 l 1 2l 2 l 3 l 2 m 3 + l 3 m2 m 2 n 3 + m 3 n 2 2m 2 m 3 2n 2 n 3 n 2 l 3 + n 3 l 2 2l 3 l 1 l 3 m 1 + l 1 m3 m 3 n 1 + m 1 n 3 2m 3 m 1 2n 3 n 1 n 3 l 1 + n 1 l 3               (60)
The concrete strains in principal direction the previous iteration. The concrete stresses in principal directions σ 1 , σ 2 , σ 3 are then deduced from the uniaxial stress-strain relations in Figure 11. According to Vecchio & Selby [START_REF] Vecchio | Toward Compression-Field Analysis of Reinforced Concrete Solids[END_REF], in direction 3 (the direction of principal compressive strain), the compressive stress σ 3 depends on the strain compressive strain ε 3 and the tensile strain ε 1 : σ 3 = f (ε 3 , ε 1 ). In direction 1 (the direction of principal tensile strain), the tensile stress σ 1 depends only on the tensile strain ε 1 : σ 1 = f (ε 1 ). Finally, in the intermediate direction 2 (that can be compressive or tensile), the stress σ 2 depends on ε 2 and ε 1 :

ε 1 , ε 2 , ε 3 (ε 1 > ε 2 > ε 3 )
σ 2 = f (ε 2 , ε 1 ).
Then, for the reinforcement's contribution, the steel stiffness matrix is expressed in the local coordinate system as:

D 3D f,s =               ρ sx E sl 0 0 0 0 0 0 ρ sy E st 0 0 0 0 0 0 ρ sz E st 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0               (61)
Where E sl is the secant modulus of longitudinal reinforcements and E st the secant modulus of transverse steels, which can be determined from the steel behavior law in figure 5; ρ sx , ρ sy and ρ sz are the reinforcement ratio in directions x, y, z respectively, which have been determined as follows:

ρ sx = A sl A 3D c ( 62a
)
ρ sy = A st P y st A 3D c s (62b)
ρ sz = A st P z st A 3D c s (62c)
where A sl is the total area of longitudinal reinforcing steel situated in 3D-zone, A 3D c is the area of 3D-zone,

A st is the area of one leg of a transverse steel bar, P y st is the length of the stirrup centerline disposed in horizontal direction and P z st is the length of the stirrup centerline disposed in vertical direction in 3D-zone; s denotes the average spacing of stirrups. It should be noted that the reinforcement percentages are evaluated with respect to the discretized cross-section area.

23

A c c e p t e d M a n u s c r i p t

Now we have a full material constitutive relation with six-components stress and strain vector and a (6 × 6) stiffness matrix:

s 3D f = D 3D f e 3D f ⇔   s 3D f s 3D of   =   D 3D f,11 D 3D f,12 D 3D f,21 D 3D f,22     e 3D f e 3D of   (63)
where:

s 3D f =      σ xx τ xy τ xz      ; s 3D of =      σ yy σ zz τ yz      ; e 3D f =      ε xx γ xy γ xz      ; e 3D of =      ε yy ε zz γ yz     
As mentioned above, the 3 stresses σ yy , σ zz and τ yz must be zero, so a static condensation needs to be realized:

s 3D of = D 3D f,21 D 3D f,22   e 3D f e 3D of   = 0 3 ⇒ e 3D of = D 3D f,22 -1 D 3D f,21 e 3D f ( 64
)
The non-zero terms e 3D of in the strain vector e 3D f , which can not be determined from the kinematic condition, are calculated from the equation 64. The stress vector s 3D f used in the sectional analysis is therefore expressed by:

s 3D f = D 3D f,11 e 3D f + D 3D f,12 e 3D of = D 3D f,11 e 3D f + D 3D f,12 D 3D f,22 -1 D 3D f,21 e 3D f (65a) ⇒ s 3D f = k 3D f e 3D f (65b)
with the material stiffness matrix:

k 3D f = D 3D f,11 + D 3D f,12 D 3D f,22 -1 D 3D f,21 (66)
Using equation ( 17), the regional stiffness matrix K 3D s can then be determined.

Proposed material model for concrete

In compression, the relationship proposed by Stevens et al. [START_REF] Stevens | Constitutive model for reinforced concrete finite element analysis[END_REF], which is based on the original MCFT, is employed in the model including the softening effect of concrete, which can be understood as the fact that the principal compressive stress is softened due to principal tensile stresses. This phenomenon is represented by a softening parameter SF depending on the tensile strain, which is applied for "shrinking" the stressstrain diagram (Figure 11a). In tension, as mentioned above in the section (1), the tensile relation proposed in the classical MCFT is based on tests of shear panels, so not really suitable for the torsional behavior.

By conducting a systematic parametric study with several experimental tests in torsion, another tensile constitutive law is proposed by Jeng [START_REF] Jeng | Softened membrane model for torsion in reinforced concrete members[END_REF]:

E cr = 5620 f c (f c in MPa) (67a)
ε cr = 0.000116 (67b)

f cr = E cr ε cr (67c)

24

A c c e p t e d M a n u s c r i p t

In this relation, the main difference from the original MCFT is a higher cracking strength of concrete. The cracking strain comes from observations of tests in torsion for plain concrete members by Hsu [START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF], in which he cited that failure of plain concrete members under torsion occur at a strain of approximately 0.0001 [6].

Although this new model for torsion showed some correct responses in the simulation compared to experimental tests, the choice of cracking strain is a bit unnatural to the authors. This is why, based on this idea, it is proposed another tensile constitutive law for torsion. Knowing that the cracking strength of RC f cr under torsion is a function of concrete strength and reinforcement ratio, the authors developed a new relationship between f cr and the cracking strength of plain concrete members under torsion f P C cr , which was widely proposed in literatures, for example in [START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF]. A parametric calibration was studied in 61 specimens of experimental tests, carried out by Hsu [START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF] (41 specimens), Csikos & Hegedus [38] (12 specimens) and Lee et al. [START_REF] Lee | Torsional Strength and Failure Modes of Reinforced Concrete Beams Subjected to Pure Torsion[END_REF] (8 specimens), with various concrete strengths, section dimensions and reinforcement ratios.

In the present model, the longitudinal reinforcement bars only have an impact on the normal stress state, so do not have any contribution to the torsional response, which depends on the shear stresses.

Experimental tests show that, for concrete members with longitudinal steel only, the effect of the longitudinal reinforcements is small, and the cracking torque is very similar to that of a plain concrete member. This ineffectiveness can be explained by the location of longitudinal steel bars, which are usually placed at the corners of the section where the shear stress is zero according to Saint-Venant's stress distribution (Figure 10a). For concrete members with longitudinal steel and stirrups, experimental tests by Hsu [START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF] show that the cracking torque might be a function of the total steel percentage (Figure 10b). Hsu also remarked that a better repartition of reinforcement will somewhat increase the cracking torque: it is the case of beam G6, G7 and G8 in which longitudinal bars are also located at the center of the section. This remark is also noted by the authors during the calibration process, indeed a repartition with additional longitudinal reinforcement steel bars along the perimeter of the cross-section can make influence on the value of cracking torque. For this reason, the reinforcement percentage ρ s must become an essential parameter for the cracking strength f cr , and its formulation should be distinguished for different cases of longitudinal steel repartitions.

For a rectangular section, the calibration study gives the following formulation of f cr :

f cr = f P C
With the expression of f P C cr [START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF]:

f P C cr = 0.76 1 + 6452 b 2 f 1/3 c (MPa) for b > 100 (mm) (69a)
f P C cr = 6.13

f c b 1/3 (MPa) for b ≤ 100 (mm) (69b)
Where f c is the cylinder compressive strength, h the section height, b the section width. The cracking of reinforcement concrete members under torsional effect is therefore depending on both the concrete strength, the reinforcement's contribution and the dimensions of the cross-section. [START_REF] Vecchio | The Modified Compression-Field Theory for Reinforced Concrete Elements Subjected to Shear[END_REF] Vecchio model for shear It is worth to note that this proposed constitutive model is the improved version of the constitutive model proposed by Jeng [START_REF] Jeng | Softened membrane model for torsion in reinforced concrete members[END_REF]. A comparison between two versions will be presented below to validate the pertinence of the new model. the torsional moment which depends on the section dimensions [START_REF] Rahal | Effect of thickness of concrete cover on shear-torsion interaction -an experimental investigation[END_REF]; the formulation used in the model of Valipour & Foster is a function of stirrup's spacing, reinforcement bar disposition and concrete section dimension [START_REF] Valipour | Nonlinear reinforced concrete frame element with torsion[END_REF]; while Hsu proposed another formulation for design based on the torsional strength [START_REF] Hsu | Shear Flow Zone in Torsion of Reinforced Concrete[END_REF].

' c f SF ' c f softening of ' ε c ε c 0 < σ = c SF ' c f ε ( ) 2 - ε ( ) 2 [ ] σ c 0 < 1 0.8-0.34 c SF ε = ' ε c t concrete 0 ' ε c ' ε c concrete strain in compression c c c c c ε t c ε c concrete strain in tension SF ' c f η ε cu (a) Compression relationship
σ c = cr f 0.33 ' c f ( ) MPa σ = c cr f ε > c : ε cr ( ) 500 ε c 1+ = cr E 3925 ' c f ( ) MPa / = cr f cr E ε cr ε cr cr f cr f 0 ε cr ε c 0 > 0 > = cr E 5620 ' c f ( ) MPa / =
Although the section discretization in this model was proposed by Navarro-Gregori et al., these authors did not validate a rule to define the width of the 2D zone, but suggested a width equal to the double of concrete cover over the transverse reinforcement. Taking care of this indication, during the construction of the model the authors realized that a so simple formulation can not represent accurately enough the torsional response, especially in the case of pure torsion. The difference when calculating for ultimate strength are important. After consulting the proposed formulations cited above and a formulation calibrated by Hsu [6], the authors noted that the effective wall thickness should be a function of the section width and of the reinforcement ratio. Another parametric study was also investigated from the 61 specimens below and gives the following formulation:

t
where ρ s is the reinforcement percentage; b, h the section dimension; m is the ratio between longitudinal and transversal reinforcement and m ≤ 1.5

Mechanical model for Plain concrete members

The torsional behavior of plain concrete members can be represented by a typical torque-twist curve as shown in Figure 12a: at low torque, the behavior is linear elastic, then becomes curved at high torque until a brittle failure just after the first crack. The torsional rigidity can be related to the stress-strain relation in uni-axial compression and tension. Knowing this, to describe the behavior of plain concrete members subjected to pure torsion, the constitutive model proposed above can be used with a little improvements in the tensile relation in order to obtain the exact shape of torque-twist curve. For this, the authors proposed to divide the stress-strain curve into 2 straight lines (Figure 12b). For rectangular sections, the cracking stress f cr is calculated from equation (69), while the first slope is taken from the proposed formulation in tension by Jeng [START_REF] Jeng | Softened membrane model for torsion in reinforced concrete members[END_REF].

Elastic torsion with small rotations

In this section, the numerical results calculated by the proposed model will be compared to the theoretical formulations and other model's results. The first example is a cantilever beam subjected to pure torsion at the free end of the beam under small displacement hypothesis (Figure 13), which was also simulated by Le

Corvec in her PhD thesis [START_REF] Corvec | Nonlinear 3d frame element with multi-axial coupling under consideration of local effects[END_REF]. The end twist angle given by the model correlates very well with the one given by the theoretical formulation (Table 1). Moreover, it confirms that neglecting the effect of warping will strongly influence on the twist angle. This result is obtained with a FE model constructed from 4 elements with 5 points Gauss-Lobatto along the element axis and a system of 100 × 50 square fibers using mid-point integration rule over the cross-section. The execution time takes about 9 seconds. In a multi-fiber model, the more number of fibers, better the results obtained, however calculation time becomes higher. In a simple elastic

28

A c c e p t e d M a n u s c r i p t model, the difference is not very important, but in order to apply the model in more complex problems, an adequate size of fiber needs to be fixed. From figure 15a, a number than 30 fibers through the section depth is a good choice to ensure a good balance between accuracy and calculation time. Here, as in the sequel, the relative error is computed with respect to the theoretical (or experimental) values:

RE = V numerical -V theoretical/experimental V theoretical/experimental 100%
The warping profile of the cross-section is represented in Figure 14. It gives a good illustration of the phenomenon. As mentioned in the expression of warping functions, we do not need an infinite (or great) number of items to obtain a satisfying result: in figure 15b, from a number n ≥ 5, the relative difference is less than 0.01%. It is worth to notice that in this multi-fiber finite element model subjected to pure torsion, it is impossible to get an exact result compared to the theoretical formulation, because of two reasons:

1. The warping function is calculated with limited number of items n, instead of infinity.

2. Sectional integration points situate in the center of fiber, not on the border.

Despite this fact, the use of multi-fiber approach is very effective and highly recommended by the authors, considering its good balance between the result obtained (relative error less than 0.5%) and the calculation time.

Elastic torsion with large rotations

To investigate the geometrical non-linearity in pure torsion, a cantilever beam is modeled (Figure 16).

Mx L=100 E=2100000 G=787500 h=10 b=0.5 Compared to the analytical solution based on Vlasov's beam theory and the numerical one of Battini's model, the results obtained by the proposed model show a very good correlation. (Figure 17). As explained in section 3.2, without the Wagner term the result will be purely linear. As the Wagner term expresses the non linearity in torsion, when looking at the analytical and numerical results, we can see that at low torsional moments there is not much difference between the linear and non-linear solution. To further Figure 18a shows that significant differences between linear and nonlinear solutions were always recorded when the torsional moment becomes bigger. However, in a range of twist angle less 1 radian per meter, the relative differences between linear and non-linear solution are very small (< 0.5%) (Figure 18b).

Considering that cracking and failure of concrete (and also reinforced concrete) member under pure torsion happens frequently at a twist rate between 0.1e-3 to 10e-3 radian per meter, the effect of large displacement for concrete member under pure torsion in the elastic material range can be considered as negligible. In this section, the behavior of plain concrete member under pure torsion will be represented and compared to a series of test carried out by Hsu [START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF], Csikos & Hegedus [38] and Lee et al. [START_REF] Lee | Torsional Strength and Failure Modes of Reinforced Concrete Beams Subjected to Pure Torsion[END_REF]. The details of beam sections and material properties can be found in Table 2. For a plain concrete member, without contribution of reinforcement, it's obvious that the whole section is considered as a 3D-zone (section 4.3), but with a constitutive law as proposed in section 5. The beam is modeled with 4 elements, and the cross-section is discretized into a system of 30 × 20 fibers. The torque-twist curves obtained by the proposed model shows a good agreement with the experimental results (Figure 19a). Moreover, in Figure 19b, we can see that the proposed model can predict very well the cracking torsional strength of plain concrete members. In the literature, there are three basic methods to determine the cracking load of plain concrete members under pure torsion: the elastic theory presented by Saint-Venant, the plastic theory developed by Nylander [START_REF] Nylander | Vridning och vridningsinspanning vid betongkonstruktioner (Torsion and torsional restraint by concrete Different types of truss model[END_REF] and the skew-bending theory based on observations of torsion tests. When consulting all the results calculated by these three theoretical formulations and the numerical results, then comparing to the experimental tests (Table 3), we can see that the cracking strength obtained by the elastic theory is always conservative, the plastic theory gives unconservative results, while the results given by skew-bending theory are more reasonable, but generally the proposed numerical model can give the best results: in most cases the relative difference compared to experimental test is less than 5%.

Reinforced concrete members

In this section, the behavior of RC members under pure torsion will be represented and compared to a series of test carried out by Hsu [START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF], Csikos & Hegedus [38], Lee et al. [START_REF] Lee | Torsional Strength and Failure Modes of Reinforced Concrete Beams Subjected to Pure Torsion[END_REF] and Lee et al. [START_REF] Lee | Torsional Behavior of Reinforced Concrete Beams Predicted by a Compatibility-Aided Truss Model[END_REF]. The details of beam section and reinforcement repartition can be found in Figure 20, with two cases of longitudinal reinforcement repartition as mentioned in section 4.

Cracking Torque

For RC members with a rectangular cross-section, Hsu proposed a new formulation for cracking load based on T P C cr :

T cr = (1 + 4ρ s )T P C cr ( 71
)
where ρ s is the total percentage of reinforcement with respect to the whole section. The cracking torque T cr , calculated by three analytical methods as mentioned in section 7.1, are then compared to the numerical results obtained by the proposed model and Jeng's model [START_REF] Jeng | Softened membrane model for torsion in reinforced concrete members[END_REF]. In Table 4, the relative differences (calculated with respect to the experimental tests) are indicated. In each series of tests the smallest average relative error is shown in bold and highlighted in red. We can see that the proposed modifications are very relevant as the smallest average relative error can be found in most cases. The numerical results given by Jeng's model are also satisfactory, although some large average relative errors (more than 10%) were found in series C and N of Hsu's test. When considering the analytical formulations, while the elastic theory always gives the unsatisfactory results, reasonable cracking loads can be sometimes obtained by the plastic theory (series C and N of Hsu's test). The skew-bending theory shows the best ability of predicting the cracking torsional moment between the three formulations, but its results are still far from satisfactory when comparing to the experimental results.

Another strong aspect of the proposed model is the ability of calculating for the case of repartition with additional longitudinal reinforcement bars along the perimeter. In cracking torque for this case of reinforcement repartition are all indicated, with a highlight in bold and red for the model with the smallest relative difference compared to the experimental result. Unconservative results can be found in many cases, except with the plastic theory, but the best agreement is obtained with the proposed model. Although giving the smallest differences in two specimens (G8 and N3), the predictions of the plastic theory are too conservative, which is not suitable for the design process. The first formulation to calculate the ultimate torsional moment was proposed by Rausch in his space truss theory [START_REF] Raush | Design of reinforced concrete in torsion[END_REF]. The main hypotheses are that both longitudinal and transversal steels are yielding at the ultimate torque and the inclination of concrete cracks is fixed at 45˚. This formulation is an ingenious concept which gives torsional strength as a function of reinforcements and concrete:

T uR = min 2A c A sl f sly u ; 2A c A st f sty s (72)
where A c is the area enclosed by the centerline of stirrups; A sl and f sly are the total area and the yield strength of longitudinal bars; A st is the cross-section area of one hoop bar and f slt are the yield strength of stirrups; u is the perimeter of the area bounded by the centerline of a completed hoop bar; s is the stirrup's spacing. Then, Cowan modified Rausch's equation and proposed a new formulation based on Saint-Venant's stress and strain distribution for rectangular cross-section [47]:

T uC = T e + 1.6 A c A st f sty s (73)
where T e is the torsional cracking strength calculated by the elastic theory as follows:

T e = αb 2 hf t
where b is the section width, h is the section height and b ≤ h, α is the Saint-Venant's coefficient depending on the ratio of h/b, f t is the tensile strength of concrete. According to Cowan, the torsional failure mode depends on the yielding of stirrups. Hsu proposed another formulation based on experimental tests' observations [6]:

T uH = b 2 h 3 2.4 f c + √ m f sly f sty 1 + 0.2 b 1 h 1 b 1 h 1 A st f sty s (74)
where f c is the compressive strength of concrete; b 1 and h 1 is the width and the height of the area enclosed by the centerline of stirrups; m is the ratio between longitudinal and transversal steel percentage.

In the ACI standard, torsional strength of RC members is calculated as [49]:

T uS = 2A o A sl f sly s cot θ (75)
where A o is the gross area enclosed by the shear flow path, which can be taken equal to 0.85A c , with A c the area enclosed by the centerline of stirrups; θ is the cracks angle which can be taken as 45˚for non-prestressed or low-prestressed members. In the European standard Eurocode 2, three different values are calculated depending on the torsional failure modes and the minimum one has to be chosen [50]. The first value is

36

A c c e p t e d M a n u s c r i p t related to the stirrups' yielding, the second one corresponds to the longitudinal bars capacity and the third one is related to the torsional capacity of concrete struts:

T uE = min 2A k A sl f sly s cot θ; 2A c A sl f sly u k tan θ; 2νf c A k t ef sin θ cos θ ( 76
)
A k is the area enclosed by the centre-lines of the effective wall thickness; t ef is the effective wall thickness and can be calculated as A/u with A the total area and u the perimeter of cross-section. The cracks angle is variable but can be taken as 45˚.

In A more detailed comparison will be examined for the case of repartition with additional longitudinal reinforcement bars along the perimeter. Similar to section 7.2.1, in Table 7, the numerical and analytical ultimate torque for this case of reinforcement repartition are all indicated, with a highlight in bold and red for the cases with the smallest relative difference compared to the experimental results. In Hsu's test, the two numerical models show their advantage by their ability to predict the torsional strength with a very reasonable error in most cases (except specimen G8), while in the test of Lee et and Relative error with respect to experimental test.

Torque-twsit curve

The torque-twist curve obtained by the proposed model is compared to the experimental results of specimen M2 and M3 in the test of Hsu (Figure 21), with and without the proposed modification of the tensile behavior. A very good agreement is achieved: the cracking and the ultimate torsional moment have the same magnitude, the slopes in the post-cracking phase are identical, and the featured horizontal plateau manifested by the transition between the two phases before and after the cracking is well represented.

Moreover, we can conclude that without the proposed modification of tensile behavior, the cracking moment achieved is about half of the one of the experimental test, and therefore the torque twist curve can not be similar to the experimental one.

Conclusion

In this present work, a non-linear multi-fiber finite element model is developed for RC element under pure torsion. In the elastic domain, the comparison of numerical results with the analytical formulation is excellent. The warping phenomenon is fully represented at sectional level as well as its influence on element deformations. Large rotation conditions have also been investigated and included successfully in the proposed numerical model for RC members subjected to torsion.

Regarding the research significance, the modification proposed by the authors for the tensile behavior of concrete and the parametric formulation developed by the authors for effective wall thickness shows a reasonable and correct influence in predicting the torsional response in the inelastic domain. Although the model is developed with a displacement-based formulation, the sectional model proposed is general enough to apply in any finite element model or formulation, such as force-based or flexibility-based. Finally, with the proposed model, the torsional response can be obtained in a simple manner with reasonable computational efficiency.

Moreover, the same sectional discretization and constitutive material approach has to be applied in a more general model, which is able to consider the effect of combined loading shear-bending-torsion. For future studies, this model may be further developed with a larger calibration study and other new formulations for different reinforcement repartitions may be proposed.

D 2D f = T T ip   D 2D ip 0 3 0 3 0 3   T ip (A.1a) s 2D f = T T ip s 2D ip (A.1b) s 2D f = D 2D f e 2D f (A.1c)
Where 0 3 is a (3 × 3) zero matrix, resulting from the exclusion of unnecessary strain components in equation

D 2D f -h =               D 2D ip,11 D 2D ip,13 0 D 2D ip,12 0 0 D 2D ip,31 D 2D ip,33 0 D 2D ip,32 0 0 0 0 0 0 0 0 D 2D ip,21 D 2D ip,23 0 D 2D ip,22 0 0 0 0 0 0 0 0 0 0 0 0 0 0               or D 2D f -v =               D 2D ip,11 0 D 2D ip,13 0 D 2D ip,12 0 0 0 0 0 0 0 D 2D ip,31 0 D 2D ip,33 0 D 2D ip,32 0 0 0 0 0 0 0 D 2D ip,21 0 D 2D ip,23 0 D 2D ip,22 0 0 0 0 0 0 0               (A.2)
and the full stress vector:

s 2D f -h = σ xx τ xy 0 σ yy 0 0 T = σ l τ lt 0 σ t 0 0 T s 2D f -v = σ xx 0 τ xz 0 σ zz 0 T = σ l 0 τ lt 0 σ t 0 T (A.3)
Knowing that the value of σ t must be zero to satisfy the internal equilibrium, the following static condensation can be established with the aid of equation (A.1c), for the 2D-horizontal zone first:

     σ yy 0 0      2D f -h =      0 0 0      =      D 2D ip,21 ε xx + D 2D ip,23 γ xy + D 2D ip,22 ε yy 0 0      ⇒      ε yy 0 0      2D f -h = - 1 D 2D ip,22      D 2D ip,21 D 2D ip,23 0 0 0 0 0 0 0           ε xx γ xy 0      (A.4)
and for the 2D-vertical zone:

     0 σ zz 0      2D f -v =      0 0 0      =      0 D 2D ip,21 ε xx + D 2D ip,23 γ xz + D 2D ip,22 ε zz 0      ⇒      0 ε zz 0      2D f -v = - 1 D 2D ip,22      0 0 0 D 2D ip,21 0 D 2D ip,23 0 0 0           ε xx 0 γ xz      (A.5)
Then, from equation (A.1c), (A.3), (A.4) and (A.5), the stress vector s 2D f used in the sectional analysis can

46

A c c e p t e d M a n u s c r i p t now be expressed as follows:

s 2D f -h =      σ xx τ xy 0      =      D 2D ip,11 D 2D ip,13 0 D 2D ip,31 D 2D ip,33 0 0 0 0           ε xx γ xy 0      +      D 2D ip,12 0 0 D 2D ip,32 0 0 0 0 0           ε yy 0 0      (A.6a) s 2D f -v =      σ xx 0 τ xz      =      D 2D ip,11 0 D 2D ip,13 0 0 0 D 2D ip,31 0 D 2D ip,33           ε xx 0 γ xz      +      0 D 2D ip,12 0 0 0 0 0 D 2D ip,32 0           0 ε zz 0      (A.6b)
From equation (A.4), (A.5) and (A.6), the following material constitutive relation can be obtained in 2Dzone:

s 2D f = k 2D f e 2D f (A.7)
with the expression of the material stiffness matrix: • Cracking strength of plain concrete members under torsion f P C cr .

k 2D f -h =      D 2D ip,11 D 2D ip,13 0 D 2D ip,31 D 2D ip,33 0 0 0 0      - 1 D 2D ip,22      D 2D ip,12 D 2D ip,21 D 2D ip,12 D 2D ip,23 0 D 2D ip,21 D 2D ip,32 D 2D ip,23 D 2D ip,32 0 0 0 0      (A.8a) k 2D f -v =      D 2D ip,11 0 D 2D ip,13 0 0 0 D 2D ip,31 0 D 2D ip,33      - 1 D 2D ip,22      D 2D ip,
• Reinforcement percentages (longitudinal, transversal steels and total ratios: ρ sl , ρ st and ρ s ).

• The ratio between ρ sl and ρ st : m = ρ sl ρ st .

• The ratio between section height and width: h b .

The concrete tensile strength f cr was determined in the phase before cracking of concrete, knowing that in this phase the torsional behavior of RC member is similar to that of a corresponding plain concrete member. The cross-section discretization contains therefore only the 3D-zone of concrete fiber. The determination of the effective wall thickness t e is carried out in the phase after cracking of concrete, knowing that the torsional behavior of RC member in this phase depends strongly on this parameter. In this phase, the discretization of cross-section contains 3 zones as described in section 4. The same process has been applied to the determination of the effective wall thickness t e , but this time the value of ultimate cracking torque T u is used to evaluate the value exact of t e , the relative difference between experimental and numerical values is therefore: Appendix C. Section dimension, reinforcement arrangement and material properties for experimental test

In Hsu test, the concrete cover thickness is 19 mm for all specimens, while the materials properties and reinforcement's dimensions are varied and cited as follows (

51

A c c e p t e d M a n u s c r i p t

In the test of Lee et al. [START_REF] Lee | Torsional Strength and Failure Modes of Reinforced Concrete Beams Subjected to Pure Torsion[END_REF], the concrete cover thickness is 20 mm for all specimens, while the compressive strength of concrete f c = 42.6 (MPa), the material properties and steel's dimension are cited as in Table C.9). In the test of Csikos & Hegedus [38], the concrete cover thickness is 15 mm for all specimens, while the concrete used was C20 and the yielding strength of steel is 240 MPa in all cases (

4

 4

Figure 1

 1 Figure 1 Ritter's and Mörsch's original truss model

Figure 2

 2 Figure 2 Mörsch's truss analogy model However, experimental studies carried out in Stuttgart during 1960-s [3] indicated that the stresses in shear reinforcement were considerably lower than those predicted by the truss analogy model. This is due to the contribution of other components to the shear carrying mechanism, among which the most significant are: contribution of concrete in the compression zone, aggregate interlock along inclined cracks and dowel action of the longitudinal reinforcement crossing the

Figure 1 :

 1 Figure 1: Different types of truss model. In order to generalize the fixed-angle truss model, Lampert & Thurlimann developed a model in which the fixed inclined angle of 45˚was derived, called the variable-angle truss model [9]. This model can be used for elements subjected to torsion as well as to a combination of torsion and bending, with the application of the theory of plasticity for RC members. Elfgren extended the use of variable-angle truss model to members

Figure 2 :

 2 Figure 2: Warping for several non circular section [4].

Figure 3 :

 3 Figure 3: Multifiber approach for a RC member and local element frame coordinate

Figure 4 :

 4 Figure4: Discretization of cross-section following the material stress state[START_REF] Navarro-Gregori | A 3D numerical model for reinforced and prestressed concrete elements subjected to combined axial, bending, shear and torsion loading[END_REF]

15

 15

Figure 6 :

 6 Figure 6: In-plane frame system

Figure 7 :

 7 Figure7: Original Modified Compression Field Theory for reinforced concrete members[START_REF] Vecchio | The Modified Compression-Field Theory for Reinforced Concrete Elements Subjected to Shear[END_REF]

Figure 8 :

 8 Figure 8: Iteration process satisfying internal equilibrium between concrete and transverse reinforcement .

f 20 A

 20 c c e p t e d M a n u s c r i p t used in the sectional analysis is now expressed by:

  one normal σ xx and two transverses τ xy and τ xz , included in a stress vector s 3D f = σ xx τ xy τ xz T . Same as in the 2D-zone, in the 3D-zone even though only 3 stress components are considered in the sectional analysis, the fully stress vector has a total of 6 components s 3D f = σ xx τ xy τ xz σ yy σ zz τ yz T . The three components σ yy , σ zz and τ yz will then be set to zero in order to satisfy the internal equilibrium. The corresponding strain vector has therefore six components: e 3D f = ε xx γ xy γ xz ε yy ε zz γ yz T . Unlike the stress vector, three strains ε yy , ε zz and γ yz are not zero, but they are not taken into account in the sectional analysis, as the distortion phenomenon is not considered in this work. Similar to the constitutive relation between s 3D f

22 A

 22 and the direction cosines are calculated from the eigenvalues and eigenvector of the strain tensor, in which three strains (ε xx ; γ xy ; γ xz ) are obtained from the kinematic condition in equation (8), while three others (ε yy ; ε zz ; γ yz ) are taken from

Figure 9 :

 9 Figure 9: Angle between principal direction i and the local coordinate system.

  Distribution of Saint-Venant stresses along edges, the shear stress at each point is calculated by τ 2 xy + τ 2 xz [6]. (b) Cracking torque as a function of reinforcement ratio inHsu tests (1968)[START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF]

Figure 10 25A

 10 Figure 1025

cr 0.38 b h ρ s + 1 ( 3 (

 13 For usual/conventional repartition of reinforcement bars at corners) (68a) f cr = f For repartition with additional reinforcement bars along the perimeter)

  Proposed tensile relation for torsion compared to the shear model of Vecchio[START_REF] Vecchio | The Modified Compression-Field Theory for Reinforced Concrete Elements Subjected to Shear[END_REF]

Figure 11 :

 11 Figure 11: Concrete constitutive relations for proposed model based on the MCFT

4. 5 .

 5 Effective wall thickness As mentioned above, determining an effective wall thickness becomes an obvious question to investigate the post-cracking behavior of RC members under torsional effect. Several formulations were proposed in the literature: Rahal & Collins proposed an average value of the effective thickness of concrete in resisting 26 A c c e p t e d M a n u s c r i p t

  e = b 0.013 h b mρ s + 0.1 (For usual/conventional repartition of reinforcement bars at corners) (70a) t e = b 0.0088 h b mρ s + 0.1 (For repartition with additional reinforcement bars along the perimeter)

27

 27

Figure 12 :

 12 Figure 12: Behavior of plain concrete member subjected to pure torsion

Figure 13 :

 13 Figure 13: Cantilever beam subjected to pure torsion.

Figure 14 :

 14 Figure 14: Warping of cross-section under torsional effect.

29 A

 29 Necessary of items in warping function

Figure 15

 15 Figure 15

Figure 16 :

 16 Figure 16: Cantilever beam subjected to pure torsion.

30

 30

A c c e p t e d M a n u s c r i p t

  investigate this observation, we have tried another example of a concrete cantilever beam, with 1000 mm length, section dimension of 200 mm × 100 mm and elastic modulus of 25000 MPa, subjected to pure torsion, in order to view when the difference becomes not negligible. The torsional moment -twist rate diagram in

Figure 17 :

 17 Figure 17: Elastic torsional response under large displacements hypothesis.

5 (

 5 b) Relative difference between two solutions.

Figure 18 :

 18 Figure 18: Investigation of relative difference between linear and nonlinear geometric conditions of beam under pure torsion.

31

 31

  Torque-twist curve for series A of Hsu's test. Experimental and numerical cracking strength.

Figure 19 32A

 19 Figure 19

4 : 1 . 2 . 4 ,

 4124 Usual/conventional repartition of reinforcement bars at corners: Series B (B1 to B10), C (C1 to C6), G (G1 to G5), I (I2 to I6), M (M1 to M5), N (N1, N2) of Hsu's test; Series B (B1 to B3), C (C1 to C3) of Csikos & Hegedus's test; Series T (Beam T1-1, T2-1, T2-2) of Lee et al. test. Repartition with additional reinforcement bars along the perimeter: Series T (Beam T1-2, T1-3, T1-T2-3, T2-4) of Lee et al. test; Beam M6, N3, N4, G6, G7, G8 of Hsu's test; Series D (D1 to D3), E (E1 to E3) of Csikos & Hegedus's test.The cracking and ultimate torsional moment calculated for 61 specimens of tests cited above are compared to theoretical formulations, numerical results of Jeng's model and 2 international standards for torsional design of RC members (Eurocode 2 and ACI).

33

 33

Figure 20 :

 20 Figure 20: Section dimension in torsional tests [36] [38] [39] .

38

 38

Figure 21

 21 Figure 21Another comparison of the torque-twist curves is made for series B of Hsu's test and series T1, T2 of Lee et al.'s test: once again the proposed model gives a agreement with the experimental results (Figure22, 23). In these figures, it can be seen that the ultimate strength and the post cracking torsional stiffness are very strongly a function of steel percentage. The cracking torsional moments are also different in each specimen and depend on the reinforcement ratio.

Figure 22 : 39 A

 2239 Figure 22: Torque-twist curve for series B of Hsu's test.

Figure 23 :

 23 Figure 23: Torque-twist curve for series T of Lee et al.'s test.

40

 40

B. 24

 24 Linear fitting in the calibration process of determining f cr , with b h ρ s in the horizontal axis. . 48 B.25 Linear fitting in the calibration process of determining f cr , with h b mρ s in the horizontal axis. 49 B.26 Linear fitting in the calibration process of determining f cr , with h b ρ st in the horizontal axis. . 49 B.27 Linear fitting in the calibration process of determining t e . . . . . . . . . . . . . . . . . . . . . 50 List of Tables 1 Mid span twist angle for cantilever beam under mid span torque . . . . . . . . . . . . . . . . 29 44 A c c e p t e d M a n u s c r i p t Appendix A. Development of material stiffness matrix in 2D-zone

( 33 )

 33 ; s 2D ip = s 2DT ip 0 0 0 T is the full in-plane stress vector. The 2D material stiffness matrix in 45 A c c e p t e d M a n u s c r i p t equation (A.1a) can be expressed as follows:

  All the experimental data, including the beam dimensions, disposition and dimension of reinforcement steel, material properties of concrete and steel, experimental values of cracking and ultimate torque, etc., have been collected for the calibration process. Then, according to the observation in section 4.4, the following variables have been chosen for the parametric formulation of determining the concrete tensile strength under torsion f cr and the effective wall thickness t e :

  47

A c c e p t e d M a n u s c r i p t

  In the first time, for each specimen of test, an initial and arbitrary value of f cr (f o cr ) is chosen to evaluate the initial value of cracking torque T cr . Then, considering the experimental values, this initial value f o cr is modified/calibrated in order to obtain the new values of T cr . This step must be performed until the values of T cr are close to the experimental ones, in other word the relative difference between experimental and numerical values are as small as possible:Rcr = T cr -T exp cr T exp cr 100% The value of f cr at this moment is called as calibrated value f c cr . Next, all the values of f c cr and the variables mentioned above are plotted in a diagram in order to determine a reasonable relationship between them. The proportion of f c cr and f P C cr is put in the vertical axis, while in the horizontal axis, several possibilities have been tested, such as h b ρ s , h b ρ st or b h mρ s , etc., in order to select the best trend fitting for two set of data. For the shake of simplicity, a linear fitting is applied for the data set, and the best result has been selected between various possibilities in the horizontal axis (Figure B.24, B.25 and B.26). The best result can be found for the set of b h ρ s as in equation (68), and the final expression could be therefore deduced. Repartition with additional longitudinal reinforcement along the perimeter.

Figure B. 24 :

 24 Figure B.24: Linear fitting in the calibration process of determining fcr, with b h ρs in the horizontal axis.

Figure B. 25 :

 25 Figure B.25: Linear fitting in the calibration process of determining fcr, with h b mρs in the horizontal axis.

Figure B. 26 :

 26 Figure B.26: Linear fitting in the calibration process of determining fcr, with h b ρst in the horizontal axis.

A

  is also applied for the data set, with the proportion of t e and the section width b in the vertical axis and the best fitting set values ( h b mρ s as in equation (70)) in the horizontal axis (Figure B.27), and the final expression could be deduced. 49 Repartition with additional longitudinal reinforcement along the perimeter.

Figure B. 27 :

 27 Figure B.27: Linear fitting in the calibration process of determining te.

Table 1 :

 1 Mid span twist angle for cantilever beam under mid span torque

Table 2 :

 2 Details of specimen of tests in torsion

	Test	Beam Dimension (mm) h/b ratio	f c (MPa)
		A2	381 × 254	1.5	31.2
	Hsu [36]	A4	254 × 254	1	32.49
		A8	280 × 152	1.84	31.2
	Csikos &	A-1	130 × 130	1	23.7
	Hegedus [38]	A-2	130 × 130	1	30
	Lee et al. [39]	T0	350 × 300	1.17	42.6

Table 3 :

 3 Cracking torsional model for Plain concrete member: Comparison between experimental, numerical result and theo-

	retical formulations

Table 5

 5

	, the numerical and analytical
	34

Table 4 :

 4 Cracking torsional moment: average relative error with respect to experimental test calculated.

Table 5 :

 5 Repartition with additional longitudinal reinforcement bars along the perimeter: Cracking torsional moment (KNm) and Relative error with respect to experimental test.In this section, the ultimate torque calculated by the proposed model is compared to the numerical results of Jeng's model, the analytical formulations proposed by Rausch, Cowan and Hsu and two standards for torsional design: Eurocode 2 and ACI.

					a n u s c r i p t
	Test	Beam G6	EXP 30.96	Proposed Model 33.2 / 7.2% M Jeng's Model -	Elastic Theory 18.47 / -40.3%	Plastic Theory 33.61 / 8.6%	Skew-bending Theory 25.27 / -18.4%
		G7	33.67	34.2 / 1.6%	31.24 / -7.2%	19.53 / -42.0%	35.53 / 5.5%	26.21 / -22.2%
	Hsu [36] Lee et al. [39] Csikos & Hegedus [38] A c c e p t e d G8 33.67 32.6 / -3.2% 30.32 / -9.9% 18.53 / -45.0% M6 22.71 21.9 / -3.9% 21.81 / -4% 15.83 / -30.3% N3 7.41 6.80 / -8.3% 6.37 / -14.1% 3.80 / -48.7% N4 7.60 7.50 / -1.4% 6.45 / -15.1% 3.94 / -48.1% T1-2 31.4 30.1 / -4.1% -22.95 / -26.9% T1-3 31.8 31.5 / -0.9% -23.52 / -26.0% T1-4 33.7 34.6 / 2.7% -24.42 / -27.5% T2-3 29.8 29.0 / -2.7% -23.36 / -21.6% T2-4 29.9 29.5 / -1.3% -23.75 / -20.6% D2 1.69 1.70 / 0.6% -1.49 / -11.6% D3 1.69 1.68 / -0.6% -1.47 / -13.0% E1 2.25 2.09 / -7.1% -1.55 / -30.9% E3 2.06 2.05 / -0.5% -1.52 / -26.1%	33.72 / 0.1% 26.18 / -22.3% 26.76 / 17.8% 21.79 / -4% 6.91 / -6.7% 6.36 / -14.1% 7.18 / -5.6% 6.61 / -13.1% 39.03 / 24.3% 27.59 / -12.1% 40. 1/ 25.8% 28.28 / -11.1% 41.53 / 23.2% 29.35 / -12.9% 39.73 / 33.3% 28.08 / -5.8% 40.40 / 35.1% 28.55 / -4.5% 2.40/ 41.7% 2.31 / 36.6% 2.36 / 39.4% 2.29 / 35.6% 2.49 / 10.8% 2.35 / 4.2% 2.44 / 18.4% 2.33 / 12.9%
					35

Table 6

 6 , the relative differences (calculated with respect to the experimental tests) are indicated. In each series of tests, the smallest average relative error is shown in bold and highlighted in red. Similar to the cracking torque, in most cases the smallest average relative error was obtained by the proposed model. The numerical results given by Jeng's model are also satisfactory with the best average relative error in series M and N of Hsu's test. When considering the analytical formulations, Rausch's and Cowan's formulations often give too high average relative error in all cases, while Hsu's formulation performs better, with one best result in series C of Hsu's test. Finally, the torsional strength given by the design standards are far from satisfactory.

	Test	Beams	Proposed	Jeng's	Rausch's	Cowan's	Hsu's	Eurocode 2	ACI
			Model	Model	Formulation	Formulation	Formulation
		Series B	4.9%	7.9%	56.0%	54.8%	16.1%	23.2%	41.4%
		Series C	12.5%	5.6%	81.3%	69.9%	5.1%	14.4%	52.2%
		Series G	3.8%	4.8%	11.1%	21.0%	20.5%	29.5%	10.5%
	Hsu [36]
		Series I	2.8%	6.1%	18.0%	29.9%	24.1%	14.8 %	9.6 %
		Series M	5.9%	3.2%	22.6%	18.7%	17.1%	20.7%	19.2%
		Series N	2.9%	2.5%	12.0%	24.1%	22.1%	22.2%	10.8%
	Lee et al. [39]	Series T	5.5%	-	38.9%	60.2%	5.6%	17.7%	23.2%
	Csikos &	Series B,C	3.7%	-	20.2%	35.5%	30.8%	46.0%	27.2%
	Hegedus [38]	Series D,E	7.0%	-	36.1%	8.9%	34.7%	56.9%	45.6 %
	Average		5.4%	5.1%	32.9%	35.9%	19.5%	27.3%	26.6%

Table 6 :

 6 Ultimate torsional moment: average relative error with respect to experimental test calculated.

  al., Hsu's formulation also gives very good results. In the test of Csikos & Hegedus, the proposed model is still giving the best

	predictions.
	Test	Beams	EXP	Proposed	Jeng's	Hsu's	Eurocode 2	ACI
				Model	Model	Formulation
		G6	39.09	40.30 / 3.1%	-	30.52 / -21.9%	27.98 / -28.4%	33.68 / -13.8%
	Hsu [36] Lee et al. [39] Csikos & Hegedus [38]	G7 G8 M6 N3 N4 T1-2 T1-3 T1-4 T2-3 T2-4 D2 D3 E1 E3	52.65 73.44 60.11 12.20 15.70 42.9 54.1 62.4 50.2 56.4 2.25 2.06 3.38 3.23	50.3 / -4.5% 61.4 / -16.4% 60.10 / 0.1% 12.05 / -1.2% 15.75 / 0.3% 44.5 / 3.7% 52.5 / -3.0% 63.4 / 1.6% 49.4 / -1.6% 55.2 / -2.1% 2.38 / 6.7% 2.20 / 6.8% 3.10 / -9.5% 3.06 / -5.3%	55.87 / 6.1% 70.49 / -4.0% 55.29 / -3.2% 12.49 / 2.4% 15.08 / -4.0% --------a n u s c r i p t 41.07 / -22.0% 40.14 / -23.8% 48.32 / -8.2% 54.26 / -26.1% 57.02 / -22.4% 68.64 / -6.5% 60.85 / 1.2% 64.22 / 6.8% 77.82 / 29.5% 9.37 /-23.2% 9.27 / -24.0% 11.14 / -8.7% 12.99 / -17.3% 14.76 / -6.0% 17.75 / 13.0% 41.86 / -2.4% 47.52 / 10.8% 55.16 / 28.6% 51.84 / -4.2% 62.15 / 14.9% 72.14 / 33.3% 79.51 / 27.42% 92.29 / 47.91% 107.13 / 71.68% 49.95 / -0.5% 45.90 / -8.6% 53.28 / 6.1% 56.75 / 0.6% 53.86 / -4.5% 62.52 / 10.9% 1.51 / -33.1% 0.99 / -55.9% 1.02 / -54.8% 1.50 / -27.3% 0.99 / -51.8% 1.02 / -50.7% 2.01 / -4.0% 1.32 / -60.9% 2.03 / -39.9% -2.00 / -38.0% 1.32 / -59.1% 2.03 / -37.1%
					M
	A c c e p t e d
					37

Table 7 :

 7 Repartition with additional longitudinal reinforcement bars along the perimeters: Ultimate torsional moment (KNm)

Table C

 C

	Beam	fc (MPa)	Longitudinal rebar (mm)	f sl (MPa)	(mm) and (mm) Stirrup	fst (MPa)
	B1	27.58	12.7	314	9.5 at 152	341
	B2	28.61	15.9	316	12.7 at 181	320
	B3	28.06	19.1	328	12.7 at 127	320
	B4 B5 B6 B7 B8 B9 B10 C1 C2 C3 C4 C5 C6 G1 G2 G3 G4 G5 G6	30.54 29.03 28.82 25.99 26.75 28.82 26.48 26.34 26.54 26.89 26.48 27.23 27.58 26.34 26.54 26.89 26.48 27.23 27.58	.8): 12.7 at 92 12.7 at 70 12.7 at 57 12.7 at 127 12.7 at 57 9.5 at 152 9.5 at 152 9.5 at 216 9.5 at 117 12.7 at 140 12.7 at 98 12.7 at 73 12.7 at 54 9.5 at 216 9.5 at 117 a n u s c r i p t 320 323 332 321 332 323 320 319 322 320 319 343 3334 342 341 341 334 345 331 330 336 328 328 329 316 328 341 341 334 345 331 12.7 at 140 330 336 12.7 at 98 328 328 12.7 at 73 329 M 22.2 25.4 28.7 12.7 12.7 19.1 28.7 9.5 12.7 15.9 19 22.2 25.4 9.5 12.7 15.9 19 22.2 25.4 316 12.7 at 54 328
	G7 G8 N1a N2 N2a N3 N4 M2 M3 M4 M5 M6 I1 A c c e p t e d 27.23 22.2 27.58 25.4 26.34 9.5 26.54 12.7 26.89 15.9 26.48 19 27.23 22.2 27.58 25.4 26.34 9.5 26.54 12.7 26.89 15.9 26.48 19 27.23 22.2 27.58 25.4 26.34 9.5 I2 26.54 12.7	328 316 341 334 331 336 328 316 341 334 331 336 328 316 341 334	12.7 at 73 12.7 at 54 9.5 at 216 9.5 at 117 12.7 at 140 12.7 at 98 12.7 at 73 12.7 at 54 9.5 at 216 9.5 at 117 12.7 at 140 12.7 at 98 12.7 at 73 12.7 at 54 9.5 at 216 9.5 at 117	329 328 341 345 330 328 329 328 341 345 330 328 329 328 341 345
	I3	26.89	15.9	331	12.7 at 140	330
	I4	26.48	19	336	12.7 at 98	328
	I5	27.23	22.2	328	12.7 at 73	329
	I6	27.58	25.4	316	12.7 at 54	328
	J1	26.34	9.5	50 341	9.5 at 216	341
	J2	26.54	12.7	334	9.5 at 117	345
	J3	26.89	15.9	331	12.7 at 140	330
	J4	26.48	19	336	12.7 at 98	328

Table C

 C

.8: Details of specimen of torsion tests carried by Hsu

[START_REF] Hsu | Torsion of Structural Concrete-Plain Cocnrete Rectangular Sections[END_REF]

.

Table C .

 C 9: Details of specimen of torsion tests carried by Lee et al.[START_REF] Lee | Torsional Strength and Failure Modes of Reinforced Concrete Beams Subjected to Pure Torsion[END_REF].

	Beam	Longitudinal rebar (mm)	f sl (MPa)	(mm) and (mm) Stirrup	fst (MPa)
	T1-1	4 × 12.7	410	9.5 at 130	370
	T1-2	6 × 12.7	410	9.5 at 85	370
	T1-3	8 × 12.7	410	9.5 at 65	370
	T1-4	6 × 15.9	510	12.7 at 75	355
	T2-1	4 × 12.7	410	9.5 at 225	370
	T2-2	4 × 15.9	510	9.5 at 130	370
	T2-3	6 × 15.9	510	9.5 at 88	370
	T2-4	2 × 12.7	512.4	9.5 at 75	370
		4 × 19

Table C .

 C Table C.10). 10: Details of specimen of torsion tests carried by Csikos & Hegedus. [38].

	Series	Longitudinal	Stirrup
		rebar (mm)	(mm) and (mm)
	B	4 × 6	6 at 130
	C	4 × 6	6 at 65
	D	8 × 6	6 at 130
	E	8 × 6	6 at 65

Appendix D. Experimental values and relative error for experimental tests

[47] Cowan, H. J., An elastic theory for the torsional strength of rectangular reinforced concrete beams, Magazine of Concrete Research 2 (1950) 3-8.

[48] ACI Committee 318, Building code requirements for reinforced concrete (ACI 318-71), Detroit, Mich. : American Concrete Institute. Reinforced concrete -Standards; Reinforced concrete construction -Standards (1971).

[49] ACI Committee 318, Building code requirements for reinforced concrete (ACI 318-95) and commentary (ACI318 R-95), Detroit, Mich. : American Concrete Institute. Reinforced concrete -Standards; Reinforced concrete construction -Standards (1995).

[50] European Committee, Eurocode 2: Design of concrete structures, Reinforced concrete -Standards; Reinforced concrete construction -Standards (1992).

57

A c c e p t e d M a n u s c r i p t