# Adaptive, Small-Rotation-Based, Corotational Technique for Analysis of 2D Nonlinear Elastic Frames

Authors: Jaroon Rungamornrat; Saethapoom Sihanartkatakul; Pattawee Kanchanakitcharoen
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper presents an efficient and accurate numerical technique for analysis of two-dimensional frames accounted for both geometric nonlinearity and nonlinear elastic material behavior. An adaptive remeshing scheme is utilized to optimally discretize a structure into a set of elements where the total displacement can be decomposed into the rigid body movement and one possessing small rotations. This, therefore, allows the force-deformation relationship for the latter part to be established based on small-rotation-based kinematics. Nonlinear elastic material model is integrated into such relation via the prescribed nonlinear moment-curvature relationship. The global force-displacement relation for each element can be derived subsequently using corotational formulations. A final system of nonlinear algebraic equations along with its associated gradient matrix for the whole structure is obtained by a standard assembly procedure and then solved numerically by Newton-Raphson algorithm. A selected set of results is then reported to demonstrate and discuss the computational performance including the accuracy and convergence of the proposed technique.


---

## Body

## 1. Introduction

It has been well recognized that various existing conventional techniques based either on linearized or simplified kinematics are generally incapable of modeling several crucial structural behaviors (e.g., buckling and postbuckling responses,P-delta effects, stability characteristics and bifurcations, etc.) or yield inaccurate results when they are applied to the modeling of structures accounted for highly geometric nonlinearities. For this reason, integration of a proper nonlinear kinematics to gain a more physically suitable mathematical model that can simulate practical situations to the level of complexity involved becomes essential. Analysis software packages fully equipped with sophisticated mathematical models have been emerged continuously with the primary objectives to enhance the capability of the prediction of structural responses under various circumstances, for example, presence of strong nonlinearities, various types of practical excitations, and intricate interactions with surrounding environments.Treatment of geometric nonlinearities in structural modeling has a very long history since the age of well-known mathematicians such as Euler, Lagrange, and Kirchoff (see extensive historical review in [1]). In that era, a problem of finding exact deformed shape of beams and columns, commonly known as an elastica problem, became very popular and mathematically challenging. However, due to lack of powerful, computationally aided tools and complexity of governing differential equations posed by highly nonlinear kinematics, solutions of most problems were derived based on analytical techniques and, therefore, limited only to simple structural configurations and loading conditions.As a result of the invention of computational devices and breakthrough in the field of numerical analysis, various powerful and robust numerical techniques have been continuously proposed, for past several decades, and employed to solve more practical and large-scale elastica problems. Some of those investigations relevant to the present study are briefly summarized below not only to indicate the historical background and the recent advances of computational techniques in this specific area but also to identify the gap of knowledge and the novel aspect of the current investigation. Lewis and Monasa [2] investigated the deformed shape of a cantilever beam which is made of nonlinear elastic materials of the Ludwik type (i.e., the uniaxial stress-strain relationship follows a pure power law) and subjected to a concentrated force by using the fourth-order Runge-Kutta method. Later, Monasa and Lewis [3] generalized their previous work to study a Ludwik-type cantilever beam under multiple concentrated forces. Prathep and Varadan [4] developed an iterative numerical scheme to investigate the inelastic postbuckling behavior of a cantilever column. In their work, the material nonlinearity was modeled by using Ramberg-Osgood stress-strain relation. Kounadis and Mallis [5] studied the influence of the slenderness ratio and a material parameter contained in the nonlinear elastic constitutive model on the buckling behavior of a bar. In 1989, a large deflection of a linearly elastic cantilever beam under a rotational, arbitrarily distributed load was studied by B. N. Rao and G. V. Rao [6] using the fourth-order Runge-Kutta technique. Wang [7] proposed a numerical procedure based on a perturbation technique to investigate the postbuckling behavior of a prismatic, roller-clamped, linearly elastic column under the applied end load. Later, Lee and Oh [8] used Runge-Kutta and Regula-Falsi methods to solve the elastica and buckling of a simply supported tapered column subjected to the end load. The postbuckling behavior of a linearly elastic cantilever column under the combined action of a uniformly distributed load and a concentrated force at the tip was also studied by Lee [9]. In the analysis, a technique based on Butcher’s fifth-order Runge-Kutta method was employed to obtain numerical results. Later, Lee [10] employed the shooting method and bisection algorithm to determine the large deflection of a nonlinear elastic cantilever beam under a uniformly distributed load and a concentrated force at its free end. Similar to the work of Lewis and Monasa [2] and Monasa and Lewis [3], the material model of Ludwik-type was assumed in the analysis. Dado et al. [11] investigated the computational performance of various techniques such as a semianalytical method, a numerical integration scheme, and a finite element method in the postbuckling analysis of a column consisting of two segments linked by a rotational spring. In 2005, Dado and Al-sadder [12] studied the large deflection of prismatic and nonprismatic cantilever beams under various types of applied loads. In their work, a technique based on the minimization of the residual error of the governing differential equation along with the representation of the beam rotation in terms of polynomial functions was employed. Vaz and Patel [13] applied both the shooting and Runge-Kutta methods to study the postbuckling behavior of a simply supported column. The material nonlinearity was handled by assuming a bilinear moment-curvature relationship. The influence of the rotational spring at the base of a cantilever column and the tip follower force on its deformed shape was also investigated by Shvartsman [14] by using a technique requiring no iterative procedure. Lacarbonara [15] investigated the postbuckling behavior of a nonprismatic nonlinearly elastic rod by using the higher-order perturbation technique. Later, Wang et al. [16] proposed an efficient and accurate technique, known as Homotopy analysis method (HAM), to solve a large deflection of a cantilever beam subjected to a concentrated force at the free end. Banerjee et al. [17] employed the shooting method along with the adomian decomposition to study a linear elastic cantilever beam under arbitrary loading conditions and containing an interior inflection point. The influence of discontinuity conditions and initial displacements on responses of frames undergoing large displacement and rotation was studied by Hu et al. [18]. In their analysis, a technique called the differential quadrature element method (DQEM) was utilized to construct numerical solutions. Lanc et al. [19] performed a large-displacement analysis of beam-type structures using the one-dimensional nonlinear finite element technique. In their work, the material nonlinearity was treated by using an elastic-plastic constitutive model. Later, Shatarat et al. [20] investigated a nonlinear elastic, rhombus frame subjected to a pair of opposite forces acting along its diagonal. In their study, the relationship between the displacement and applied forces at the corner was obtained by using a semianalytical technique. Chen [21] proposed a technique based on the moment integral treatment to determine the large deflection of a cantilever beam and results were found to be in good agreement with those obtained by the elliptic integral method and the commercial finite element package ANSYS©.Recently, Rungamornrat and Tangnovarad [22] proposed a semianalytical technique based on a direct stiffness procedure to analyze two-dimensional, linearly elastic, inextensible frames of arbitrary geometry and undergoing large displacement and rotation. In their formulation, the element tangent stiffness matrix and the residual load vector were derived in an exact form using the conventional elliptic integral technique. Later, Douanevanh et al. [23] generalized the work of Rungamornrat and Tangnovarad [22] to remove the inextensibility assumption. It is important to remark that while their technique yielded highly accurate results comparable to analytical solutions without the requirement of mesh refinement, it still requires solving a system of exact nonlinear algebraic equations in the calculation of the element tangent stiffness matrix and evaluation of the load residuals. In addition, that system of equations is strongly dependent on the presence of inflection points within elements. More recently, Sinsamutpadung et al. [24] developed a simplified numerical technique to solve the large displacement and rotation of two-dimensional frames. Their proposed strategy was based primarily on the displacement decomposition and a linearized kinematics for the displacement measured in the corotational coordinate system of an element. While their technique was found promising and robust, the formulation is still limited to frames made of linearly elastic materials.For various structures found in practices, as their displacement induced by excitations becomes sufficiently large (e.g., structures in the vicinity of collapses), it generally introduces large deformations or strains in certain regions and, as a result, linear constitutive models cannot be utilized to accurately represent the material behavior in this regime. In this paper, the technique proposed by Sinsamutpadung et al. [24] is enhanced by incorporating the material nonlinearity in the analysis of frames undergoing large displacement and rotation. The following sections present basic equations and assumptions, element force-displacement relations in different reference coordinate systems, load-displacement relationship for a structure, solution procedures, a selected set of numerical results, and conclusions and remarks.


