A 3 D mixed frame element with multi-axial coupling for thin-walled structures with damage

A 3D mixed beam finite element is presented, modeling the warping of the cross-sections as an independent kinematic field. The beam formulation is derived on the basis of the Hu-Washizu variational principle, expressed as function of four independent fields: the standard displacements, strains and stresses and the additional warping displacement. This is interpolated along the beam axis and on the cross-section, by placing on it a regular grid of interpolation points and adopting Lagrange polynomials. The warping degrees of freedom defined at the cross-section interpolation points are condensed, thus preserving the element matrix and vector sizes. A fiber discretization of the cross-sections is adopted. The constitutive relationship at the midpoint of each fiber is based on an isotropic damage model for brittle-like materials, distinguishing between the damage variables in tension and in compression to properly describe the unilateral effect. An efficient algorithm is formulated for the element state determination, based on a consistent linearization of the governing equations. A simple numerical application on a cantilever beam with torsion in the linear elastic range is presented and two torsion tests on plain concrete beams are performed, by comparing the numerical results with the experimental outcomes.


INTRODUCTION
he development of accurate and efficient finite element (FE) codes, based on enhanced beam formulations, is a significant challenge in many engineering fields and, in particular, in structural engineering.In today's professional structural applications it is often required to analyze large scale structures with irregular geometry, made from innovative composite materials, under severe loading conditions, especially in high seismicity areas.Hence, to accurately describe the global nonlinear structural response, as well as the local distribution of stresses and damaging paths, it is of great interest to formulate enhanced beam FEs, taking into account nonlinear geometric and constitutive behavior.To this end, it was widely demonstrated that the fiber approach is an efficient tool for introducing sufficiently general 3D nonlinear constitutive relationships, including plasticity, damage, crushing and so on [1,2].Moreover, the multi-axial coupling between the beam stress resultants is described.

T
The classical beam formulations, standing on Euler-Bernoulli or Timoshenko theory, assume that beam sections remain plane during the loading process.Although this results in very simple and effective formulations for describing the coupling between axial force and biaxial bending moment, an enriched kinematic description at the section level is required to include shear and torsion [3,4].It is worthwhile noting that the models proposed in literature, although including variable shear distributions and section warping, usually ignore the interaction between the beam element sections.Hence, these fail in describing the local response near the boundaries and the shear lag phenomenon, widely observed and discussed in literature [5,6], which significantly influence the normal stress distributions on the crosssections.The adoption of these formulations can lead to an incorrect underestimated evaluation of the stress values under warping constraints.In this work a 2-node 3D beam element, derived on the basis of the Hu-Washizu variational potential, is presented.A four-field mixed formulation is adopted, assuming as independent fields the standard displacements, the stresses, the strains and the warping displacements.Indeed, to take into account the warping of the cross-sections related to shear and torsion, an enriched kinematic description is introduced, as proposed in [7].The warping field is interpolated at two levels: along the element axis and at each section, where a variable number of warping degrees of freedom (DOFs) is defined.In literature, the adoption of interpolation functions derived, analytically or numerically, by solving the linear elastic problems for specific cross-section geometries under torsion was proposed.Herein, Lagrange polynomials are adopted both for the interpolation along the axis and for describing the warping profile on the section.These result more general to describe inelastic material response and suitable for performing the required numerical integration on the cross-section area.Aiming at reproducing the degrading mechanisms typical of many engineering materials, and in particular the brittle-like materials, a nonlinear stress-strain relationship with damage is introduced [8].Two distinct scalar damage variables are defined, reproducing the damaging mechanisms in tension and compression, and evolving according to laws defined in terms of an equivalent strain measure.The irreversibility of the two damaging processes is assumed.Furthermore, the unilateral effect is taken into account assuming the tensile damage greater than the compressive one, during the overall loading process.A fiber discretization of the cross-sections is adopted.Hence, the stress-strain constitutive relationship is defined at the material point level and integrated over the cross-section to derive the generalized stress measures.The FE formulation is defined in the basic system, i.e. in the reference system where element rigid body motions are removed, and under the assumption of small strains.The proposed beam element is implemented in the FE program FEAP [10], which is used to perform the numerical analyses.The beam FE is validated by performing a simple test on a cantilever beam in the linear elastic range, and two applications on plain concrete beams with damage experimentally analyzed, by comparing the numerically obtained results with the experimental outcomes.

BEAM FORMULATION
3D beam FE formulation, based on the Hu-Washizu mixed variational principle, including material nonlinear behavior, is introduced.Warping effects of the cross-sections are accounted for, by introducing independent DOFs.Small displacement and strain assumption holds.In Fig. 1 the FE reference configuration in the global reference system (O, X, Y, Z) is shown.The local coordinate system (I, x, y, z) is also defined (Fig. 2).The nodal displacement vector U is defined as: where U I / J , Φ I / J are the global translation and rotation vectors at node I/J (Fig. 1).Similarly, the nodal force vector results as: with P I / J , M I / J being the force and moment vectors at node I/J.By restraining the element nodal displacements, so as to remove all possible rigid body motions, the basic displacement vector D is introduced as: A Figure 1: FE global reference system: nodal displacement vector components.
where , x J U is the translation of node J parallel to the local axis x, z,I / J Φ and y,I / J Φ are the rotations at node I/J around z and y axis, respectively, and x,J Φ is the rotation at node J around x axis (Fig. 2).The corresponding basic force vector Q is defined as: where , x J P is the force parallel to the local axis x at node J, z,I / J M and y,I / J M are the moments at node I/J around z and y axis, respectively, and x, J M is the moment at node J around x axis.The element deformation D is derived from the global nodal displacement U, using the following transformation: where the kinematic matrix is introduced as:

Kinematic and static description
The evaluation of the element stiffness matrix K and force vector Q, associated to a given element deformation D, is formulated with respect to the basic coordinate system.According to the standard beam formulation, based on the assumption of rigid plane sections, the section displacement vector u(x) is defined as: x u x v x w x θ x θ x θ x  u (7) where u(x), v(x) and w(x) are the translation components of the beam axis, and  x (x),  y (x) and  z (x) are the rotations of the cross-section (Fig. 3).To describe the warping of the beam cross-sections, the classical assumption of rigid plane sections is partially removed, by assuming that these remain rigid in their plane, but can undergo out of plane deformations.Indeed, the displacement at the typical point P of the cross-section is expressed as the composition of the rigid part u r (x,y,z) and the displacement associated to the warping uw(x,y,z) (Fig. 3).As a consequence of the assumption of undeformable sections in their plane, the warping displacement has non-zero values only in the x direction, i.e.: Figure 3: Cross-section rigid displacement and cross-section warping displacement.
Hence, the displacement of the point is expressed as: By applying the compatibility operator, the deformation vector at the generic point of the cross-section is evaluated as: ( 1 0 ) d(x) being the generalized section deformation vector, resulting as: where  G (x) is the axial strain,  z (x) and  y (x) the flexural curvatures,  x (x) the torsional curvature,  y (x) and  z (x) the shear strains.The vector w(x,y,z) describes the deformation due to the warping displacement uw(x,y,z), i.e.: ( The stress components work conjugated with the deformation quantities in (x,y,z) are collected in the stress vector (x,y,z), defined as: (13) where  x is the normal stress along the beam axis direction, and  xy and  xz are the shear stresses in the cross-section plane parallel to the y and z axes, respectively.By applying the virtual work equivalence, the following definition of the generalized stress vector q(x) arises: with A being the cross-section area, N(x) the axial stress, M z (x) and M y (x) the flexural moments, M x (x) the torsional moment, Ty(x) and Tz(x) the generalized shear stresses.According to the equilibrated formulation [1,11,12], the stress vector q(x) is expressed in function of the basic element force vector Q, as: with b(x) being the equilibrium matrix and q p (x) the generalized section stresses due to the loads distributed along the element axis.Finally, the force field p w (x,y,z) is introduced.This is work conjugated with the warping displacement u w (x,y,z) and arises when in the cross-section the warping displacement is constrained.

Warping displacement interpolation
The warping displacement field uw(x,y,z) is interpolated according to the classical approach based on the use of shape functions.In particular, the interpolation is provided at two levels: along the element axis x and over the cross-section.
Along the element axis l w points are defined, whose distribution is arbitrary; on these points 1D Lagrange polynomials Nw,i(x) are defined.In this work, the distribution related to the Gauss-Lobatto numerical integration rule is adopted to easily impose warping restraints at the element ends.At each of the points located along the beam axis corresponds a cross-section on which the warping displacements are interpolated.In particular, assuming that the typical cross-section can be decomposed in a set of rectangular portions, in each of these a regular distribution of mw points is adopted (Fig. 4).Hence, 2D Lagrange polynomials M w,j (y,z) are defined, each of which corresponding to an additional warping displacement DOF for the element.The warping displacement field u w (x,y,z) can thus be written firstly considering the interpolation at the section level: , , , , where   , , , w i i u x y z is the warping displacement field over the section i, u w,i is the vector collecting all the warping DOFs u w,ij of the points located on the section i and M w (y,z) is a matrix containing all the shape functions M w,j (y,z) defined for the generic cross-section.Then, the interpolation along the beam axis is considered: The vector pw,i, collecting the warping forces pw,ij at the points located on the section i, can be defined, which is work conjugated with u w,i .
As the displacement field uw(x,y,z) does not contain the cross-section rigid body motions, the interpolation functions Mw,j(y,z) need to be modify to ensure the elimination of these rigid body motions.This elimination follows a specific procedure based on the definition of a projection matrix Pr and of another matrix k  , used to modify the warping stiffness matrix, as shown in Element variational formulation.

Element variational formulation
The equations governing the element state determination are derived on the basis of the Hu-Wahizu variational principle.
The following functional, depending on the four independent fields u(x), d(x),  (x,y,z) and uw(x,y,z), is introduced: where W(d, w ) is the internal potential energy and p contains the loads distributed along the element axis.Enforcing the stationarity of the functional  with respect to the four independent fields, the following governing equations are derived: Stationarity with respect to Governing equation Eq. ( 19) is the element equilibrium equation, rp P being the element nodal forces that arise from the distributed loads, Eq. (20) represents the constitutive law characterizing the response at the generic point of the cross-section, while Eq. ( 21) is the element compatibility equation, here enforced in weak form.Note that these are the classical equations derived from a standard three-field mixed FE formulation.Eq. ( 22), instead, represents the section equilibrium condition related to the warping effects, i.e., it requires that the warping loads p w,i at each cross-section are equal to the integral of the stresses x w q and yz w q due to the section warping, with: ) being matrices with dimensions 3  m w , composed as follows: To determine the solution of the nonlinear structural problems, the governing Eq. (19)(20)(21)(22) need to be linearized, resulting as: After some manipulations, the following set of equations is obtained: Eq. (30) represents the generalized section constitutive law and contains two contributions.The first term takes into account the standard deformation d associated to the cross-section rigid motion, depending on the section stiffness matrix, defined as: The second term is related to the section deformations associated to the warping.The section warping stiffness matrices are defined as follows: where: ( ) and identity matrix and R a 3 w m  matrix that represents the rigid body motions of the element section, i.e. a matrix containing the coordinates of the w m warping points as follows: Finally, V is a matrix containing the average value and the first moments of the shape functions over the cross section: According to the equilibrated beam formulation, the flexibility matrix F can be evaluated, by manipulating Eq. (29-32), as the sum of two symmetric matrices as: where S F is the classical flexibility matrix, while w F represents the contribution due to the warping effects.The matrices w K , ww K and B w are derived by substituting in the warping equilibrium Eq. (32) the expression of d obtained by Eq. (30).Hence, in compact form, Eq. (32) results as: where w U and w P are the increments of the two vectors collecting the warping displacements of all points at all sections used for the interpolation and the warping forces at the same points, respectively, which are defined as: he proposed FE is implemented in the numerical code FEAP [10], which is used to perform all the numerical analyses.Based on a step-by-step time discretization, the nonlinear problem in each time step is solved by adopting a classical iterative Newton-Raphson algorithm.At each iteration FEAP requires at the element level the computation of the element stiffness matrix +1 i K and of the force vector +1 i P , that are used to assembly the global stiffness matrix and the residual vector.In this section, the proposed algorithm for the element state determination is described.It is a generalization of the one presented in [7] and can be summarized as follows: 1.At the global Newton-Raphson current iteration ' +1' i , the nodal displacement vector +1 i