## 2. Basic Equations for Small-Rotation-Based Kinematics

Consider a straight and prismatic member of lengthl, subjected to applied loads at both ends as shown schematically in Figure 1(a). The member is made of a homogeneous, hypothetical material with the moment-curvature and axial-force-strain relationships at any cross section prescribed a priori. An undeformed state of this member, represented by a line connecting the centroid of all cross sections and its subsequent deformed state resulting from the action of applied loads, is also indicated in Figure 1(a).(a) Schematic of undeformed and deformed members and (b) free-body diagram of an infinitesimal elementdx.
(a)
(b)By assuming that the member undergoes small displacement and rotation relative to its undeformed state, the static equilibrium equations can be formulated based on the known undeformed geometry (see Figure1(b) for a free-body diagram of the undeformed element dx) and it finally leads to
(1)dfdx=0;dVdx=0;dmdx=V(x),
where f, V, and m are the axial force, shear force, and bending moment at any cross section, respectively. Components of the displacement in the x- and y-directions, denoted by u=u(x) and v=v(x), are related to the axial strain ε and curvature κ via the small-rotation-based kinematics:
(2)ε=dudx;κ=d2vdx2.
It should be remarked that the relation (2) results directly from the approximations sinθ=θ and cos⁡θ=1 when the rotation θ is sufficiently small. The axial strain and curvature of any cross section are assumed to be related to the axial force f and the bending moment m via the fully uncoupled relations:
(3)f=AEε;κ=κ^(m),
where A denotes the cross sectional area of a member, E is Young’s modulus, and κ^=κ^(m) is any prescribed, sufficiently smooth function. It is worth pointing out that a form of the function κ^(m) generally depends on the stress-strain relationship and shape of the cross section. In the present study, the axial force is assumed to be related linearly to the axial strain and its influence on the moment-curvature relationship is neglected. It should be noted that the use of this simplified constitutive relation for the cross section does not aim to represent the real material behavior but is only for convenience purpose to demonstrate the computational performance of the proposed technique when dealing with nonlinear elasticity. The technique can be further extended to incorporate more complex material models by only changing the constitutive relation (3).


## 3. Element Force-Displacement Relationship

In this section, the force-displacement relationship and its best linear approximation for a generic member undergoing large displacement and rotation are developed. By following a well-known corotational technique, the total displacement of a member (measured from the undeformed stateS0 to its deformed state S) can be decomposed into a deformation-free displacement (measured from S0 to a state S^ which results directly from the rigid translation and rotation) and a displacement measured from S^ to S as illustrated in Figure 2. Three different reference coordinate systems, that is, a global system {0g;xg,yg}, a local system {0l;xl,yl}, and a corotational system {0*;x*,y*} are introduced for further reference. It should be remarked that {0*;x*,y*} is a moving coordinate system with its origin 0* always located at the end of a deformed member and the x*-axis coincident with a line connecting the two end points (see Figure 2).Figure 2
Schematic of a member in three different states, that is,S0,  S^, and S, and three different coordinate systems.A vector containing the total displacements and rotations at both ends of the member in the global coordinate system is defined byug=[u1gv1gθ1gu2gv2gθ2g]T, where the subscript {1,2} indicates the member end, the symbols u and v denote the displacement along the xg-and yg-axes, and θ denotes the end rotation. Similarly, a vector of the total end displacements and rotations of the member in the local coordinate system can be defined by ul=[u1lv1lθ1lu2lv2lθ2l]T. From a standard law of coordinate transformations, ug can readily be related to ul by
(4)ul=RT(γ)ug,
where R denotes a conventional transformation matrix which is a function of a member orientation angle γ. From the geometric consideration, a vector containing nonzero components of the end displacement and rotations of the member measured relative to the moving coordinate system {0*;x*,y*}, denoted by u*=[θ1*θ2*e*]T where θ1*, θ2*, and e* are clearly defined in Figure 2, can be expressed, in terms of ul and the angle β (i.e., an angle between the xl-axis and the x*-axis), as
(5)u*=u*(β,ul)=[θ1l-βθ2l-β(l+u2l-u1l)cos⁡β+(v2l-v1l)sinβ-l],
where l is the initial length of the member and the angle β can be obtained in terms of ul (i.e., β=β(ul)) from the following relation:
(6)(l+u2l-u1l)sinβ-(v2l-v1l)cos⁡β=0.
Vectors containing the end forces and end moments of each member in the global and local coordinate systems, denoted, respectively, by fg=[fx1gfy1gm1gfx2gfy2gm2g]T and fl=[fx1lfy1lm1lfx2lfy2lm2l]T, where all entries are ordered in a fashion consistent with ug and ul, can also be related, again, via the standard law of coordinate transformations:
(7)fg=R(γ)fl.
Similarly, a member force vector measured in the corotational coordinate system, denoted by f*=[fx1*fy1*m1*fx2*fy2*m2*]T, can be related to fl by
(8)fl=T(β)f*,
where T is a standard transformation matrix which is clearly dependent on the angle β.


### 3.1. Force-Displacement Relation in Corotational Coordinate System

To establish the relationship between the force vectorf*and the vector u*, it is assumed that the member length l resulting from the discretization is sufficiently small such that the rotation at any point in the current state S relative to the state S^ is infinitesimal. As a result, the governing equations (1) and (2) apply. It is worth noting that the movement of a member perceived by an observer in the corotational system {0*;x*,y*} is identical to that of a member restrained by a pinned support at the origin 0* and the roller support at the other end and subjected to applied end loads {m1*,m2*,fx2*} as shown in Figure 3.Figure 3
Schematic of a member restrained by a pinned support at the left end and a roller support at the right end and subjected to a set of end loads{m1*,m2*,fx2*}.By forming equilibrium equations for the entire member in its deformed state, all support reactive forces{fx1*,fy1*,fy2*} can readily be obtained in terms of the applied loads {m1*,m2*,fx2*} and the elongation e* by
(9)fx1*=-fx2*,fy1*=-fy2*=m1*+m2*l+e*.
By employing the first equation of (1), (2), and (3) along with the end boundary conditions, the applied force fx2* can be related linearly to e* by
(10)fx2*=AEe*l.
To establish the relation between the applied end moments {m1*,m2*} and the end rotations {θ1*,θ2*}, a technique based on the generalized moment area or curvature area equations is employed as briefly summarized below (also see details in the work of Pinyochotiwong et al. [25]). First, by using the equilibrium equations (1), the bending moment at any cross section can be expressed in terms of the applied moments {m1*,m2*} as
(11)m(x*;m1*,m2*)=m1*ψ1(x*)+m2*ψ2(x*),
where ψ1(x*) and ψ2(x*) are linear functions of x* defined by
(12)ψ1(x*)=x*l-1;ψ2(x)=x*l.
The curvature of any cross section can then be obtained by substituting (11) into the second equation of (3). Finally, the end rotations {θ1*,θ2*} are obtained from the following two curvature area equations:
(13)θ1*=∫0lψ1(x*)κ(m(x*;m1*,m2*))dx*,θ2*=∫0lψ2(x*)κ(m(x*;m1*,m2*))dx*.
Equations (13) implicitly define the end moments {m1*,m2*} in terms of the end rotations {θ1*,θ2*}; that is, m1*=m1*(θ1*,θ2*) and m2*=m2*(θ1*,θ2*). From Taylor series expansion of (13), the best linearized relationship between {θ1*,θ2*} and {m1*,m2*} is given by
(14){θ1*θ2*}={θ10*θ20*}+st*{m1*-m10*m2*-m20*},
where {m10*,m20*}T and {θ10*,θ20*}T are end moments and end rotations at a reference equilibrium state and st* is a flexural tangent flexibility matrix of a member with its entries explicitly defined by
(15)[st*]ij=∫0Lψi(x*)ψj(x*)∂κ∂mdx.
Relations (9), (10), and (13) implicitly define the force vector f* as a nonlinear function of the displacement vector u*; that is, f*=f*(u*). The best linear approximation of this nonlinear function is given by
(16)f*=f0*+kt*(u*-u0*),
where u0* and f0*=f*(u0*) are displacement and force vectors at a reference equilibrium state and kt*=∂f*/∂u* is the gradient matrix evaluated at u*=u0*. Entries of kt* can be computed as follows: ∂fx2/∂u* is obtained directly from (10), ∂{m1*,m2*}/∂{θ1*,θ2*} is clearly the inverse of the matrix st*, ∂{m1*,m2*}/∂e*=0, and ∂{fx1*,fy1*,fy2*}/∂u* can readily be obtained in terms of other computed entries via the relation (9).


### 3.2. Force-Displacement Relation in Local Coordinate System

By combining (5), (6), (8), and the implicit function f*=f*(u*) constructed in Section 3.1, it leads to the member force-displacement relationship in the local coordinate system
(17)fl=fl(ul)=T(β(ul))f*(u*(β(ul),ul)).
The best linear approximation of (17) can be obtained, again, by Taylor series expansion with the final result:
(18)fl=f0l+ktl(ul-u0l),
where u0l and f0l=fl(u0l) are local displacement and force vectors at a reference equilibrium state and ktl is a local element tangent stiffness matrix evaluated at ul=u0l. The form of ktl can be obtained from the chain rule for differentiations and the final expression can be written in a concise form as
(19)ktl=∂T∂βf*∂β∂ul+Tkt*∂u*∂ul,where the involved submatrices∂T/∂β,  ∂β/∂ul, and ∂u*/∂ul can readily be computed using the relations (4)-(5) and the fact that the transformation matrix T is given explicitly in terms of elementary functions of β (see explicit expression of these submatrices in the work of Sinsamutpadung et al. [24]). It can also be verified that ktl is essentially symmetric.


### 3.3. Force-Displacement Relation in Global Coordinate System

By combining (4), (7), and (17), it yields the member force-displacement relationship in the global coordinate system:
(20)fg=fg(ug)=R(γ)fl(RT(γ)ug).
Clearly, this relationship allows the end forces and end moments to be determined for any prescribed end displacements and rotations ug. Again, by performing Taylor series expansion about any admissible vector u0g, a linearized relationship between fg and ug in the neighborhood of u0g can be obtained as
(21)fg=f0g+ktg(ug-u0g),
where f0g=fg(u0g) and ktg is a global element tangent stiffness matrix evaluated at ug=u0g. From the chain rule for differentiations and the fact that γ is a constant, it can be verified that
(22)ktg=R(γ)ktlRT(γ).


## 3.1. Force-Displacement Relation in Corotational Coordinate System