U
and its increment +1 i U with respect to the previous iteration ' ' i are given; 2. The basic element deformation increment +1 i D with respect to the previous iteration is computed, removing from +1 i U the element rigid body motions through the matrix T : The element basic forces are updated, using the element basic stiffness matrix at the previous iteration: To enforce the satisfaction of both the element equilibrium and compatibility conditions, a nested iterative procedure is performed (the superscript ' ' j and ' +1' j denote the previous and current internal iteration, respectively): a.By means of Eq. ( 35), the increment of the warping displacement vector at each section, where it is not constrained, is computed: and the increment of the warping force vector at the sections, where the corresponding displacement is constrained, is evaluated: At each section, the increment of the deformation vector is evaluated: At each point over the cross-section, the strains are computed and the resulting stress vector and material stiffness matrix are determined through the constitutive law: The stress vector +1 j


and the material stiffness matrix are integrated over the cross-section to obtain the generalized section stress vector and the section stiffness matrix:

B
are computed using Eq.(37); f.The section deformation vector is updated by adding to it the residual that arises from the difference between the balanced and the constitutive generalized section stress vectors: g.The compatible basic element deformation vector, including the above residual, and the basic element stiffness matrix are computed: h.The basic element force vector is computed as follows: The following residual internal energy is evaluated: If this results less then a specified tolerance or max j j  (with max j the maximum number of internal iterations), the nested iterative procedure is stopped and the superscript j is substituted with the superscript i , otherwise the procedure returns to step 4.a.It is worth noting that the non-iterative form of the proposed algorithm, i.e. without the nested iterative procedure of the step 4, can be obtained simply setting max 1 j  . 5. The section deformation vector is again updated: as well as the generalized stress vector, according to the equilibrium equation: The element stiffness matrix K and the force vector P are evaluated:

NUMERICAL APPLICATIONS
n this section some examples are reported to show the capability of the proposed FE in representing structural responses where the warping of the cross-sections and its interaction with the damaging mechanisms are relevant.The first one concerns the analysis of a channel-shaped cantilever beam subjected to a torsional load, under the hypotesis of linear elastic material; in this example, the influence of the uniform and non uniform warping distribution along the element axis is investigated.Other two applications are shown concerning the analysis of the damaging mechanisms due to the torsion in plain concrete beams.In all the applications, the Gauss-Lobatto integration rule is adopted to evaluate the integrals along the element axis together with a fiber model discretization of the cross-section, where the mid-point rule is used [9].

Twist of a channel-shaped cantilever beam
The first application considered was studied by Vlasov [16], Tralli [17], Capurso [18] and Back and Will [19].It is a cantilever beam, whose geometry is represented in Fig. 5, and an isotropic linear elastic constitutive law is adopted.The beam is modeled by a single FE with six warping points along the axis, 6 w l  , i.e. quintic shape functions, and with the distribution of the warping points over the cross-section shown in the same figure.The case of warping displacements restrained at the fixed end is considered and the rotations of the cross-sections of the beam are determined.The results are illustrated in Fig. 6.In the same figure the warping profile of the free end section is represented.evaluated on the basis of the thin-walled beams theory, is really more flexible than the other ones, as standard FEs are not able to take into account the contribution to the stiffness offered by the constraints applied to the warping displacements.The same solution can be obtained with the proposed FE, removing the warping constraints at the fixed end of the beam.As it is known, in thin-walled beams the effect of the warping restraints at the element ends does not vanish quickly along the beam axis, but it involves the whole element.The proposed FE is able to correctly describe this phenomenon, as it is shown in Fig. 7, where the distribution along the beam of the axial stress x  due to the shear-lag effect and of the axial displacement of the lower right corner of the section is depicted.Moreover, the x  distribution on the entire section is shown in Fig. 8 for different values of the abscissa x.  (x,y,z) and axial displacement u(x,y,z) at the lower right corner of the cross-section.