To establish the relationship between the force vectorf*and the vector u*, it is assumed that the member length l resulting from the discretization is sufficiently small such that the rotation at any point in the current state S relative to the state S^ is infinitesimal. As a result, the governing equations (1) and (2) apply. It is worth noting that the movement of a member perceived by an observer in the corotational system {0*;x*,y*} is identical to that of a member restrained by a pinned support at the origin 0* and the roller support at the other end and subjected to applied end loads {m1*,m2*,fx2*} as shown in Figure 3.Figure 3
Schematic of a member restrained by a pinned support at the left end and a roller support at the right end and subjected to a set of end loads{m1*,m2*,fx2*}.By forming equilibrium equations for the entire member in its deformed state, all support reactive forces{fx1*,fy1*,fy2*} can readily be obtained in terms of the applied loads {m1*,m2*,fx2*} and the elongation e* by
(9)fx1*=-fx2*,fy1*=-fy2*=m1*+m2*l+e*.
By employing the first equation of (1), (2), and (3) along with the end boundary conditions, the applied force fx2* can be related linearly to e* by
(10)fx2*=AEe*l.
To establish the relation between the applied end moments {m1*,m2*} and the end rotations {θ1*,θ2*}, a technique based on the generalized moment area or curvature area equations is employed as briefly summarized below (also see details in the work of Pinyochotiwong et al. [25]). First, by using the equilibrium equations (1), the bending moment at any cross section can be expressed in terms of the applied moments {m1*,m2*} as
(11)m(x*;m1*,m2*)=m1*ψ1(x*)+m2*ψ2(x*),
where ψ1(x*) and ψ2(x*) are linear functions of x* defined by
(12)ψ1(x*)=x*l-1;ψ2(x)=x*l.
The curvature of any cross section can then be obtained by substituting (11) into the second equation of (3). Finally, the end rotations {θ1*,θ2*} are obtained from the following two curvature area equations:
(13)θ1*=∫0lψ1(x*)κ(m(x*;m1*,m2*))dx*,θ2*=∫0lψ2(x*)κ(m(x*;m1*,m2*))dx*.
Equations (13) implicitly define the end moments {m1*,m2*} in terms of the end rotations {θ1*,θ2*}; that is, m1*=m1*(θ1*,θ2*) and m2*=m2*(θ1*,θ2*). From Taylor series expansion of (13), the best linearized relationship between {θ1*,θ2*} and {m1*,m2*} is given by
(14){θ1*θ2*}={θ10*θ20*}+st*{m1*-m10*m2*-m20*},
where {m10*,m20*}T and {θ10*,θ20*}T are end moments and end rotations at a reference equilibrium state and st* is a flexural tangent flexibility matrix of a member with its entries explicitly defined by
(15)[st*]ij=∫0Lψi(x*)ψj(x*)∂κ∂mdx.
Relations (9), (10), and (13) implicitly define the force vector f* as a nonlinear function of the displacement vector u*; that is, f*=f*(u*). The best linear approximation of this nonlinear function is given by
(16)f*=f0*+kt*(u*-u0*),
where u0* and f0*=f*(u0*) are displacement and force vectors at a reference equilibrium state and kt*=∂f*/∂u* is the gradient matrix evaluated at u*=u0*. Entries of kt* can be computed as follows: ∂fx2/∂u* is obtained directly from (10), ∂{m1*,m2*}/∂{θ1*,θ2*} is clearly the inverse of the matrix st*, ∂{m1*,m2*}/∂e*=0, and ∂{fx1*,fy1*,fy2*}/∂u* can readily be obtained in terms of other computed entries via the relation (9).


## 3.2. Force-Displacement Relation in Local Coordinate System

By combining (5), (6), (8), and the implicit function f*=f*(u*) constructed in Section 3.1, it leads to the member force-displacement relationship in the local coordinate system
(17)fl=fl(ul)=T(β(ul))f*(u*(β(ul),ul)).
The best linear approximation of (17) can be obtained, again, by Taylor series expansion with the final result:
(18)fl=f0l+ktl(ul-u0l),
where u0l and f0l=fl(u0l) are local displacement and force vectors at a reference equilibrium state and ktl is a local element tangent stiffness matrix evaluated at ul=u0l. The form of ktl can be obtained from the chain rule for differentiations and the final expression can be written in a concise form as
(19)ktl=∂T∂βf*∂β∂ul+Tkt*∂u*∂ul,where the involved submatrices∂T/∂β,  ∂β/∂ul, and ∂u*/∂ul can readily be computed using the relations (4)-(5) and the fact that the transformation matrix T is given explicitly in terms of elementary functions of β (see explicit expression of these submatrices in the work of Sinsamutpadung et al. [24]). It can also be verified that ktl is essentially symmetric.


## 3.3. Force-Displacement Relation in Global Coordinate System

By combining (4), (7), and (17), it yields the member force-displacement relationship in the global coordinate system:
(20)fg=fg(ug)=R(γ)fl(RT(γ)ug).
Clearly, this relationship allows the end forces and end moments to be determined for any prescribed end displacements and rotations ug. Again, by performing Taylor series expansion about any admissible vector u0g, a linearized relationship between fg and ug in the neighborhood of u0g can be obtained as
(21)fg=f0g+ktg(ug-u0g),
where f0g=fg(u0g) and ktg is a global element tangent stiffness matrix evaluated at ug=u0g. From the chain rule for differentiations and the fact that γ is a constant, it can be verified that
(22)ktg=R(γ)ktlRT(γ).


## 4. Structure Load-Displacement Relation

To establish the structure load-displacement relationship, the given structure is first discretized inton straight and prismatic members. By defining U and F as a vector containing all nodal degrees of freedom (i.e., nodal displacements and rotations) and a vector of applied nodal loads, respectively, the force vector F can then be obtained by assembling all global member force vectors fg(1),fg(2),…,fg(n) and the displacement vector U can directly be related to all global member displacement vectors ug(1),ug(2),…,ug(n); that is,
(23)F=Assembly{fg(1),fg(2),…,fg(n)},(24)ug(i)=SelectionfromU,∀i∈{1,2,…,n}.
By using the fact that fg(i) is a function of ug(i) for all i∈{1,2,…,n} along with the relations (23)-(24), the force vector F is therefore a function of the displacement vector U; that is,
(25)F=F(U).
It should be remarked that, for any prescribed displacement vector U, the force vector F can be determined by applying the following equations consecutively: (24), (4), (5), the relation f*=f*(u*) developed in Section 3.1, (8), (7), and (23). The best linear approximation of the relation (25) in the neighborhood of U0 can also be obtained as
(26)F=F(U0)+Kt(U-U0),
where Kt=∂F/∂U is the structure tangent stiffness matrix evaluated at U=U0. By using the relations (23) and (24) and the fact that ktg is available for all members through the relations (21) and (22), Kt can therefore be obtained from a standard assembly procedure of all ktg.


## 5. Solution Procedure

It should be remarked that the nonlinearity of the relationship between the force vectorF and the displacement vector U stems directly from two main sources, the material nonlinearity present in the moment-curvature model (3) and the nonlinearity associated with the superposed large rigid body rotation β in the total displacement decomposition scheme.Although the procedure for determining the force vectorF for a given displacement vector U is clearly stated in the previous section, it still requires solving the force vector f* for a given displacement vector u* in the corotational coordinate system for all members. This nontrivial task can be achieved as follows: (i) the force fx2* is computed directly from (10); (ii) the end moments {m1*,m2*} are obtained for any prescribed {θ1*,θ2*} by solving nonlinear equations (13) using Newton-Raphson scheme; and (iii) the reactive forces {fx1*,fy1*,fy2*} are obtained from (9). It is noted that the linearized relation (14) along with (15) is useful for the iterative procedure carried out in step (ii) and the algorithm indicated above is restricted only to the constitutive relation (3).A more challenging and more practical task is associated with the determination of the deformed state of a structure under a specified set of applied loads. Since the relationshipF=F(U) cannot readily be inversed, the solution for the displacement vector U for a specified force vector F in each load step is obtained numerically via Newton-Raphson iteration. Details of the iterative procedure are given below.(a)
Choose an initial guess of the displacement vector,U(j) (j=0).(b)
Obtainug(i), ul(i), and u*(i) from (24), (4), and (5), respectively, for all i∈{1,2,…,n}.(c)
Obtainf*(i) for each u*(i) using procedure indicated above and then compute kt*(i) for all i∈{1,2,…,n} using procedure described in Section 3.1.(d)
Obtainfl(i) from (8) and ktl(i) from (19) for all i∈{1,2,…,n}.(e)
Obtainfg(i) from (7) and ktg(i) from (22) for all i∈{1,2,…,n}.(f)
ObtainF(j)and Kt(j) from direct assembly procedure.(g)
Check convergence condition||F-F(j)||/||F||<tol where tol is a specified tolerance. If such condition is satisfied, iteration is terminated and U(j) becomes an approximate solution; otherwise, proceed to the next step.(h)
Update the displacement vectorU(j+1) from U(j+1)=U(j)+Kt(j)(F-F(j)) and then return to step (b).It is remarked that, to accelerate the iterative procedure for the current load step, the information from the previous load step is utilized to provide the suitable initial guess. In addition, the nonuniform step size is employed to obtain results for the given loading history and this proves to enhance the computational efficiency of the proposed technique especially when the structure undergoes very large displacement and rotation, and the bending moment within members is far beyond the linear regime.