Damage model
To reproduce the degrading processes of the concrete material under several loading conditions, a constitutive damage model based on the introduction of a scalar damage variable D is adopted [8].The following constitutive law, relating the strain and stress tensors E and Σ , is defined at the fiber level: where E and  are the material Young's modulus and Poisson ratio, and I is the 6 6  identity matrix.
The damage variable is bounded in the range   0 ,1 , where D = 0 correponds to the initial undamaged state of the material and D = 1 to the complete degraded state. (x,y,z) over the cross-section in the channel-shaped cantilever beam with warping restraints at the fixed end.
For a given strain state of the material, the value assumed by D is computed by defining a damage associated variable, which governs its evolution.In particular, an equivalent strain eq,t  is defined on the basis of the positive part of the principal strains k  as: To model the unilateral effect, two different damage variables are considered: D t describing the damage related due to tensile deformation states, and D c reproducing the damage in presence of compressive deformation states.The following evolution laws are stated for both: where 0  is the initial damage threshold, and A i and B i are material parameters.The resulting damage variable D is a combination of the two variables D t and D c defined as: The weighting coefficients i  are defined as follows: with: where t,k  and c,k  are the principal strains evaluated on the basis of the two tensors t E and c E , defined as: -1 -1 and where i Σ is the stress tensor containing only the positive or the negative part of the principal stresses, and Λ is the secant stiffness matrix, depending on the damage D. The parameter  is a shear factor assumed equal to 1.06.
The adopted stress-strain law is presented in Fig. 9.a in the case of uniaxial tension and compression, by adopting the material parameters reported in Fig. 11.In Fig. 9.b the damage evolution law is presented.

Non local formulation
As known, in presence of a strain-softening constitutive law, localization problems and related mesh-dependency of the FE solution arise, thus requiring the adoption of a regularization technique.In the proposed FE, the localization of the generalized deformations at the integration points along the beam axis can appear, as typical in the force-based element [11].To overcome this problem, a nonlocal integral regularization technique is introduced.It is based on the nonlocal definition of the generalized section deformation d , computed in each cross-section as: where G x x  represents the distance between the quadrature point located at G x and all the others, and the parameter c L influences the extension of the averaging process.Then, at each point of the cross-section the total strains are evaluated as: and on the basis of these the damage parameters used for the constitutive law in the step 4.c are computed.