## 6. Adaptive Remeshing Algorithm

To ensure that the small-rotation-based kinematics employed in the formulation is suitable for modeling the relative displacement of any member in the corotational coordinate system{0*;x*,y*}, meshes used in the analysis must be sufficiently refined. In general, segments or regions whose curvature is significantly large must be discretized into several members whereas those possessing small curvature or being nearly straight require less number of members in the discretization. To achieve this crucial task efficiently and conveniently, a simple automatic remeshing scheme is integrated into the implemented numerical procedure as briefly described below.For the current load step, the analysis is first carried out using the mesh from the previous load step (the initial mesh is required and must be supplied for the first load step). Then, the maximum relative rotation of all members (i.e., maximum rotation measured in{0*;x*,y*} system) is calculated by using the curvature area equations and then is checked with the specified rotation limit. All members whose maximum relative rotation exceeds the given rotation limit are subdivided into two members by adding a new node at their center as shown in Figure 4. Once the new mesh is generated, the analysis is carried out again for this load step to obtain the new equilibrium state (following the procedure described in Section 5). The process of checking the requirement for member subdivision and performing the analysis of the current load step with the new mesh is repeated until the maximum relative rotation of all members in the current mesh is less than the rotation tolerance. It should be remarked that the information of all new nodes and members can be obtained by using results of old members before subdivision. This implemented adaptivity not only reduces effort for mesh generation but also provides a means to obtain a suitable mesh where the refinement is carried out only for regions needed.Figure 4
Schematic indicating a subdivision of a member into two members when its maximum relative rotation exceeds the rotation limit.


## 7. Numerical Results