Plain concrete beams subjected to torsional loads
The applications presented in this section were introduced by Karayannis [14], based on experimental tests, and were also analyzed by Mazars [15].It is a beam, whose total length is 160 cm, divided into three parts: two end parts, properly reinforced, so as to remain elastic, and a plain concrete middle part, characterized by a cracking behavior (Fig. 10).A torque x M is applied at both the ends with a particular system that allows all the cross-sections of the beam to undergo warping displacements.Two shapes are considered for the beam cross-section: a rectangular one and a T-shaped one.For the rectangular beam, the adopted mesh is made from three FEs, one for each part of the specimen.Only one warping interpolation point is located along the axis, 1 w l  , while 4 4 16 w m    warping points are used over the crosssection.In fact, in this case the warping distribution along the element axis is uniform and has a cubic distribution in the principal directions over the section.Four Gauss-Lobatto integration points are used in each FE and a grid of 8 16  fibers is adopted.The material parameters are given in Fig. 11.The modulus E is not specified in [14], thus it is taken so as to reproduce the experimental initial stiffness.In Fig. 11, the global response curves and the warping displacement distribution are plotted, comparing the results obtained with the proposed FE with the experimental and analytical outcomes presented in [14].The curves show that the proposed FE is able to correctly predict the strain-softening behavior of the material, experimentally observed.Instead, the same accuracy can not be accomplished with standard FE formulations.In fact, the lower curve represents the response obtained under the hypothesis of rigid cross-sections, but correcting the section torsional inertia J with a torsion factor, that is J =  H B 3 with  = 0.229, determined by using the semi-analytical solution based on the Fourier series (according to the ratio H/B = 2.0).In this case, the peak load is significantly underestimated.This is due to the incorrect distribution of the damage variable D over the cross-section, that without warping is similar to that of a circular section, where no warping occurs [15].Comparing the damage distribution in the two cases (Fig. 12), i.e., that obtained with the proposed FE and with rigid section hypothesis, it is evident that, for the same value of rotation x  , in the element without warping a larger part of the section is characterized by high values of damage, thus the torsional stiffness of the beam is lower than that of the element with warping.For higher values of the rotation x  , the distribution of the damage in the FE with warping tends to show a pseudo-circular shape, similar to that of the element without warping.Thus, it seems that the two responses tend to coincide.However, the experimental tests in [14] show that the post-peak part of the response is usually characterized by abrupt failure of the specimen, so high values of the rotation are not reached.The case of the T-shaped section is then analyzed, using the same mesh and the same warping point distribution along the axis as in the case of the rectangular section.In this case, the fiber mesh grid and the warping points distribution shown in Fig. 13 are used over the cross-section, in order to have cubic warping interpolation functions in the principal directions over each rectangular part of the section.Four Gauss-Lobatto integration points are used in each FE.The material parameters are given in Fig. 14, where also in this case the modulus E is set so as to reproduce the experimental initial stiffness.In the same figure, the global response curves and the warping displacement profile are plotted, while the damage distributions are shown in Fig. 15.For this section, the response obtained under the hypothesis of rigid cross-sections is related to the value of the torsional inertia J = 4010 cm 4 , derived by using the semi-analytical solution based on the Fourier series.Similar considerations as in the case of the rectangular section can be made.

CONCLUSIONS
n efficient 3D beam FE formulation was presented, aiming at reproducing frame structural responses with multiaxial coupling between the stress resultants, including shear and torsion.To describe the response of structural elements, where warping mechanisms are relevant, the out of plane deformation of the beam cross-sections was A introduced, partially removing the standard rigid section assumption.The warping displacement field has been independently interpolated along the beam axis and on the cross-sections by using Lagrange polynomials, thus allowing to model cross-sections with arbitrary geometries.The proposed enhanced beam formulation overcomes some limitations of the existing models, thus allowing to describe the influence of the boundary conditions on the warping distribution, as well as the shear lag phenomenon.The material nonlinear behavior was taken into account, by adopting a fiber discretization of the cross-sections and considering a 3D constitutive relationship at the fiber level.Thanks to the adoption of a scalar damage model for concrete-like materials, the global response curves of experimental concrete beams under torsional loads have been correctly reproduced, as well as the damage distributions over the cross-sections.To this end, it was very interesting to underline the relevant influence of the warping on the beams damaging behavior.The implemented algorithm adopted for the element state determination has resulted robust, computationally efficient and accurate in reproducing well known reference examples, as well as more complex nonlinear structural responses.

Figure 4 :
Figure 4: Warping interpolation points and Lagrange shape functions along beam axis (left) and on the cross-section (right).

Figure 6 :
Figure 6: Torsional rotation x  of the cross-sections of the channel shaped cantilever beam and warping profile of the free end

Figure 10 :
Figure 10: Plain concrete rectangular beam subjected to end torsion load: specimen configuration and FE discretization.

Figure 11 :
Figure 11: Response of the rectangular beam subjected to concentrated torque: moment vs rotation curve and warping profile.

Figure 12 :
Figure 12: Rectangular beam: evolution of damage over the element cross-section without (first raw) and with (second raw) warping.

Figure 13 :
Figure 13: Plain concrete T-shaped beam subjected to end torsion load, warping points distribution and fiber discretization.

Figure 14 :
Figure 14: Response of the T-shaped beam subjected to concentrated torque: moment vs rotation curve and warping profile.

Figure 15 :
Figure 15: T-shaped beam: evolution of damage over the cross-section of the element without (first raw) and with (second raw) warping.