To examine both the accuracy and convergence behavior of numerical solutions generated by the proposed technique, various structures under different loading conditions are fully investigated. A simple problem whose solution exists analytically is considered first not only to verify the formulation and implementation but also to demonstrate the relationship between solution error and the level of mesh refinement. The latter should provide the guideline for setting the rotation limit in the adaptive remeshing scheme to control the solution error. In addition, for problems whose analytical solution does not exist, benchmark results obtained from a reliable finite element package are utilized in the verification procedure. Note in particular that all the benchmark solutions used in the verification procedure are obtained based on the same constitutive relation (3) while the geometric nonlinearity is handled by exact kinematics for large displacement and rotation (i.e., the approximations sinθ=θ and cos⁡θ=1 are not removed). Once the implemented algorithm is tested, it is then applied to analyze more complex and practical structures and a selected set of results is reported and discussed. In the numerical study, the following nonlinear moment-curvature relationship is utilized:
(27)κκ0={MM0;MM0≤1,(1-α)+α(MM0)n;MM0>1,
where {α,n,κ0,M0} are model constants. It should be noted that M0 and κ0 denote the bending moment and curvature at the boundary of a linear regime and {α,n} are parameters used to describe the nonlinearity of the model. Since the relation (27) contains both linear and nonlinear regimes with the nonlinearity being controlled by {α,n}, it can be used to model a wide range of moment-curvature relationship of the beam cross section. When {α,n} are taken to be {1,1}, {α∈(1,∞),1}, and {α→∞,1}, the model recovers the linear, bilinear hardening, and bilinear without hardening M-κ relationship, respectively (see Figure 5(c)). Graphs of M-κ relationship of the family (27) for various values of {α,n} are shown in Figure 5; clearly, as α and n increase, the slope of the nonlinear part is descending.M
-
κ relationship for a family described by (27) for various parameters {α,n}.
(a)
(b)
(c)


### 7.1. Cantilever Beam Subjected to End Moment

Consider a cantilever beam of lengthl and cross-sectional area A. The beam is made of a homogeneous material with EAl/M0=1000, κ0l=0.3, α=1.25, and n=2 and it is subjected to a concentrated, counterclockwise moment M at the tip as shown in Figure 6(a). First, the analysis is carried out without the automatic remeshing algorithm and five meshes containing 1, 2, 4, 8, and 16 members shown in Figure 6(b) are utilized.(a) Schematic of a cantilever beam subjected to a concentrated momentM at the tip and (b) five meshes adopted in the analysis.
(a)
(b)To demonstrate the capability of the current technique to analyze structures undergoing large displacement and rotation, values of the applied momentM are chosen to be sufficiently large in the numerical experiments. The computed vertical and horizontal displacements at the tip of the beam from all five meshes are normalized by the length l and then reported in Table 1 for M/M0=3 along with the maximum relative rotation (measured in the corotational coordinate system). The corresponding deformed shapes are shown in Figure 7. It is evident that the numerical results converge to the exact solution as the mesh is refined. This is due to the fact that as the number of members in the mesh increases, the rotation in each member measured in the corotational coordinate system becomes smaller and this, as a result, improves the suitability of using small-rotation-based governing equations (1) and (2). Figure 8 shows quantitatively the relationship between the percent error of both the horizontal and vertical displacements at the tip of the beam and the maximum rotation of all members measured in the corotational coordinate system as the mesh is refined. Since the maximum rotation is a one-to-one correspondence to the number or length of discretized members, results indicated in Figure 8 can be used to obtain the rotation limit to achieve a required level of accuracy of an approximate solution.Table 1
Comparison of normalized computed horizontal and vertical displacements at the tip of cantilever beam and available analytical solution.
Mesh
u
/
l
% error
v
/
l
% error
Maximum relative rotation (deg)
1
−1.0791
2.989
0.9969
65.519
94.538
2
−1.0537
0.562
0.6764
12.314
47.269
3
−1.0492
0.132
0.6197
2.893
23.635
4
−1.0481
0.032
0.6066
0.713
11.817
5
−1.0479
0.008
0.6033
0.177
5.909Figure 7
Deformed shape of a cantilever beam subjected to tip momentM with M/M0=3 obtained from five different meshes.Figure 8
Relationship between the maximum rotation measured in corotational coordinate system and the percent error of horizontal displacement (u) and vertical displacement (v) at the tip.To ensure that the maximum rotation measured in the corotational coordinate system of each member is sufficiently small (i.e., not exceeding a specified rotation limit), the automatic remeshing scheme described in Section6 is integrated in the analysis procedure. Comparison of results with the exact solution for the same applied load level (i.e., M/M0=3), when such adaptive strategy is utilized, is shown in Table 2 for different values of the rotation limit. According to this set of results, it can be seen that as the rotation limit decreases, the number of members, together with the level of accuracy, increases. Deformed shapes of the beam for M/M0∈{2,3,4} and different rotation limits are also reported in Figure 9. The number of members in the final mesh for the given rotation limit and load level is also provided. It is apparent that numerical results, for this particular case, are almost indistinguishable from the exact solution although the large rotation limit (i.e., 20 degree) is used in the adaptive scheme. Note in addition to the fact that the number of members in the final mesh resulting from the adaptive scheme also depends on the curvature of the member; for this particular problem, as the applied end moment M increases, the curvature of the beam increases and, as a result, the mesh must be more refined to ensure that the relative rotation does not exceed the rotation limit.Table 2
Error of normalized vertical and horizontal tip displacement of cantilever beam using adaptive re-meshing scheme with various rotation limits.
Rotation limit (degree)
Number of members
% error ofu/l
% error ofv/l
20
8
0.0324
0.7125
10
16
0.0080
0.1774
5
32
0.0020
0.0444
2
64
0.0005
0.0110Figure 9
Schematic of deformed shapes of cantilever beam using adaptive remeshing scheme forM/M0∈{2,3,4}. The number of members in the final mesh for M/M0=2,3,4 is given by {32,16,8,4}, {64,32,16,8}, and {128,64,32,16} for rotation limit of 2, 5, 10, and 20 degrees, respectively.


### 7.2. Diamond-Shaped Frame Subjected to Vertical Load

Consider a diamond-shaped frame with a fixed support at the bottom corner and subjected to a vertical loadP at the top corner as shown schematically in Figure 10(a). All four sides of the frame have the same length 2l and the same properties, that is, with EAl/M0=1000, κ0l=0.2, α=1.25, and n=2. Note in particular that the moment-curvature relationship chosen in the numerical study is bilinear. In the analysis, three different load levels associated with Pl/M0∈{10,20,30} are considered, the initial mesh consisting of four members as shown in Figure 10(b) is utilized, and the rotation limit is set equal to 2 degrees to ensure the high level of accuracy.(a) Schematic of diamond-shaped frame subjected to vertical loadP at its top corner and (b) initial mesh used in the analysis.
(a)
(b)The deformed shapes obtained from the current technique are reported in Figure11 along with the benchmark solution generated by a reliable finite element analysis package. The final mesh for each load level resulting from the adaptive scheme is shown in Figure 12. It is evident that as the load increases, the mesh is more refined. In particular, most of members are located in the corner region since the corresponding curvature is relatively large. In addition, results shown in Figure 11 indicate that numerical solutions obtained from the proposed technique are highly accurate and almost identical to the benchmark solutions.Figure 11
Deformed shapes of diamond-shaped frame subjected to vertical loadP at the top corner.Final mesh resulting from automatic remeshing algorithm for (a)Pl/M0=10 (48 elements), (b) Pl/M0=20 (88 elements), and (c) Pl/M0=30 (152 elements).
(a)
(b)
(c)Note also that the deformed shape of the frame can still be accurately captured if the rotation limit changes from 2 to 20 degrees as shown in Figure13. In particular, the percent difference of the vertical displacement at the top corner, the horizontal displacement at point A (indicated in Figure 10), and the vertical displacement at point A obtained by using two different rotational limits (i.e., 2 and 20 degrees) are 0.323%, 1.235%, and 0.322%, respectively.Figure 13
Deformed shapes of diamond-shaped frame subjected to vertical load at the top corner. Results are reported forPl/M0=30 and two different rotation limits.


### 7.3. Cantilever Beam Subjected to End Loads

To demonstrate the influence of large displacement and rotation on values of internal forces such as the shear force, axial force, and bending moment at any cross section, let us consider a cantilever beam of lengthl and cross-sectional area A, made of a material with EAl/M0=10000, κ0l=0.3, α=1.25, and n=1. The beam is subjected to a horizontal force P, a vertical force Q, and a moment M at the tip as shown in Figure 14(a). In the analysis, an automatic remeshing scheme is used with the rotation limit equal to 2 degrees and an initial mesh containing only one member (see Figure 14(b)).(a) Schematic of a cantilever beam subjected to forces and moment at the tip and (b) initial mesh adopted in the analysis.
(a)
(b)In the numerical study, the end loads{P,Q,M} are applied proportionally such that 2Pl/M0=Ql/M0=2M/M0=λ, where λ is a load parameter. The deformed shapes obtained from the analysis for λ=10,20, and 50 are reported in Figure 15. It is remarked that three different values of the loading parameter λ are chosen to represent the beam undergoing small, intermediate, and large deflections, respectively. The axial force, the shear force, and the bending moment at any cross section of the beam for above proportional loads are obtained and reported along with those generated by the linear elastic analysis in Figures 16, 17, and 18. It is evident that results obtained from the proposed technique show significant deviation from the linear solution and such difference increases as the magnitude of applied loads becomes larger. This is mainly due to that the deformed configuration of the member at larger loads is significantly different from the undeformed configuration and this therefore boosts the influence of the axial-shear-bending interaction. In particular, the axial force and shear force obtained from the linear analysis are constant throughout the member (due to that the equilibrium is enforced in the undeformed state) whereas those obtained from the large displacement and rotation analysis vary nonlinearly along the member due to the significant change of the member axis. This similar trend was also observed in the case of the bending moment. However, the deviation from the linear elastic solution for the bending moment is less significant since the influence of the P-delta effect is reduced by the considerable reduction of the moment arm of the force Q when the beam undergoes very large displacement and rotation. Clearly, the increase in moment due to the P-delta effect is less than the decrease in moment due to the reduction of the moment arm for λ=50 and vice versa for λ=10,20. It can be noted also that as λ increases, the bending stiffness of the beam is further reduced since the region where the bending moment is larger than M0 almost occupies the entire beam.Figure 15
Deformed shapes of cantilever beam under various end forces and moments.Figure 16
Normalized axial force diagram of cantilever beam under end forces and moments.Figure 17
Normalized shear force diagram of cantilever beam under end forces and moments.Figure 18
Normalized bending moment diagram of cantilever beam under end forces and moments.


### 7.4. Portal Frame Subjected to Horizontal and Vertical Forces

As a final example, let us consider a problem of finding the deformed shape of a portal frame clamped at its bases and subjected to horizontal and vertical loads as shown in Figure19(a). The given frame is made of a material with Young’s modulus E and moment-curvature relationship described by (27) and all members have the same cross-sectional area A. An initial mesh containing only three members, indicated in Figure 19(b), is adopted in the analysis and the rotation limit equal to 2 degrees is utilized in the automatic remeshing scheme. To demonstrate the influence of material nonlinearity on the deformed shape of the frame, various values of material parameters {α,n} are considered along with EAl/M0=10000 and κ0l=0.2.(a) Schematic of a portal frame clamped at its bases and subjected to horizontal and vertical forces and (b) initial mesh adopted in the analysis.
(a)
(b)The normalized computed deformed shapes of the frame are reported in Figures20 and 21 for Pl/M0=20, Ql/M0=2, and various material parameters {α,n}. The load level considered in the analysis is chosen to ensure that the total displacement and rotation of the frame are sufficiently large to demonstrate the capability of the proposed technique to capture highly geometric and material nonlinearities. It is evident from this set of results that the frame becomes softer as α increases for a fixed n and, similarly, as n increases for a fixed α. This is clearly consistent with the moment-curvature relationship; that is, the larger the value of n and α, the higher the material nonlinearity and lesser material stiffness.Figure 20
Deformed shape of portal frame clamped at its bases and subjected to horizontal and vertical loads. Results are reported forn=2 and various values of α.Figure 21
Deformed shape of portal frame clamped at its bases and subjected to horizontal and vertical loads. Results are reported forα=2 and various values of n.


## 7.1. Cantilever Beam Subjected to End Moment

Consider a cantilever beam of lengthl and cross-sectional area A. The beam is made of a homogeneous material with EAl/M0=1000, κ0l=0.3, α=1.25, and n=2 and it is subjected to a concentrated, counterclockwise moment M at the tip as shown in Figure 6(a). First, the analysis is carried out without the automatic remeshing algorithm and five meshes containing 1, 2, 4, 8, and 16 members shown in Figure 6(b) are utilized.(a) Schematic of a cantilever beam subjected to a concentrated momentM at the tip and (b) five meshes adopted in the analysis.
(a)
(b)To demonstrate the capability of the current technique to analyze structures undergoing large displacement and rotation, values of the applied momentM are chosen to be sufficiently large in the numerical experiments. The computed vertical and horizontal displacements at the tip of the beam from all five meshes are normalized by the length l and then reported in Table 1 for M/M0=3 along with the maximum relative rotation (measured in the corotational coordinate system). The corresponding deformed shapes are shown in Figure 7. It is evident that the numerical results converge to the exact solution as the mesh is refined. This is due to the fact that as the number of members in the mesh increases, the rotation in each member measured in the corotational coordinate system becomes smaller and this, as a result, improves the suitability of using small-rotation-based governing equations (1) and (2). Figure 8 shows quantitatively the relationship between the percent error of both the horizontal and vertical displacements at the tip of the beam and the maximum rotation of all members measured in the corotational coordinate system as the mesh is refined. Since the maximum rotation is a one-to-one correspondence to the number or length of discretized members, results indicated in Figure 8 can be used to obtain the rotation limit to achieve a required level of accuracy of an approximate solution.Table 1
Comparison of normalized computed horizontal and vertical displacements at the tip of cantilever beam and available analytical solution.
Mesh
u
/
l
% error
v
/
l
% error
Maximum relative rotation (deg)
1
−1.0791
2.989
0.9969
65.519
94.538
2
−1.0537
0.562
0.6764
12.314
47.269
3
−1.0492
0.132
0.6197
2.893
23.635
4
−1.0481
0.032
0.6066
0.713
11.817
5
−1.0479
0.008
0.6033
0.177
5.909Figure 7
Deformed shape of a cantilever beam subjected to tip momentM with M/M0=3 obtained from five different meshes.Figure 8
Relationship between the maximum rotation measured in corotational coordinate system and the percent error of horizontal displacement (u) and vertical displacement (v) at the tip.To ensure that the maximum rotation measured in the corotational coordinate system of each member is sufficiently small (i.e., not exceeding a specified rotation limit), the automatic remeshing scheme described in Section6 is integrated in the analysis procedure. Comparison of results with the exact solution for the same applied load level (i.e., M/M0=3), when such adaptive strategy is utilized, is shown in Table 2 for different values of the rotation limit. According to this set of results, it can be seen that as the rotation limit decreases, the number of members, together with the level of accuracy, increases. Deformed shapes of the beam for M/M0∈{2,3,4} and different rotation limits are also reported in Figure 9. The number of members in the final mesh for the given rotation limit and load level is also provided. It is apparent that numerical results, for this particular case, are almost indistinguishable from the exact solution although the large rotation limit (i.e., 20 degree) is used in the adaptive scheme. Note in addition to the fact that the number of members in the final mesh resulting from the adaptive scheme also depends on the curvature of the member; for this particular problem, as the applied end moment M increases, the curvature of the beam increases and, as a result, the mesh must be more refined to ensure that the relative rotation does not exceed the rotation limit.Table 2
Error of normalized vertical and horizontal tip displacement of cantilever beam using adaptive re-meshing scheme with various rotation limits.
Rotation limit (degree)
Number of members
% error ofu/l
% error ofv/l
20
8
0.0324
0.7125
10
16
0.0080
0.1774
5
32
0.0020
0.0444
2
64
0.0005
0.0110Figure 9
Schematic of deformed shapes of cantilever beam using adaptive remeshing scheme forM/M0∈{2,3,4}. The number of members in the final mesh for M/M0=2,3,4 is given by {32,16,8,4}, {64,32,16,8}, and {128,64,32,16} for rotation limit of 2, 5, 10, and 20 degrees, respectively.


## 7.2. Diamond-Shaped Frame Subjected to Vertical Load

Consider a diamond-shaped frame with a fixed support at the bottom corner and subjected to a vertical loadP at the top corner as shown schematically in Figure 10(a). All four sides of the frame have the same length 2l and the same properties, that is, with EAl/M0=1000, κ0l=0.2, α=1.25, and n=2. Note in particular that the moment-curvature relationship chosen in the numerical study is bilinear. In the analysis, three different load levels associated with Pl/M0∈{10,20,30} are considered, the initial mesh consisting of four members as shown in Figure 10(b) is utilized, and the rotation limit is set equal to 2 degrees to ensure the high level of accuracy.(a) Schematic of diamond-shaped frame subjected to vertical loadP at its top corner and (b) initial mesh used in the analysis.
(a)
(b)The deformed shapes obtained from the current technique are reported in Figure11 along with the benchmark solution generated by a reliable finite element analysis package. The final mesh for each load level resulting from the adaptive scheme is shown in Figure 12. It is evident that as the load increases, the mesh is more refined. In particular, most of members are located in the corner region since the corresponding curvature is relatively large. In addition, results shown in Figure 11 indicate that numerical solutions obtained from the proposed technique are highly accurate and almost identical to the benchmark solutions.Figure 11
Deformed shapes of diamond-shaped frame subjected to vertical loadP at the top corner.Final mesh resulting from automatic remeshing algorithm for (a)Pl/M0=10 (48 elements), (b) Pl/M0=20 (88 elements), and (c) Pl/M0=30 (152 elements).
(a)
(b)
(c)Note also that the deformed shape of the frame can still be accurately captured if the rotation limit changes from 2 to 20 degrees as shown in Figure13. In particular, the percent difference of the vertical displacement at the top corner, the horizontal displacement at point A (indicated in Figure 10), and the vertical displacement at point A obtained by using two different rotational limits (i.e., 2 and 20 degrees) are 0.323%, 1.235%, and 0.322%, respectively.Figure 13
Deformed shapes of diamond-shaped frame subjected to vertical load at the top corner. Results are reported forPl/M0=30 and two different rotation limits.


## 7.3. Cantilever Beam Subjected to End Loads

To demonstrate the influence of large displacement and rotation on values of internal forces such as the shear force, axial force, and bending moment at any cross section, let us consider a cantilever beam of lengthl and cross-sectional area A, made of a material with EAl/M0=10000, κ0l=0.3, α=1.25, and n=1. The beam is subjected to a horizontal force P, a vertical force Q, and a moment M at the tip as shown in Figure 14(a). In the analysis, an automatic remeshing scheme is used with the rotation limit equal to 2 degrees and an initial mesh containing only one member (see Figure 14(b)).(a) Schematic of a cantilever beam subjected to forces and moment at the tip and (b) initial mesh adopted in the analysis.
(a)
(b)In the numerical study, the end loads{P,Q,M} are applied proportionally such that 2Pl/M0=Ql/M0=2M/M0=λ, where λ is a load parameter. The deformed shapes obtained from the analysis for λ=10,20, and 50 are reported in Figure 15. It is remarked that three different values of the loading parameter λ are chosen to represent the beam undergoing small, intermediate, and large deflections, respectively. The axial force, the shear force, and the bending moment at any cross section of the beam for above proportional loads are obtained and reported along with those generated by the linear elastic analysis in Figures 16, 17, and 18. It is evident that results obtained from the proposed technique show significant deviation from the linear solution and such difference increases as the magnitude of applied loads becomes larger. This is mainly due to that the deformed configuration of the member at larger loads is significantly different from the undeformed configuration and this therefore boosts the influence of the axial-shear-bending interaction. In particular, the axial force and shear force obtained from the linear analysis are constant throughout the member (due to that the equilibrium is enforced in the undeformed state) whereas those obtained from the large displacement and rotation analysis vary nonlinearly along the member due to the significant change of the member axis. This similar trend was also observed in the case of the bending moment. However, the deviation from the linear elastic solution for the bending moment is less significant since the influence of the P-delta effect is reduced by the considerable reduction of the moment arm of the force Q when the beam undergoes very large displacement and rotation. Clearly, the increase in moment due to the P-delta effect is less than the decrease in moment due to the reduction of the moment arm for λ=50 and vice versa for λ=10,20. It can be noted also that as λ increases, the bending stiffness of the beam is further reduced since the region where the bending moment is larger than M0 almost occupies the entire beam.Figure 15
Deformed shapes of cantilever beam under various end forces and moments.Figure 16
Normalized axial force diagram of cantilever beam under end forces and moments.Figure 17
Normalized shear force diagram of cantilever beam under end forces and moments.Figure 18
Normalized bending moment diagram of cantilever beam under end forces and moments.


## 7.4. Portal Frame Subjected to Horizontal and Vertical Forces

As a final example, let us consider a problem of finding the deformed shape of a portal frame clamped at its bases and subjected to horizontal and vertical loads as shown in Figure19(a). The given frame is made of a material with Young’s modulus E and moment-curvature relationship described by (27) and all members have the same cross-sectional area A. An initial mesh containing only three members, indicated in Figure 19(b), is adopted in the analysis and the rotation limit equal to 2 degrees is utilized in the automatic remeshing scheme. To demonstrate the influence of material nonlinearity on the deformed shape of the frame, various values of material parameters {α,n} are considered along with EAl/M0=10000 and κ0l=0.2.(a) Schematic of a portal frame clamped at its bases and subjected to horizontal and vertical forces and (b) initial mesh adopted in the analysis.
(a)
(b)The normalized computed deformed shapes of the frame are reported in Figures20 and 21 for Pl/M0=20, Ql/M0=2, and various material parameters {α,n}. The load level considered in the analysis is chosen to ensure that the total displacement and rotation of the frame are sufficiently large to demonstrate the capability of the proposed technique to capture highly geometric and material nonlinearities. It is evident from this set of results that the frame becomes softer as α increases for a fixed n and, similarly, as n increases for a fixed α. This is clearly consistent with the moment-curvature relationship; that is, the larger the value of n and α, the higher the material nonlinearity and lesser material stiffness.Figure 20
Deformed shape of portal frame clamped at its bases and subjected to horizontal and vertical loads. Results are reported forn=2 and various values of α.Figure 21
Deformed shape of portal frame clamped at its bases and subjected to horizontal and vertical loads. Results are reported forα=2 and various values of n.


## 8. Conclusions and Remarks

The numerical technique which is capable of modeling two-dimensional, nonlinearly elastic, extensible frames undergoing large displacement and rotation has been developed. The force-deformation relationship for all members in the corotational coordinate system has been established based on the small-rotation-based kinematics. The large displacement and rotation feature has been taken into account in the analysis by superposing the displacement measured relative to the corotational coordinate system by the large rigid body displacement. The material nonlinearity has been integrated in the modeling via the prescribed arbitrary nonlinear moment-curvature relationship. The final system of nonlinear algebraic equations has been solved by Newton-Raphson iterative scheme. To ensure the validity of using the small-rotation-based governing equations in the corotational system, an adaptive remeshing algorithm has been implemented to automatically obtain an optimal mesh for a specified rotation limit. The benefit of this proposed technique is to avoid the direct treatment of highly nonlinear differential equations in the development of force-displacement relationship and the corresponding element tangent stiffness matrix.From extensive numerical experiments, it has been found that the proposed technique yields very accurate numerical results compared with the benchmark solutions provided that the mesh used in the analysis is sufficiently refined. The adaptive remeshing scheme is proved to provide a convenient means to obtain an optimal mesh used to acquire accurate results; only meshes of regions or segments whose curvature is sufficiently large are refined. It is worth noting that the present study is still restricted to a simple model for the treatment of material nonlinearity. The moment-curvature relationship utilized is free of the axial interaction and applies only to history-independent materials. The former can potentially lead to an inaccurate prediction when applied to axially dominated structures whereas the latter renders the technique incapable of modeling inelastic structural responses. One possible extension of the current work to better model material nonlinearity is to incorporate a fiber model to directly construct the cross-sectional behavior from a basic constitutive model in the level of material points. This additionally eases the treatment of the axial-bending coupling behavior.



---

*Source: 640983-2014-04-03.xml*