# Analysis of Linearly Elastic Inextensible Frames Undergoing Large Displacement and Rotation

Authors: Jaroon Rungamornrat; Peerasak Tangnovarad
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

This paper presents an efficient semi-analytical technique for modeling two-dimensional, linearly elastic, inextensible frames undergoing large displacement and rotation. A system of ordinary differential equations governing an element is first converted into a system of nonlinear algebraic equations via appropriate enforcement of boundary conditions. Taylor's series expansion is then employed along with the co-rotational approach to derive the best linear approximation of such system and the corresponding exact element tangent stiffness matrix. A standard assembly procedure is applied, next, to obtain the best linear approximation of governing nonlinear equations for the structure. This final system is exploited in the solution search by Newton-Ralphson iteration. Key features of the proposed technique include that (i) exact load residuals are evaluated from governing nonlinear algebraic equations, (ii) an exact form of the tangent stiffness matrix is utilized, and (iii) all elements are treated in a systematic way via direct stiffness strategy. The first two features enhance the performance of the technique to yield results comparable to analytical solutions and independent of mesh refinement whereas the last feature allows structures of general geometries and loading conditions be modeled in a straightforward fashion. The implemented algorithm is tested for various structures not only to verify its underlying formulation but also to demonstrate its capability and robustness.


---

## Body

## 1. Introduction

It is well known that a small-deformation analysis of flexure-dominating structures (e.g., beams and frames) based primarily on linearized kinematics and fully decoupled axial-bending interactions (e.g., [1, 2]) can lead to results that are of insufficient accuracy, especially when the displacement and rotation of a structure are relatively large and the axial-bending coupling becomes significant (e.g., [3]). In addition, such so-called linear analysis provides limited information about either the stability of the structure (e.g., bifurcation loads and identification of the stability status of structures) or the behavior beyond a point of bifurcation (i.e., postbuckling behavior). Besides mathematical curiosity and computational challenge, the necessity to incorporate proper nonlinear kinematics in the mathematical model is obligatory and arises naturally in numerous practical applications, for example, modeling of beam-columns where the axial-bending interaction is crucial, analysis and design of flexible components of machines and systems vulnerable to postbuckling, collapse analysis of structures under severe loading conditions, and modeling of very slender and flexible structures where the displacement and rotation are substantial under service conditions.One simple approach that has been widely used to model geometric nonlinearity of structures is known as the second-order analysis (e.g., [3–7]). The key improvement of this approach from the linear analysis stems from the use of simplified nonlinear kinematics along with forming equilibrium equations based on a deformed state. The integration of this level of geometric nonlinearity enables the mathematical model to explore additional structural responses such as critical or bifurcation loads (e.g., [3, 4, 6]) and the interaction between the bending and axial effects (e.g., [3, 5, 7]). Nevertheless, the second-order analysis still possesses limited capabilities due to the use of low-order approximate kinematics. For instance, it still provides limited information on behavior of the structure beyond points of bifurcation (i.e., postbuckling behavior) and provides results of insufficient accuracy when the structure undergoes very large displacement and rotation relative to its dimensions. As a result, modeling of geometric nonlinearity based on exact kinematics has become an attractive alternative to circumvent all those limitations.A more sophisticated mathematical model incorporating exact kinematics (i.e., exact relationship between the displacement, rotation, and curvature) was introduced more than two centuries ago by Euler (1774) and Lagrange (1770–1773) in their study of finding an exact, elastic, or deformed curve of beams, known as an elastica problem (see also [8] for an extensive historical discussion). Later, Kirchhoff [9] made a significant progress by establishing an analogy between a problem of finding elastica of a cantilever column and a problem associated with the oscillation of a pendulum. With such simple analogy, a closed-form solution could be constructed using a so-called, elliptic integral method. Due to complexity posed by the exact kinematics, solutions of elastica problems in its toddler age, based purely on analytical techniques, were very limited to structures of simple geometries and loading conditions.Due to the rapid growth of powerful computing devices and robust numerical techniques, the analysis capability has been significantly enhanced and a much broader class of complex and more practical elastica problems can be solved. In past decades, the large displacement and rotation analysis based on exact kinematics has gained significant attention from various researchers and been used extensively to predict complex structural responses such as the postbuckling behavior. Some of those relevant studies are briefly summarized here not only to demonstrate the chronological progress and the recent advances in the area but also to indicate the key contribution of the current study. B. N. Rao and G. V. Rao [10] employed the fourth-order Runge-Kutta technique to solve for the large deflection of a cantilever beam subjected to either a rotating distributed load or a rotating concentrated load. Wang [11] applied a perturbation technique to investigate the postbuckling behavior of a single column with one of its ends being clamped and the other end being simply supported and subjected to an axial load. The postbuckling behavior of a prismatic cantilever column under the combined action between a uniformly distributed load and a concentrated load at the tip was later examined by Lee [12]. In such analysis, the numerical integration scheme based on Butcher’s fifth-order Runge-Kutta method was utilized to construct the numerical solutions. In 2002, Phungpaingam and Chucheepsakul [13] applied both the elliptic integral technique and the shooting method to analyze a simply supported beam of variable arc-length and subjected to an inclined follower force at any location within the member. Later, Vaz and Silva [14] utilized a two-parameter shooting method to generalize the work of Wang [11] by replacing the clamped end of a column by a rotational spring. Results from their study revealed that the rotational restraint at the end of the column significantly influences the postbsuckling configuration. Madhusudan et al. [15] extended the work of Lee [12] to explore the influence of a nonuniform cross-ssection on the postbuckling behavior of a cantilever column. The problem was formulated within the dynamic context, and the resulting nonlinear equations were solved by a fourth-order Runge-Kutta scheme. In 2007, Shvartsman [16] employed a technique of changing variables along with a solution scheme requiring no iteration to examine the influence of a rotational spring at the base of a cantilever beam and a tip-concentrated follower force on its deformed shape. Lacarbonara [17] applied the higher-order perturbation technique to determine postbuckling solutions for nonsprismatic nonlinearly elastic rods. Wang et al. [18] reexamined a cantilever beam for an explicit solution of the displacement and rotation at the free end by using a homotopy analytical method. Banerjee et al. [19] analyzed a cantilever beam subjected to arbitrary loading conditions and containing an interior inflection point by using a nonlinear shooting method and an adomain decomposition. Shvartsman [20] generalized the work of [16] to explore the influence of a variable cross-section and two follower forces on the behavior of a cantilever beam undergoing large displacement and rotation. Recently, Chen [21] proposed a moment integral scheme to solve for a large deflection of a cantilever beam. It was found from this study that the technique is computationally efficient, yields very accurate results, and is applicable to beams under various loading conditions and varying material and member properties. While there have been extensive studies related to large displacement and rotation analysis, all those mentioned above [10–21] are restricted mainly to structures consisting of only a single member.A comprehensive literature survey reveals that work focusing on the large displacement and rotation analysis or elastica of structures consisting of multiple members excluding those based on finite element approximations is still very limited. Some of those studies are summarized as follows. Ohtsuki and Ellyin [22] obtained an analytical solution in terms of elliptic integrals for flexural quantities (e.g., displacements, curvature, bending moment, etc.) of a square frame subjected to a pair of opposite nodal concentrated forces. Dado et al. [23] studied the postbuckling behavior of a cantilever beam consisting of two segments with different properties connected by an elastic rotational spring. Several methods including a semi-analytical technique, a numerical integration scheme, and a finite element method capable of modeling large displacement and rotation were employed and their performance was compared. Suwansheewasiri and Chucheepsakul [24] utilized the elliptic integral method to investigate both buckling and postbuckling of a two-member frame under nodal loads at its apex. Both symmetric and nonsymmetric postbuckling shapes were examined in their study. Dado and Al-Sadder [25] proposed a robust numerical technique based on an approximation of an angular deflection along the beam axis by a polynomial function to analyze a rhombus frame consisting of nonprismatic members and subjected to a pair of opposite nodal forces along its diagonal. Hu et al. [26] employed the differential quadrature element method (DQEM) to perform large displacement analysis of frames containing discontinuity conditions and initial displacements. While its computational efficiency and applicability to structures with general geometries were concluded, the technique itself is still an approximate scheme, and, as a result, the accuracy of numerical solutions depends primarily on the level of mesh refinement. Recently, Shatarat et al. [27] reinvestigated a rhombus nonlinear elastic frame under a pair of opposite nodal forces along its diagonal. In their study, a semi-analytical technique was utilized to obtain the relation between the displacement and applied forces at its corners. Based on existing literatures in the area as clearly indicated, the development of efficient and accurate techniques that are capable of modeling large displacement and rotation of structures with general geometries and loading conditions still deserves further investigations.In this paper, we propose a systematic and simple semi-analytical technique based on a direct stiffness method for large displacement and rotation analysis of linearly elastic, flexure-dominating skeleton structures ofarbitrary geometries and subjected to general nodal loads. The crucial feature of the current technique is the use of an exact element tangent stiffness matrix to form the best linear approximation of the governing nonlinear equations. Such local linear approximation along with the Newton-Ralphson iteration via the exact evaluation of residuals allows a semi-analytical solution to be obtained. It is worth noting that while the current approach and various existing techniques based on corotational finite element formulations (e.g., [28–32]) are closely related, the present study offers an alternative in which the exact form of governing equations is employed throughout the solution procedure, and this, as a result, yields results comparable to analytical solutions without refining the discretization. The accuracy and capability of the proposed technique are demonstrated via extensive numerical experiments.


## 2. Basic Equations

Basic assumptions and key components used to form a mathematical model for an individual member and to derive key differential equations governing its behavior include that (i) the member is prismatic and made of a homogeneous, isotropic, linearly elastic material; (ii) the curvature, displacement, and rotation are related by an exact kinematics; (iii) static equilibrium equations are enforced in the deformed state; (iv) member loads are absent; (v) the member is inextensible; (vi) shear deformation is negligible.Let us consider a straight, prismatic member of lengthL and moment of inertia I and made of an elastic material with Young’s modulus E. An undeformed configuration of this member occupies a straight line y=0,0≤x≤L, and its subsequent deformed state (resulting from a particular motion) is shown in Figure 1(a). It is important to remark first that the Lagrangion description is utilized throughout the following development, and, by supplying simplified kinematics of the cross-section (e.g., plane section always remains plane before and after undergoing deformation), the entire (three-dimensional) member can be sufficiently and completely represented by a single line connecting the centroid of all cross-sections (i.e., the neutral axis). From here to what follows, the phrases “cross-section at (S,0)” and “point (S,0)” are often utilized for convenience in this paper to refer to a cross-section at any state with its centroid located at a point (S,0) in the undeformed state.  During a particular motion, the cross-section at (S,0) in the undeformed state displaces to a point (S+u,v) in the deformed state where u=u(S) and v=v(S) denote the x-component and the y-component of the displacement at point (S,0), respectively. Let fx=fx(S),fy=fy(S), and m=m(S) denote a resultant internal force in x-direction, a resultant internal force in y-direction, and a resultant bending moment at the cross-section (S,0), respectively.(a) Schematic of deformed and undeformed configurations of member and (b) free body diagram of infinitesimal elementdS.
(a)(b)By enforcing static equilibrium of an infinitesimal elementdS in the deformed state (see its free body diagram in Figure 1(b)), it leads to three differential equations(2.1)dfxdS=0,(2.2)dfydS=0,(2.3)dmdS=fxsin⁡⁡θ+fycos⁡⁡θ,
where θ is the rotation of any cross-section. Clearly indicated by (2.1) and (2.2), the internal forces fx and fy must be constant for the entire member. It is evident from geometry of the element dS in the deformed state that the rotation θ can be related to the two components of the displacement u and v by(2.4)sin⁡⁡θ=dvdS,cos⁡⁡θ=1+dudS.
From the kinematics assumption (ii), the curvature κ and the rotation θ are related through κ=dθ/dS, and, by using assumptions (i) and (vi), we then obtain the linear moment-curvature relationship: m=EIκ. By using these two relations, (2.3) becomes(2.5)d2θdξ2=f̂xsin⁡⁡θ+f̂ycos⁡⁡θ,
where nondimensional parameters appearing in above equation are defined by ξ=S/L, f̂x=fxL2/EI, and f̂y=fyL2/EI. By performing a direct integration of (2.5), it leads to (2.6)(dθdξ)2=C-2f̂xcos⁡⁡θ+2f̂ysin⁡⁡θ,
where C is a constant of integration that can be determined from boundary conditions. It is worth noting that, from the moment-curvature relationship, sign of the normalized curvature dθ/dξ is uniquely dictated by that of the bending moment m. As a consequence, dξ/dθ can readily be solved from (2.6) to obtain(2.7)dξdθ=ϑ(m̂)C-2f̂xcos⁡⁡θ+2f̂ysin⁡⁡θ,
where ϑ(m̂) is a moment-dependent function defined such that ϑ(m̂)=1 if m̂=mL/EI>0 and ϑ(m̂)=-1 if m̂<0. By combining (2.4) and (2.7), it leads to following two relations:(2.8)dv̂dθ=ϑ(m̂)sin⁡⁡θC-2f̂xcos⁡⁡θ+2f̂ysin⁡⁡θ,dûdθ=ϑ(m̂)(cos⁡⁡θ-1)C-2f̂xcos⁡⁡θ+2f̂ysin⁡⁡θ,
where û=u/L and v̂=v/L. The three differential relations (2.7)-(2.8) constitute a basis for the development described below.


## 3. Governing Equations and Tangent Stiffness Matrix for Elements

In this section, we apply basic differential equations derived above to form a set of governing nonlinear algebraic equations and the exact expression of the tangent stiffness matrix of a two-dimensional member. To aid such development and attain anticipated outcomes, we first establish some useful results for a simply supported member and then use them along with the geometric consideration and the coordinate transformation.


### 3.1. Results for Simply Supported Member

Now, let us consider a simply-supported member (with a pinned support at its left end(S=0) and a roller support at its right end (S=L)) subjected to end loads {m1,m2,fx2} as shown in Figure 2. By imposing the moment boundary condition at the right end, that is, dθ/dξ(θ2)=m̂2=m2L/EI, we obtain the constant C as(3.1)C=m̂22+2f̂xcos⁡⁡θ2-2f̂ysin⁡⁡θ2.
By inserting the constant C into the relations (2.7)-(2.8), it yields(3.2)dξdθ=ϑ(m̂)F(θ,θ2;f̂x,f̂y,m̂2),dv̂dθ=ϑ(m̂)sin⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2),dûdθ=ϑ(m̂)(cos⁡⁡θ-1)F(θ,θ2;f̂x,f̂y,m̂2),
where the function F is defined by(3.3)F(θ,θ2;f̂x,f̂y,m̂2)=1m̂22+2f̂x(cos⁡⁡θ2-cos⁡⁡θ)-2f̂y(sin⁡⁡θ2-sin⁡⁡θ).
By imposing the remaining two natural boundary conditions at the left and right ends of a member, we obtain(3.4)m̂22-m̂12+2f̂x(cos⁡⁡θ2-cos⁡⁡θ1)-2f̂y(sin⁡⁡θ2-sin⁡⁡θ1)=0,(3.5)f̂x2-f̂x=0,
where m̂1=m1L/EI and f̂x2=fx2L2/EI. Normalized support reactions {f̂x1,f̂y1,f̂y2} can readily be computed from equilibrium of the entire member, and final results are given by(3.6)f̂x1=-f̂x,f̂y1=m̂1+m̂2d̂,f̂y2=-m̂1+m̂2d̂,
where d̂=1+û2 and û2=u2/L. Next, by integrating (3.2) from θ=θ1 to θ=θ2, it leads to a system of three nonlinear algebraic equations  (3.7)∫θ1θ2ϑ(m̂)F(θ,θ2;f̂x,f̂y,m̂2)dθ=1,(3.8)∫θ1θ2ϑ(m̂)sin⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ=0,(3.9)∫θ1θ2ϑ(m̂)cos⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ=d̂.
For a given set of end loads {f̂x2,m̂1,m̂2}, the end displacement and rotations {û2,θ1,θ2} can be solved from a system of nonlinear equations (3.7)–(3.9) with use of (3.4) and (3.5) to eliminate the internal forces {f̂x,f̂y}. On the other hand, the problem finding the end loads {f̂x2,m̂1,m̂2} for a particular set of prescribed displacement and rotations {û2,θ1,θ2} can possess no solution. Because of the geometric constraint imposed by the member inextensibility, the boundary value problem indicated above is not well-posed or, in other words, {û2,θ1,θ2} cannot be specified arbitrarily. However, if {f̂x,θ1,θ2} are prescribed instead, the end loads {f̂x2,m̂1,m̂2} can be obtained by solving (3.4), (3.5), (3.7), and (3.8) simultaneously, and the end displacement û2 can subsequently be computed from (3.9). However, lack of the displacement component û2 renders a set {f̂x,θ1,θ2} not well suited for the development of a solution procedure by a direct stiffness method.Figure 2
Schematic of simply supported member subjected to end loads{m1,m2,fx2}.In the present study,{û2,θ1,θ2,f̂x} is chosen as a set of primary unknowns. It is worth noting that to allow û2 to be one of independent variables, the strong requirement posed by (3.9) must be relaxed via the introduction of the residual ℛ such that(3.10)R≡d̂-∫θ1θ2ϑ(m̂)cos⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ.
Supplemented by (3.10), for any given set {û2,θ1,θ2,f̂x}, the quantities {f̂x2,m̂1,m̂2,R} can be obtained from (3.5), (3.7), (3.8), and (3.10) with use of (3.4) to get rid of f̂y. It should be emphasized that {û2,θ1,θ2,f̂x} and the corresponding {f̂x2,m̂1,m̂2,R} are solutions of the boundary value problem if and only if the residual ℛ vanishes.Letf be a vector defined by f=[fpfr]T where fp={f̂x2,m̂1,m̂2,R} and fr={f̂x1,f̂y1,f̂y2} and let u be a vector defined by u={û2,θ1,θ2,f̂x}T. From the relations (3.4)–(3.8) and (3.10), it can readily be verified that f=f(u), and, from Taylor series expansion, the best linear approximation of this nonlinear function f in the neighborhood of a vector u0 takes a form(3.11)f(u)=f(u0)+g(u0)(u-u0),
where a gradient matrix g is given by(3.12)g=∂f∂u=[gpgr].
The submatrix gp is expressed explicitly in terms of differentiations as(3.13)gp=[∂f̂x2∂û2∂f̂x2∂θ1∂f̂x2∂θ2∂f̂x2∂f̂x∂m̂1∂û2∂m̂1∂θ1∂m̂1∂θ2∂m̂1∂f̂x∂m̂2∂û2∂m̂2∂θ1∂m̂2∂θ2∂m̂2∂f̂x∂R∂û2∂R∂θ1∂R∂θ2∂R∂f̂x].
By defining gij as an entry located at the ith row and jth column of the submatrix gp, the submatrix gr can readily be obtained in terms of gij as(3.14)gr=1d̂[-g11d̂-g12d̂-g13d̂-g14d̂-ŝg22+g32g23+g33g24+g34ŝ-g22-g32-g23-g33-g24-g34],
where ŝ=(m̂1+m̂2)/d̂. As clearly indicated by (3.5), (3.7), (3.8), and (3.10), entries g11,g12,g13,g21, and g31 vanish whereas g14=g41=1. The remaining entries of gp are contained in a submatrix gp given by(3.15)gp=[∂m̂1∂θ1∂m̂1∂θ2∂m̂1∂f̂x∂m̂2∂θ1∂m̂2∂θ2∂m̂2∂f̂x∂R∂θ1∂R∂θ2∂R∂f̂x].
It should be remarked that determination of all entries of the submatrix gp is nontrivial and requires implicit differentiations.


### 3.2. Derivation of Submatrixℊp

The explicit form of the submatrixgp can be derived for various cases depending on the location of an inflection point within the member. For instance, a single-curvature member contains either no inflection point or an inflection point at its end whereas a double-curvature member contains an inflection point within the member. The key differences between these two cases are associated  with the value of the moment-dependent function ϑ(m̂) and the singularity behavior of the function F defined by (3.3).


#### 3.2.1. Member Containing No Inflection Point

The normalized bending moment is strictly positive (i.e.,m̂(ξ)>0) or strictly negative (i.e., m̂(ξ)<0) for the entire beam (i.e., ξ∈[0,1]) when the applied end moments {m̂1,m̂2} are nonzero and possess different sign (i.e., m̂1m̂2<0). The deformed configuration of the member for this particular case possesses a single curvature, and, in addition, ϑ(m̂) becomes a constant function with its value equal to either 1 or −1 depending on the sign of m̂; more specifically, ϑ(m̂)=1for m̂>0 and ϑ(m̂)=-1 for m̂<0. It is worth noting that the function F, defined by (3.3), is well behaved in the sense that the quantity within the square root is always greater than zero; this results directly from that m̂≠0 for the entire member. Such desirable feature of F renders all involved integrals nonsingular and, therefore, allows a standard procedure to be employed in their treatment. For convenience in further development, the relations (3.7), (3.8), and (3.10) are re-expressed as(3.16)Γ1(θ1,θ2;f̂x,f̂y,m̂2)=∫θ1θ2ϑ(m̂)F(θ,θ2;f̂x,f̂y,m̂2)dθ-1=0,(3.17)Γ2(θ1,θ2;f̂x,f̂y,m̂2)=∫θ1θ2ϑ(m̂)sin⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ=0,(3.18)R=R̂(θ1,θ2,û2;f̂x,f̂y,m̂2)=d̂-∫θ1θ2ϑ(m̂)cos⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ.
It is noted that the relation (3.4) implicitly defines f̂y in terms of {m̂1,m̂2} and {θ1,θ2,f̂x}. As a result, by taking derivative of (3.16)–(3.18) with respect to {θ1,θ2,f̂x} along with employing the chain rule of differentiation, it leads to(3.19)Sgp=-D,
where the matrices S and D are given by(3.20)S=[∂Γ1∂f̂y∂f̂y∂m̂1∂Γ1∂m̂2+∂Γ1∂f̂y∂f̂y∂m̂20∂Γ2∂f̂y∂f̂y∂m̂1∂Γ2∂m̂2+∂Γ2∂f̂y∂f̂y∂m̂20-∂R̂∂f̂y∂f̂y∂m̂1-∂R̂∂m̂2-∂R̂∂f̂y∂f̂y∂m̂21],(3.21)D=[∂Γ1∂θ1+∂Γ1∂f̂y∂f̂y∂θ1∂Γ1∂θ2+∂Γ1∂f̂y∂f̂y∂θ2∂Γ1∂f̂x+∂Γ1∂f̂y∂f̂y∂f̂x∂Γ2∂θ1+∂Γ2∂f̂y∂f̂y∂θ1∂Γ2∂θ2+∂Γ2∂f̂y∂f̂y∂θ2∂Γ2∂f̂x+∂Γ2∂f̂y∂f̂y∂f̂x-∂R̂∂θ1-∂R̂∂f̂y∂f̂y∂θ1-∂R̂∂θ2-∂R̂∂f̂y∂f̂y∂θ2-∂R̂∂f̂x-∂R̂∂f̂y∂f̂y∂f̂x].
Entries of the matrices S and D can be expressed in a closed form and the submatrix ℊp can then be obtained by solving (3.19) (see explicit results in Appendix A).


#### 3.2.2. Member Containing Interior Inflection Point

Now, let us consider a member containing an interior inflection point atξz∈(0,1) or, equivalently, the bending moment vanishing at ξz and m̂(ξ1)m̂(ξ2)<0 for ξ1∈[0,ξz) and ξ2∈(ξz,1]. This particular case arises when the applied end moments {m̂1,m̂2} are nonzero and possess the same sign (i.e., m̂1m̂2>0). The resulting deformed configuration of the member possesses a double-curvature shape. As a result, the moment-dependent function ϑ is discontinuous at ξz and takes different values on both sides of the inflection point. For the applied end moments m̂1,m̂2>0, it results in ϑ=-1 for ξ∈[0,ξz) and ϑ=1for ξ∈(ξz,1], and for m̂1,m̂2<0, it results in ϑ=1 for ξ∈[0,ξz) and ϑ=-1 for ξ∈(ξz,1].For this special case, the constantC appearing in (2.6) is alternatively obtained by imposing a condition at the inflection point, that is, dθ/dξ(θz)=0, and this leads to(3.22)C=2f̂xcos⁡⁡θz-2f̂ysin⁡⁡θz,
where θz is the rotation at the inflection point. Upon using (3.22), the relations (2.7)-(2.8) become(3.23)dξdθ=ϑ(m̂)Fz(θ,θz;f̂x,f̂y),dv̂dθ=ϑ(m̂)sin⁡⁡θFz(θ,θz;f̂x,f̂y),dûdθ=ϑ(m̂)(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y),
where the function Fz is defined by(3.24)Fz(θ,θz;f̂x,f̂y)=12f̂x(cos⁡⁡θz-cos⁡⁡θ)-2f̂y(sin⁡⁡θz-sin⁡⁡θ).
By integrating equations (3.23) over the entire member, we obtain(3.25)ψ[-∫θ1θzFz(θ,θz;f̂x,f̂y)dθ+∫θzθ2Fz(θ,θz;f̂x,f̂y)dθ]=1,ψ[-∫θ1θzsin⁡⁡θFz(θ,θz;f̂x,f̂y)dθ+∫θzθ2sin⁡⁡θFz(θ,θz;f̂x,f̂y)dθ]=0,ψ[-∫θ1θz(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y)dθ+∫θzθ2(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y)dθ]=û2,
where a constant ψ is defined such that ψ=1 if m̂1,m̂2>0 and ψ=-1 if m̂1,m̂2<0.It is evident from (3.24) that the function Fz is weakly singular at the inflection point (i.e., at θ=θz); as a result, all singular integrals appearing in equations (3.25) need special treatment. A series of variable transformations and some identities used to remove and regularize such singularity are summarized as follows: (i) introducing identities f̂s2=f̂x2+f̂y2,cosθo=f̂x/f̂s2, and sinθo=-f̂y/f̂s2, (ii) applying change of variable θ̅=π+(θ-θo) and identity cosθ̅=1-2sin2(θ̅/2), and (iii) introducing another variable transformation p̅sinϕ=sin(θ̅/2) with p̅=sin(θ̅z/2). The function Fz simply reduces to (3.26)Fz=12f̂s2[cos⁡⁡(θz-θo)-cos⁡⁡(θ-θo)]=14f̂s2[sin2(θ̅z/2)-sin2(θ̅/2)]=12f̂sp̅cos⁡⁡ϕ,
where θ̅z=π+(θz-θo). Finally, the nonlinear algebraic equations (3.25) can be expressed in terms of integrals over a dummy variable ϕ as(3.27)Γo(θ1,θ2,θz,f̂x,f̂y)=∫ϕ1π/2fo(ϕ;p̅)dϕ+∫ϕ2π/2fo(ϕ;p̅)dϕ-f̂s=0,(3.28)Γv(θ1,θ2,θz,f̂x,f̂y)=∫ϕ1π/2fv(ϕ,θo;p̅)dϕ+∫ϕ2π/2fv(ϕ,θo;p̅)dϕ=0,(3.29)R(θ1,θ2,û2,f̂x)=R̂(θ1,θ2,θz,û2,f̂x,f̂y)=û2-1f̂s[∫ϕ1π/2fu(ϕ,θo;p̅)dϕ+∫ϕ2π/2fu(ϕ,θo;p̅)dϕ],
where p̅sinϕ1=sin(θ̅1/2), p̅sinϕ2=sin(θ̅2/2), θ̅1=π+(θ1-θo), θ̅2=π+(θ2-θo), and (3.30)fo(ϕ;p̅)=11-p̅2sin2ϕ,fv(ϕ,θo;p̅)=sin⁡⁡θo1-p̅2sin2ϕ-2sin⁡⁡θo1-p̅2sin2ϕ-2ψp̅cos⁡⁡θosin⁡⁡ϕ,fu(ϕ,θo;p̅)=cos⁡⁡θo-11-p̅2sin2ϕ-2cos⁡⁡θo1-p̅2sin2ϕ+2ψp̅sin⁡⁡θosin⁡⁡ϕ.By enforcing the remaining two moment boundary conditions at both ends of the member, it leads to(3.31)m̂12+2f̂x(cos⁡⁡θ1-cos⁡⁡θz)-2f̂y(sin⁡⁡θ1-sin⁡⁡θz)=0,m̂22+2f̂x(cos⁡⁡θ2-cos⁡⁡θz)-2f̂y(sin⁡⁡θ2-sin⁡⁡θz)=0.
By taking derivative of (3.31) with respect to θ1, θ2, and f̂x, it results in(3.32)[∂m̂1∂θ1∂m̂1∂θ2∂m̂1∂f̂x∂m̂2∂θ1∂m̂2∂θ2∂m̂2∂f̂x]=M1+M2A,
where matrices M1,M2, and A are given by(3.33)M1=[λ1m̂10cz-c1m̂10λ2m̂2cz-c2m̂2],(3.34)M2=[-λzm̂1s1-szm̂1-λzm̂2s2-szm̂2],(3.35)A=[∂θz∂θ1∂θz∂θ2∂θz∂f̂x∂f̂y∂θ1∂f̂y∂θ2∂f̂y∂f̂x],
where s1=sinθ1, s2=sinθ2, sz=sinθz, c1=cosθ1, c2=cosθ2, cz=cosθz, λ1=f̂xs1+f̂yc1, λ2=f̂xs2+f̂yc2 and λz=f̂xsz+f̂ycz. Differentiating (3.29) with respect to θ1,θ2, and f̂x yields(3.36)[∂R∂θ1∂R∂θ2∂R∂f̂x]=B+CA,
where matrices B and C are given by(3.37)B=[∂R̂∂θ1∂R̂∂θ2∂R̂∂f̂x],C=[∂R̂∂θz∂R̂∂f̂y].
Note that all entries of the matrices B and C can be obtained directly from (3.29) along with the change of variables θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix B).To compute all entries of the matrixA, we differentiate (3.27) and (3.28) with respect to θ1,θ2, and f̂x, and this results in(3.38)DA=F,
where matrices D and F are given by(3.39)D=[∂Γo∂θz∂Γo∂f̂y∂Γv∂θz∂Γv∂f̂y],F=-[∂Γo∂θ1∂Γo∂θ2∂Γo∂f̂x∂Γv∂θ1∂Γv∂θ2∂Γv∂f̂x].
Similarly, all entries of the matrices D and F can be obtained directly from (3.27) and (3.28) along with the change of variables θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix B). Once the matrix A is solved from (3.38), it is inserted into (3.32) and (3.36) to obtain all entries of the matrix gp. Due to the complexity of all functions resulting from variable transformations, the matrices B,C,D, and F are evaluated numerically by standard Gaussian quadrature.


#### 3.2.3. Member Containing Inflection Point at Member End

Finally, consider a member containing an inflection point at one of its ends, or, equivalently, the bending moment possesses the same sign throughout the member and vanishes only at one of its ends. This particular case arises when one of the applied end moments{m̂1,m̂2} vanishes. The deformed configuration of the member possesses a single-curvature shape, and, in addition, the moment-dependent function ϑ becomes a constant function with its value equal to either 1 or −1 depending on the direction of the nonzero applied moments. Without loss of generality, the development presented below focuses only on the member containing an inflection point at its right end. While results for the member containing an inflection point at its left end are also needed, the treatment of such case follows the same procedure and is not presented here for brevity.Now, let us restrict our attention to a member subjected only to a nonzerom̂1. It should be noted that this particular member can be treated as a special case of a double-curvature member discussed in the Section 3.2.2 by simply taking θz=θ2 and ξz=1. As a consequence, basic equations and procedures adopted in the previous case can, after the proper specialization, be applied to this particular case. By replacing θz=θ2 into (3.25), it leads to(3.40)-ψ∫θ1θ2Fz1(θ,θ2;f̂x,f̂y)dθ=1,(3.41)-ψ∫θ1θ2sin⁡⁡θFz1(θ,θ2;f̂x,f̂y)dθ=0,(3.42)-ψ∫θ1θ2(cos⁡θ-1)Fz1(θ,θ2;f̂x,f̂y)dθ=û2,
where a function Fz1 is given by (3.43)Fz1(θ,θ2;f̂x,f̂y)=12f̂x(cos⁡⁡θ2-cos⁡⁡θ)-2f̂y(sin⁡⁡θ2-sin⁡⁡θ).By introducing the same type of variable transformations as employed in the previous case, (3.40)–(3.42) become(3.44)Γ̅o(θ1,θ2,f̂x,f̂y)=∫ϕ1π/2fo(ϕ;p̅)dϕ-f̂s=0,(3.45)Γ̅v(θ1,θ2,f̂x,f̂y)=∫ϕ1π/2fv(ϕ,θo;p̅)dϕ=0,(3.46)R̅(θ1,û2,f̂x)=R̅̂(θ1,θ2,û2,f̂x,f̂y)=û2-1f̂s∫ϕ1π/2fu(ϕ,θo;p̅)dϕ.
Since the right-end moment is prescribed (i.e., m̂2=0), the rotation at the right end θ2 is no longer an independent quantity but can be obtained in terms of {θ1,f̂x}.Let us redefinef̅=[f̅pf̅r]T where f̅p={f̂x2,m̂1,ℛ̅} and f̅r={f̂x1,f̂y1,f̂y2} and redefine u̅={û2,θ1,f̂x}T. Consistent with these new definitions, the reduced gradient matrix takes the form g̅=[g̅pTg̅rT]T where the submatrices g̅p and g̅r are given by(3.47)g̅p=[∂f̂x2∂û2∂f̂x2∂θ1∂f̂x2∂f̂x∂m̂1∂û2∂m̂1∂θ1∂m̂1∂f̂x∂R̅∂û2∂R̅∂θ1∂R̅∂f̂x],g̅r=1d̂[-g̅11d̂-g̅12d̂-g̅13d̂-s̅̂g̅22g̅23s̅̂-g̅22-g̅23],
where g̅ij denotes any entry of the submatrix g̅p and s̅̂=m̂1/d̂. As clearly indicated in the previous section, the entries g̅11,g̅12, and g̅21 vanish and g̅13=g̅31=1. The remaining entries are contained in a submatrix g̅p given by(3.48)g̅p=[∂m̂1∂θ1∂m̂1∂f̂x∂R̅∂θ1∂R̅∂f̂x].
By imposing the remaining moment boundary conditions at the left end of the member, we obtain(3.49)m̂12+2f̂x(cos⁡⁡θ1-cos⁡⁡θ2)-2f̂y(sin⁡⁡θ1-sin⁡⁡θ2)=0.Taking derivative of (3.49) with respect to θ1 and f̂x leads to(3.50)[∂m̂1∂θ1∂m̂1∂f̂x]=M̅1+M̅2A̅,
where matrices M̅1, M̅2, and A̅ are given by(3.51)M̅1=[λ1m̂1c2-c1m̂1],M̅2=[-λ2m̂1s1-s2m̂1],A̅=[∂θ2∂θ1∂θ2∂f̂x∂f̂y∂θ1∂f̂y∂f̂x].
By differentiating (3.46) with respect to θ1 and f̂x, it results in(3.52)[∂R̅∂θ1∂R̅∂f̂x]=B̅+C̅A̅,
where matrices B̅ and C̅ are(3.53)B̅=[∂R̅̂∂θ1∂R̅̂∂f̂x],C̅=[∂R̅̂∂θ2∂R̅̂∂f̂y].
Note that all entries of the matrices B̅ and C̅ can be obtained directly from (3.46) along with the transformations θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix C). To compute all entries of the matrix A̅, (3.44) and (3.45) are differentiated with respect to θ1 and f̂x and this yields(3.54)D̅A̅=F̅,
where matrices D̅ and F̅  are given by(3.55)D̅=[∂Γ̅o∂θ2∂Γ̅o∂f̂y∂Γ̅v∂θ2∂Γ̅v∂f̂y],F̅=-[∂Γ̅o∂θ1∂Γ̅o∂f̂x∂Γ̅v∂θ1∂Γ̅v∂f̂x].
Explicit expressions for all entries of the matrices D̅ and F̅ are reported in Appendix C. Once the matrix A̅ is solved from (3.54), it is substituted into (3.50) and (3.52) to obtain all entries of the matrix g̅p. Again, due to the complexity of all functions resulting from variable transformations, the matrices B̅,C̅,D̅, and F̅ are evaluated numerically using Gaussian quadrature.


### 3.3. Local Element Tangent Stiffness Matrix

Consider now a member with general boundary conditions as shown schematically in Figure3. Let {x,y} be a local coordinate system of the undeformed member and {x*,y*} be the coordinate system of the deformed member defined such that a chord connecting its end points always lies on the x* axis. With this specific choice of {x*,y*}, behavior of the member observed from this coordinate system is identical to that of a simply supported member.Figure 3
Schematic of undeformed and deformed configurations of a member subjected to general boundary conditions.The normalized end loads and normalized end displacements and rotations are denoted by{f̂x1,f̂y1,m̂1,f̂x2,f̂y2,m̂2} and {û1,v̂1,θ1,û2,v̂2,θ2} in the {x,y} coordinate system and by {f̂x1*,f̂y1*,m̂1*,f̂x2*,f̂y2*,m̂2*} and {û2*,θ1*,θ2*} in the {x*,y*} coordinate system. From geometric consideration of the deformed configuration, {û2*,θ1*,θ2*} can readily be expressed in terms of {û1,v̂1,θ1,û2,v̂2,θ2} by(3.56)θ1=θ1*+ϕ,θ2=θ2*+ϕ,û2*=(1+û2-û1)cos⁡⁡ϕ+(v̂2-v̂1)sin⁡⁡ϕ-1,
where ϕ is a chord rotation governed by(3.57)(1+û2-û1)sin⁡⁡ϕ-(v̂2-v̂1)cos⁡⁡ϕ=0.Letf̂x* be the internal force in x* direction and R* be the residual defined in the {x*,y*} system in the same way as given by (3.10). In the {x,y} coordinate system, we choose {f̂x,ℛ} such that(3.58)f̂x=f̂x*,R=R*.By applying the coordinate transformation, a vectorf={f̂x1,f̂y1,m̂1,f̂x2,f̂y2,m̂2,ℛ} is related to a vector f*={f̂x2*,m̂1*,m̂2*,ℛ*,f̂x1*,f̂y1*,f̂y2*} by(3.59)f=R(ϕ)f*,
where R(ϕ) is a transformation matrix of dimensions 7×7 given by(3.60)R(ϕ)=[0000cϕ-sϕ00000sϕcϕ00100000cϕ00000-sϕsϕ00000cϕ00100000001000],
in which sϕ=sinϕ and cϕ=cosϕ. By defining u={û1,v̂1,θ1,û2,v̂2,θ2,f̂x} and u*={û2*,θ1*,θ2*,f̂x*} and then recalling (3.56)–(3.58), we obtain the relation u*=u*(u). From the fact that behavior of the member in the {x*,y*} system is identical to that for a simply supported member, f* and u* are therefore related by f*=f*(u*). Combining (3.59), u*=u*(u) and f*=f*(u*) leads to the relation f=f(u)=R(ϕ)f*(u*(u)), and, from Taylor series expansion, this nonlinear function possesses the best linear approximation in the neighborhood of a vector u0(3.61)f(u)=f(u0)+kl(u0)(u-u0),
where kl is a local element tangent stiffness matrix given by(3.62)kl=∂R∂ϕf*∂ϕ∂u+Rg∂u*∂u,
in which the relation g=∂f*/∂u* is utilized.For a special case in which the member contains an inflection point at its right end, the end momentm̂2 vanishes and the corresponding end rotation θ2 is eliminated from the set of unknowns. Let us first define following reduced vectors f̅={f̂x1,f̂y1,m̂1,f̂x2,f̂y2,ℛ}, f̅*={f̂x2*,m̂1*,ℛ*,f̂x1*,f̂y1*,f̂y2*}, u̅={û1,v̂1,θ1,û2,v̂2,f̂x}, and u̅*={û2*,θ1*,f̂x*}. By applying the coordinate transformation, the relationship between the reduced vectors f̅ and f̅* is given by(3.63)f̅=R̅(ϕ)f̅*,
where R̅(ϕ) is a reduced transformation matrix of dimensions 6×6 given by(3.64)R̅(ϕ)=[000cϕ-sϕ0000sϕcϕ0010000cϕ0000-sϕsϕ0000cϕ001000].
The relation u̅*=u̅*(u̅) is then obtained by recalling (3.56)–(3.58) and f̅*=f̅*(u̅*) results from the fact that the behavior of the member in the {x*,y*} system is identical to that for a simply supported member. Combining (3.63), u̅*=u̅*(u̅) and f̅*=f̅*(u̅*) yields f̅=R̅(ϕ)f̅*(u̅*(u̅)), and, from Taylor series expansion, the best linear approximation of this nonlinear function in the neighborhood of a vector u̅0 takes the form(3.65)f̅(u̅)=f̅(u̅0)+k̅l(u̅0)(u̅-u̅0),
where k̅l is a reduced local element tangent stiffness matrix given by(3.66)k̅l=∂R̅∂ϕf̅*∂ϕ∂u̅+R̅g̅∂u̅*∂u̅,
in which the relation g̅=∂f̅*/∂u̅* is utilized. Note that the local element tangent stiffness matrix for this particular case is of dimensions 6×6.


### 3.4. Global Element Tangent Stiffness Matrix

Let the orientation of a member in the undeformed state be denoted by an angleβ measured from the global X-axis to the local x-axis (defined in the Section 3.3). The element tangent stiffness matrix kg referring to the global coordinate system {X,Y} can be obtained using the following standard coordinate transformation(3.67)kg=QTklQ,
where Q is a transformation matrix of dimension 7×7 given by(3.68)Q=[cβsβ00000-sβcβ000000010000000cβsβ00000-sβcβ0000000100000001],
in which sβ=sinβ and cβ=cosβ. Similarly, for a member containing an inflection point at the right end, we obtain(3.69)k̅g=Q̅Tk̅lQ̅,
where k̅g is the reduced, global element tangent stiffness matrix and Q̅ is a reduced transformation matrix given by(3.70)Q̅=[cβsβ0000-sβcβ0000001000000cβsβ0000-sβcβ0000001].


## 3.1. Results for Simply Supported Member

Now, let us consider a simply-supported member (with a pinned support at its left end(S=0) and a roller support at its right end (S=L)) subjected to end loads {m1,m2,fx2} as shown in Figure 2. By imposing the moment boundary condition at the right end, that is, dθ/dξ(θ2)=m̂2=m2L/EI, we obtain the constant C as(3.1)C=m̂22+2f̂xcos⁡⁡θ2-2f̂ysin⁡⁡θ2.
By inserting the constant C into the relations (2.7)-(2.8), it yields(3.2)dξdθ=ϑ(m̂)F(θ,θ2;f̂x,f̂y,m̂2),dv̂dθ=ϑ(m̂)sin⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2),dûdθ=ϑ(m̂)(cos⁡⁡θ-1)F(θ,θ2;f̂x,f̂y,m̂2),
where the function F is defined by(3.3)F(θ,θ2;f̂x,f̂y,m̂2)=1m̂22+2f̂x(cos⁡⁡θ2-cos⁡⁡θ)-2f̂y(sin⁡⁡θ2-sin⁡⁡θ).
By imposing the remaining two natural boundary conditions at the left and right ends of a member, we obtain(3.4)m̂22-m̂12+2f̂x(cos⁡⁡θ2-cos⁡⁡θ1)-2f̂y(sin⁡⁡θ2-sin⁡⁡θ1)=0,(3.5)f̂x2-f̂x=0,
where m̂1=m1L/EI and f̂x2=fx2L2/EI. Normalized support reactions {f̂x1,f̂y1,f̂y2} can readily be computed from equilibrium of the entire member, and final results are given by(3.6)f̂x1=-f̂x,f̂y1=m̂1+m̂2d̂,f̂y2=-m̂1+m̂2d̂,
where d̂=1+û2 and û2=u2/L. Next, by integrating (3.2) from θ=θ1 to θ=θ2, it leads to a system of three nonlinear algebraic equations  (3.7)∫θ1θ2ϑ(m̂)F(θ,θ2;f̂x,f̂y,m̂2)dθ=1,(3.8)∫θ1θ2ϑ(m̂)sin⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ=0,(3.9)∫θ1θ2ϑ(m̂)cos⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ=d̂.
For a given set of end loads {f̂x2,m̂1,m̂2}, the end displacement and rotations {û2,θ1,θ2} can be solved from a system of nonlinear equations (3.7)–(3.9) with use of (3.4) and (3.5) to eliminate the internal forces {f̂x,f̂y}. On the other hand, the problem finding the end loads {f̂x2,m̂1,m̂2} for a particular set of prescribed displacement and rotations {û2,θ1,θ2} can possess no solution. Because of the geometric constraint imposed by the member inextensibility, the boundary value problem indicated above is not well-posed or, in other words, {û2,θ1,θ2} cannot be specified arbitrarily. However, if {f̂x,θ1,θ2} are prescribed instead, the end loads {f̂x2,m̂1,m̂2} can be obtained by solving (3.4), (3.5), (3.7), and (3.8) simultaneously, and the end displacement û2 can subsequently be computed from (3.9). However, lack of the displacement component û2 renders a set {f̂x,θ1,θ2} not well suited for the development of a solution procedure by a direct stiffness method.Figure 2
Schematic of simply supported member subjected to end loads{m1,m2,fx2}.In the present study,{û2,θ1,θ2,f̂x} is chosen as a set of primary unknowns. It is worth noting that to allow û2 to be one of independent variables, the strong requirement posed by (3.9) must be relaxed via the introduction of the residual ℛ such that(3.10)R≡d̂-∫θ1θ2ϑ(m̂)cos⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ.
Supplemented by (3.10), for any given set {û2,θ1,θ2,f̂x}, the quantities {f̂x2,m̂1,m̂2,R} can be obtained from (3.5), (3.7), (3.8), and (3.10) with use of (3.4) to get rid of f̂y. It should be emphasized that {û2,θ1,θ2,f̂x} and the corresponding {f̂x2,m̂1,m̂2,R} are solutions of the boundary value problem if and only if the residual ℛ vanishes.Letf be a vector defined by f=[fpfr]T where fp={f̂x2,m̂1,m̂2,R} and fr={f̂x1,f̂y1,f̂y2} and let u be a vector defined by u={û2,θ1,θ2,f̂x}T. From the relations (3.4)–(3.8) and (3.10), it can readily be verified that f=f(u), and, from Taylor series expansion, the best linear approximation of this nonlinear function f in the neighborhood of a vector u0 takes a form(3.11)f(u)=f(u0)+g(u0)(u-u0),
where a gradient matrix g is given by(3.12)g=∂f∂u=[gpgr].
The submatrix gp is expressed explicitly in terms of differentiations as(3.13)gp=[∂f̂x2∂û2∂f̂x2∂θ1∂f̂x2∂θ2∂f̂x2∂f̂x∂m̂1∂û2∂m̂1∂θ1∂m̂1∂θ2∂m̂1∂f̂x∂m̂2∂û2∂m̂2∂θ1∂m̂2∂θ2∂m̂2∂f̂x∂R∂û2∂R∂θ1∂R∂θ2∂R∂f̂x].
By defining gij as an entry located at the ith row and jth column of the submatrix gp, the submatrix gr can readily be obtained in terms of gij as(3.14)gr=1d̂[-g11d̂-g12d̂-g13d̂-g14d̂-ŝg22+g32g23+g33g24+g34ŝ-g22-g32-g23-g33-g24-g34],
where ŝ=(m̂1+m̂2)/d̂. As clearly indicated by (3.5), (3.7), (3.8), and (3.10), entries g11,g12,g13,g21, and g31 vanish whereas g14=g41=1. The remaining entries of gp are contained in a submatrix gp given by(3.15)gp=[∂m̂1∂θ1∂m̂1∂θ2∂m̂1∂f̂x∂m̂2∂θ1∂m̂2∂θ2∂m̂2∂f̂x∂R∂θ1∂R∂θ2∂R∂f̂x].
It should be remarked that determination of all entries of the submatrix gp is nontrivial and requires implicit differentiations.


## 3.2. Derivation of Submatrixℊp

The explicit form of the submatrixgp can be derived for various cases depending on the location of an inflection point within the member. For instance, a single-curvature member contains either no inflection point or an inflection point at its end whereas a double-curvature member contains an inflection point within the member. The key differences between these two cases are associated  with the value of the moment-dependent function ϑ(m̂) and the singularity behavior of the function F defined by (3.3).


### 3.2.1. Member Containing No Inflection Point

The normalized bending moment is strictly positive (i.e.,m̂(ξ)>0) or strictly negative (i.e., m̂(ξ)<0) for the entire beam (i.e., ξ∈[0,1]) when the applied end moments {m̂1,m̂2} are nonzero and possess different sign (i.e., m̂1m̂2<0). The deformed configuration of the member for this particular case possesses a single curvature, and, in addition, ϑ(m̂) becomes a constant function with its value equal to either 1 or −1 depending on the sign of m̂; more specifically, ϑ(m̂)=1for m̂>0 and ϑ(m̂)=-1 for m̂<0. It is worth noting that the function F, defined by (3.3), is well behaved in the sense that the quantity within the square root is always greater than zero; this results directly from that m̂≠0 for the entire member. Such desirable feature of F renders all involved integrals nonsingular and, therefore, allows a standard procedure to be employed in their treatment. For convenience in further development, the relations (3.7), (3.8), and (3.10) are re-expressed as(3.16)Γ1(θ1,θ2;f̂x,f̂y,m̂2)=∫θ1θ2ϑ(m̂)F(θ,θ2;f̂x,f̂y,m̂2)dθ-1=0,(3.17)Γ2(θ1,θ2;f̂x,f̂y,m̂2)=∫θ1θ2ϑ(m̂)sin⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ=0,(3.18)R=R̂(θ1,θ2,û2;f̂x,f̂y,m̂2)=d̂-∫θ1θ2ϑ(m̂)cos⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ.
It is noted that the relation (3.4) implicitly defines f̂y in terms of {m̂1,m̂2} and {θ1,θ2,f̂x}. As a result, by taking derivative of (3.16)–(3.18) with respect to {θ1,θ2,f̂x} along with employing the chain rule of differentiation, it leads to(3.19)Sgp=-D,
where the matrices S and D are given by(3.20)S=[∂Γ1∂f̂y∂f̂y∂m̂1∂Γ1∂m̂2+∂Γ1∂f̂y∂f̂y∂m̂20∂Γ2∂f̂y∂f̂y∂m̂1∂Γ2∂m̂2+∂Γ2∂f̂y∂f̂y∂m̂20-∂R̂∂f̂y∂f̂y∂m̂1-∂R̂∂m̂2-∂R̂∂f̂y∂f̂y∂m̂21],(3.21)D=[∂Γ1∂θ1+∂Γ1∂f̂y∂f̂y∂θ1∂Γ1∂θ2+∂Γ1∂f̂y∂f̂y∂θ2∂Γ1∂f̂x+∂Γ1∂f̂y∂f̂y∂f̂x∂Γ2∂θ1+∂Γ2∂f̂y∂f̂y∂θ1∂Γ2∂θ2+∂Γ2∂f̂y∂f̂y∂θ2∂Γ2∂f̂x+∂Γ2∂f̂y∂f̂y∂f̂x-∂R̂∂θ1-∂R̂∂f̂y∂f̂y∂θ1-∂R̂∂θ2-∂R̂∂f̂y∂f̂y∂θ2-∂R̂∂f̂x-∂R̂∂f̂y∂f̂y∂f̂x].
Entries of the matrices S and D can be expressed in a closed form and the submatrix ℊp can then be obtained by solving (3.19) (see explicit results in Appendix A).


### 3.2.2. Member Containing Interior Inflection Point

Now, let us consider a member containing an interior inflection point atξz∈(0,1) or, equivalently, the bending moment vanishing at ξz and m̂(ξ1)m̂(ξ2)<0 for ξ1∈[0,ξz) and ξ2∈(ξz,1]. This particular case arises when the applied end moments {m̂1,m̂2} are nonzero and possess the same sign (i.e., m̂1m̂2>0). The resulting deformed configuration of the member possesses a double-curvature shape. As a result, the moment-dependent function ϑ is discontinuous at ξz and takes different values on both sides of the inflection point. For the applied end moments m̂1,m̂2>0, it results in ϑ=-1 for ξ∈[0,ξz) and ϑ=1for ξ∈(ξz,1], and for m̂1,m̂2<0, it results in ϑ=1 for ξ∈[0,ξz) and ϑ=-1 for ξ∈(ξz,1].For this special case, the constantC appearing in (2.6) is alternatively obtained by imposing a condition at the inflection point, that is, dθ/dξ(θz)=0, and this leads to(3.22)C=2f̂xcos⁡⁡θz-2f̂ysin⁡⁡θz,
where θz is the rotation at the inflection point. Upon using (3.22), the relations (2.7)-(2.8) become(3.23)dξdθ=ϑ(m̂)Fz(θ,θz;f̂x,f̂y),dv̂dθ=ϑ(m̂)sin⁡⁡θFz(θ,θz;f̂x,f̂y),dûdθ=ϑ(m̂)(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y),
where the function Fz is defined by(3.24)Fz(θ,θz;f̂x,f̂y)=12f̂x(cos⁡⁡θz-cos⁡⁡θ)-2f̂y(sin⁡⁡θz-sin⁡⁡θ).
By integrating equations (3.23) over the entire member, we obtain(3.25)ψ[-∫θ1θzFz(θ,θz;f̂x,f̂y)dθ+∫θzθ2Fz(θ,θz;f̂x,f̂y)dθ]=1,ψ[-∫θ1θzsin⁡⁡θFz(θ,θz;f̂x,f̂y)dθ+∫θzθ2sin⁡⁡θFz(θ,θz;f̂x,f̂y)dθ]=0,ψ[-∫θ1θz(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y)dθ+∫θzθ2(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y)dθ]=û2,
where a constant ψ is defined such that ψ=1 if m̂1,m̂2>0 and ψ=-1 if m̂1,m̂2<0.It is evident from (3.24) that the function Fz is weakly singular at the inflection point (i.e., at θ=θz); as a result, all singular integrals appearing in equations (3.25) need special treatment. A series of variable transformations and some identities used to remove and regularize such singularity are summarized as follows: (i) introducing identities f̂s2=f̂x2+f̂y2,cosθo=f̂x/f̂s2, and sinθo=-f̂y/f̂s2, (ii) applying change of variable θ̅=π+(θ-θo) and identity cosθ̅=1-2sin2(θ̅/2), and (iii) introducing another variable transformation p̅sinϕ=sin(θ̅/2) with p̅=sin(θ̅z/2). The function Fz simply reduces to (3.26)Fz=12f̂s2[cos⁡⁡(θz-θo)-cos⁡⁡(θ-θo)]=14f̂s2[sin2(θ̅z/2)-sin2(θ̅/2)]=12f̂sp̅cos⁡⁡ϕ,
where θ̅z=π+(θz-θo). Finally, the nonlinear algebraic equations (3.25) can be expressed in terms of integrals over a dummy variable ϕ as(3.27)Γo(θ1,θ2,θz,f̂x,f̂y)=∫ϕ1π/2fo(ϕ;p̅)dϕ+∫ϕ2π/2fo(ϕ;p̅)dϕ-f̂s=0,(3.28)Γv(θ1,θ2,θz,f̂x,f̂y)=∫ϕ1π/2fv(ϕ,θo;p̅)dϕ+∫ϕ2π/2fv(ϕ,θo;p̅)dϕ=0,(3.29)R(θ1,θ2,û2,f̂x)=R̂(θ1,θ2,θz,û2,f̂x,f̂y)=û2-1f̂s[∫ϕ1π/2fu(ϕ,θo;p̅)dϕ+∫ϕ2π/2fu(ϕ,θo;p̅)dϕ],
where p̅sinϕ1=sin(θ̅1/2), p̅sinϕ2=sin(θ̅2/2), θ̅1=π+(θ1-θo), θ̅2=π+(θ2-θo), and (3.30)fo(ϕ;p̅)=11-p̅2sin2ϕ,fv(ϕ,θo;p̅)=sin⁡⁡θo1-p̅2sin2ϕ-2sin⁡⁡θo1-p̅2sin2ϕ-2ψp̅cos⁡⁡θosin⁡⁡ϕ,fu(ϕ,θo;p̅)=cos⁡⁡θo-11-p̅2sin2ϕ-2cos⁡⁡θo1-p̅2sin2ϕ+2ψp̅sin⁡⁡θosin⁡⁡ϕ.By enforcing the remaining two moment boundary conditions at both ends of the member, it leads to(3.31)m̂12+2f̂x(cos⁡⁡θ1-cos⁡⁡θz)-2f̂y(sin⁡⁡θ1-sin⁡⁡θz)=0,m̂22+2f̂x(cos⁡⁡θ2-cos⁡⁡θz)-2f̂y(sin⁡⁡θ2-sin⁡⁡θz)=0.
By taking derivative of (3.31) with respect to θ1, θ2, and f̂x, it results in(3.32)[∂m̂1∂θ1∂m̂1∂θ2∂m̂1∂f̂x∂m̂2∂θ1∂m̂2∂θ2∂m̂2∂f̂x]=M1+M2A,
where matrices M1,M2, and A are given by(3.33)M1=[λ1m̂10cz-c1m̂10λ2m̂2cz-c2m̂2],(3.34)M2=[-λzm̂1s1-szm̂1-λzm̂2s2-szm̂2],(3.35)A=[∂θz∂θ1∂θz∂θ2∂θz∂f̂x∂f̂y∂θ1∂f̂y∂θ2∂f̂y∂f̂x],
where s1=sinθ1, s2=sinθ2, sz=sinθz, c1=cosθ1, c2=cosθ2, cz=cosθz, λ1=f̂xs1+f̂yc1, λ2=f̂xs2+f̂yc2 and λz=f̂xsz+f̂ycz. Differentiating (3.29) with respect to θ1,θ2, and f̂x yields(3.36)[∂R∂θ1∂R∂θ2∂R∂f̂x]=B+CA,
where matrices B and C are given by(3.37)B=[∂R̂∂θ1∂R̂∂θ2∂R̂∂f̂x],C=[∂R̂∂θz∂R̂∂f̂y].
Note that all entries of the matrices B and C can be obtained directly from (3.29) along with the change of variables θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix B).To compute all entries of the matrixA, we differentiate (3.27) and (3.28) with respect to θ1,θ2, and f̂x, and this results in(3.38)DA=F,
where matrices D and F are given by(3.39)D=[∂Γo∂θz∂Γo∂f̂y∂Γv∂θz∂Γv∂f̂y],F=-[∂Γo∂θ1∂Γo∂θ2∂Γo∂f̂x∂Γv∂θ1∂Γv∂θ2∂Γv∂f̂x].
Similarly, all entries of the matrices D and F can be obtained directly from (3.27) and (3.28) along with the change of variables θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix B). Once the matrix A is solved from (3.38), it is inserted into (3.32) and (3.36) to obtain all entries of the matrix gp. Due to the complexity of all functions resulting from variable transformations, the matrices B,C,D, and F are evaluated numerically by standard Gaussian quadrature.


### 3.2.3. Member Containing Inflection Point at Member End

Finally, consider a member containing an inflection point at one of its ends, or, equivalently, the bending moment possesses the same sign throughout the member and vanishes only at one of its ends. This particular case arises when one of the applied end moments{m̂1,m̂2} vanishes. The deformed configuration of the member possesses a single-curvature shape, and, in addition, the moment-dependent function ϑ becomes a constant function with its value equal to either 1 or −1 depending on the direction of the nonzero applied moments. Without loss of generality, the development presented below focuses only on the member containing an inflection point at its right end. While results for the member containing an inflection point at its left end are also needed, the treatment of such case follows the same procedure and is not presented here for brevity.Now, let us restrict our attention to a member subjected only to a nonzerom̂1. It should be noted that this particular member can be treated as a special case of a double-curvature member discussed in the Section 3.2.2 by simply taking θz=θ2 and ξz=1. As a consequence, basic equations and procedures adopted in the previous case can, after the proper specialization, be applied to this particular case. By replacing θz=θ2 into (3.25), it leads to(3.40)-ψ∫θ1θ2Fz1(θ,θ2;f̂x,f̂y)dθ=1,(3.41)-ψ∫θ1θ2sin⁡⁡θFz1(θ,θ2;f̂x,f̂y)dθ=0,(3.42)-ψ∫θ1θ2(cos⁡θ-1)Fz1(θ,θ2;f̂x,f̂y)dθ=û2,
where a function Fz1 is given by (3.43)Fz1(θ,θ2;f̂x,f̂y)=12f̂x(cos⁡⁡θ2-cos⁡⁡θ)-2f̂y(sin⁡⁡θ2-sin⁡⁡θ).By introducing the same type of variable transformations as employed in the previous case, (3.40)–(3.42) become(3.44)Γ̅o(θ1,θ2,f̂x,f̂y)=∫ϕ1π/2fo(ϕ;p̅)dϕ-f̂s=0,(3.45)Γ̅v(θ1,θ2,f̂x,f̂y)=∫ϕ1π/2fv(ϕ,θo;p̅)dϕ=0,(3.46)R̅(θ1,û2,f̂x)=R̅̂(θ1,θ2,û2,f̂x,f̂y)=û2-1f̂s∫ϕ1π/2fu(ϕ,θo;p̅)dϕ.
Since the right-end moment is prescribed (i.e., m̂2=0), the rotation at the right end θ2 is no longer an independent quantity but can be obtained in terms of {θ1,f̂x}.Let us redefinef̅=[f̅pf̅r]T where f̅p={f̂x2,m̂1,ℛ̅} and f̅r={f̂x1,f̂y1,f̂y2} and redefine u̅={û2,θ1,f̂x}T. Consistent with these new definitions, the reduced gradient matrix takes the form g̅=[g̅pTg̅rT]T where the submatrices g̅p and g̅r are given by(3.47)g̅p=[∂f̂x2∂û2∂f̂x2∂θ1∂f̂x2∂f̂x∂m̂1∂û2∂m̂1∂θ1∂m̂1∂f̂x∂R̅∂û2∂R̅∂θ1∂R̅∂f̂x],g̅r=1d̂[-g̅11d̂-g̅12d̂-g̅13d̂-s̅̂g̅22g̅23s̅̂-g̅22-g̅23],
where g̅ij denotes any entry of the submatrix g̅p and s̅̂=m̂1/d̂. As clearly indicated in the previous section, the entries g̅11,g̅12, and g̅21 vanish and g̅13=g̅31=1. The remaining entries are contained in a submatrix g̅p given by(3.48)g̅p=[∂m̂1∂θ1∂m̂1∂f̂x∂R̅∂θ1∂R̅∂f̂x].
By imposing the remaining moment boundary conditions at the left end of the member, we obtain(3.49)m̂12+2f̂x(cos⁡⁡θ1-cos⁡⁡θ2)-2f̂y(sin⁡⁡θ1-sin⁡⁡θ2)=0.Taking derivative of (3.49) with respect to θ1 and f̂x leads to(3.50)[∂m̂1∂θ1∂m̂1∂f̂x]=M̅1+M̅2A̅,
where matrices M̅1, M̅2, and A̅ are given by(3.51)M̅1=[λ1m̂1c2-c1m̂1],M̅2=[-λ2m̂1s1-s2m̂1],A̅=[∂θ2∂θ1∂θ2∂f̂x∂f̂y∂θ1∂f̂y∂f̂x].
By differentiating (3.46) with respect to θ1 and f̂x, it results in(3.52)[∂R̅∂θ1∂R̅∂f̂x]=B̅+C̅A̅,
where matrices B̅ and C̅ are(3.53)B̅=[∂R̅̂∂θ1∂R̅̂∂f̂x],C̅=[∂R̅̂∂θ2∂R̅̂∂f̂y].
Note that all entries of the matrices B̅ and C̅ can be obtained directly from (3.46) along with the transformations θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix C). To compute all entries of the matrix A̅, (3.44) and (3.45) are differentiated with respect to θ1 and f̂x and this yields(3.54)D̅A̅=F̅,
where matrices D̅ and F̅  are given by(3.55)D̅=[∂Γ̅o∂θ2∂Γ̅o∂f̂y∂Γ̅v∂θ2∂Γ̅v∂f̂y],F̅=-[∂Γ̅o∂θ1∂Γ̅o∂f̂x∂Γ̅v∂θ1∂Γ̅v∂f̂x].
Explicit expressions for all entries of the matrices D̅ and F̅ are reported in Appendix C. Once the matrix A̅ is solved from (3.54), it is substituted into (3.50) and (3.52) to obtain all entries of the matrix g̅p. Again, due to the complexity of all functions resulting from variable transformations, the matrices B̅,C̅,D̅, and F̅ are evaluated numerically using Gaussian quadrature.


## 3.2.1. Member Containing No Inflection Point

The normalized bending moment is strictly positive (i.e.,m̂(ξ)>0) or strictly negative (i.e., m̂(ξ)<0) for the entire beam (i.e., ξ∈[0,1]) when the applied end moments {m̂1,m̂2} are nonzero and possess different sign (i.e., m̂1m̂2<0). The deformed configuration of the member for this particular case possesses a single curvature, and, in addition, ϑ(m̂) becomes a constant function with its value equal to either 1 or −1 depending on the sign of m̂; more specifically, ϑ(m̂)=1for m̂>0 and ϑ(m̂)=-1 for m̂<0. It is worth noting that the function F, defined by (3.3), is well behaved in the sense that the quantity within the square root is always greater than zero; this results directly from that m̂≠0 for the entire member. Such desirable feature of F renders all involved integrals nonsingular and, therefore, allows a standard procedure to be employed in their treatment. For convenience in further development, the relations (3.7), (3.8), and (3.10) are re-expressed as(3.16)Γ1(θ1,θ2;f̂x,f̂y,m̂2)=∫θ1θ2ϑ(m̂)F(θ,θ2;f̂x,f̂y,m̂2)dθ-1=0,(3.17)Γ2(θ1,θ2;f̂x,f̂y,m̂2)=∫θ1θ2ϑ(m̂)sin⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ=0,(3.18)R=R̂(θ1,θ2,û2;f̂x,f̂y,m̂2)=d̂-∫θ1θ2ϑ(m̂)cos⁡⁡θF(θ,θ2;f̂x,f̂y,m̂2)dθ.
It is noted that the relation (3.4) implicitly defines f̂y in terms of {m̂1,m̂2} and {θ1,θ2,f̂x}. As a result, by taking derivative of (3.16)–(3.18) with respect to {θ1,θ2,f̂x} along with employing the chain rule of differentiation, it leads to(3.19)Sgp=-D,
where the matrices S and D are given by(3.20)S=[∂Γ1∂f̂y∂f̂y∂m̂1∂Γ1∂m̂2+∂Γ1∂f̂y∂f̂y∂m̂20∂Γ2∂f̂y∂f̂y∂m̂1∂Γ2∂m̂2+∂Γ2∂f̂y∂f̂y∂m̂20-∂R̂∂f̂y∂f̂y∂m̂1-∂R̂∂m̂2-∂R̂∂f̂y∂f̂y∂m̂21],(3.21)D=[∂Γ1∂θ1+∂Γ1∂f̂y∂f̂y∂θ1∂Γ1∂θ2+∂Γ1∂f̂y∂f̂y∂θ2∂Γ1∂f̂x+∂Γ1∂f̂y∂f̂y∂f̂x∂Γ2∂θ1+∂Γ2∂f̂y∂f̂y∂θ1∂Γ2∂θ2+∂Γ2∂f̂y∂f̂y∂θ2∂Γ2∂f̂x+∂Γ2∂f̂y∂f̂y∂f̂x-∂R̂∂θ1-∂R̂∂f̂y∂f̂y∂θ1-∂R̂∂θ2-∂R̂∂f̂y∂f̂y∂θ2-∂R̂∂f̂x-∂R̂∂f̂y∂f̂y∂f̂x].
Entries of the matrices S and D can be expressed in a closed form and the submatrix ℊp can then be obtained by solving (3.19) (see explicit results in Appendix A).


## 3.2.2. Member Containing Interior Inflection Point

Now, let us consider a member containing an interior inflection point atξz∈(0,1) or, equivalently, the bending moment vanishing at ξz and m̂(ξ1)m̂(ξ2)<0 for ξ1∈[0,ξz) and ξ2∈(ξz,1]. This particular case arises when the applied end moments {m̂1,m̂2} are nonzero and possess the same sign (i.e., m̂1m̂2>0). The resulting deformed configuration of the member possesses a double-curvature shape. As a result, the moment-dependent function ϑ is discontinuous at ξz and takes different values on both sides of the inflection point. For the applied end moments m̂1,m̂2>0, it results in ϑ=-1 for ξ∈[0,ξz) and ϑ=1for ξ∈(ξz,1], and for m̂1,m̂2<0, it results in ϑ=1 for ξ∈[0,ξz) and ϑ=-1 for ξ∈(ξz,1].For this special case, the constantC appearing in (2.6) is alternatively obtained by imposing a condition at the inflection point, that is, dθ/dξ(θz)=0, and this leads to(3.22)C=2f̂xcos⁡⁡θz-2f̂ysin⁡⁡θz,
where θz is the rotation at the inflection point. Upon using (3.22), the relations (2.7)-(2.8) become(3.23)dξdθ=ϑ(m̂)Fz(θ,θz;f̂x,f̂y),dv̂dθ=ϑ(m̂)sin⁡⁡θFz(θ,θz;f̂x,f̂y),dûdθ=ϑ(m̂)(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y),
where the function Fz is defined by(3.24)Fz(θ,θz;f̂x,f̂y)=12f̂x(cos⁡⁡θz-cos⁡⁡θ)-2f̂y(sin⁡⁡θz-sin⁡⁡θ).
By integrating equations (3.23) over the entire member, we obtain(3.25)ψ[-∫θ1θzFz(θ,θz;f̂x,f̂y)dθ+∫θzθ2Fz(θ,θz;f̂x,f̂y)dθ]=1,ψ[-∫θ1θzsin⁡⁡θFz(θ,θz;f̂x,f̂y)dθ+∫θzθ2sin⁡⁡θFz(θ,θz;f̂x,f̂y)dθ]=0,ψ[-∫θ1θz(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y)dθ+∫θzθ2(cos⁡⁡θ-1)Fz(θ,θz;f̂x,f̂y)dθ]=û2,
where a constant ψ is defined such that ψ=1 if m̂1,m̂2>0 and ψ=-1 if m̂1,m̂2<0.It is evident from (3.24) that the function Fz is weakly singular at the inflection point (i.e., at θ=θz); as a result, all singular integrals appearing in equations (3.25) need special treatment. A series of variable transformations and some identities used to remove and regularize such singularity are summarized as follows: (i) introducing identities f̂s2=f̂x2+f̂y2,cosθo=f̂x/f̂s2, and sinθo=-f̂y/f̂s2, (ii) applying change of variable θ̅=π+(θ-θo) and identity cosθ̅=1-2sin2(θ̅/2), and (iii) introducing another variable transformation p̅sinϕ=sin(θ̅/2) with p̅=sin(θ̅z/2). The function Fz simply reduces to (3.26)Fz=12f̂s2[cos⁡⁡(θz-θo)-cos⁡⁡(θ-θo)]=14f̂s2[sin2(θ̅z/2)-sin2(θ̅/2)]=12f̂sp̅cos⁡⁡ϕ,
where θ̅z=π+(θz-θo). Finally, the nonlinear algebraic equations (3.25) can be expressed in terms of integrals over a dummy variable ϕ as(3.27)Γo(θ1,θ2,θz,f̂x,f̂y)=∫ϕ1π/2fo(ϕ;p̅)dϕ+∫ϕ2π/2fo(ϕ;p̅)dϕ-f̂s=0,(3.28)Γv(θ1,θ2,θz,f̂x,f̂y)=∫ϕ1π/2fv(ϕ,θo;p̅)dϕ+∫ϕ2π/2fv(ϕ,θo;p̅)dϕ=0,(3.29)R(θ1,θ2,û2,f̂x)=R̂(θ1,θ2,θz,û2,f̂x,f̂y)=û2-1f̂s[∫ϕ1π/2fu(ϕ,θo;p̅)dϕ+∫ϕ2π/2fu(ϕ,θo;p̅)dϕ],
where p̅sinϕ1=sin(θ̅1/2), p̅sinϕ2=sin(θ̅2/2), θ̅1=π+(θ1-θo), θ̅2=π+(θ2-θo), and (3.30)fo(ϕ;p̅)=11-p̅2sin2ϕ,fv(ϕ,θo;p̅)=sin⁡⁡θo1-p̅2sin2ϕ-2sin⁡⁡θo1-p̅2sin2ϕ-2ψp̅cos⁡⁡θosin⁡⁡ϕ,fu(ϕ,θo;p̅)=cos⁡⁡θo-11-p̅2sin2ϕ-2cos⁡⁡θo1-p̅2sin2ϕ+2ψp̅sin⁡⁡θosin⁡⁡ϕ.By enforcing the remaining two moment boundary conditions at both ends of the member, it leads to(3.31)m̂12+2f̂x(cos⁡⁡θ1-cos⁡⁡θz)-2f̂y(sin⁡⁡θ1-sin⁡⁡θz)=0,m̂22+2f̂x(cos⁡⁡θ2-cos⁡⁡θz)-2f̂y(sin⁡⁡θ2-sin⁡⁡θz)=0.
By taking derivative of (3.31) with respect to θ1, θ2, and f̂x, it results in(3.32)[∂m̂1∂θ1∂m̂1∂θ2∂m̂1∂f̂x∂m̂2∂θ1∂m̂2∂θ2∂m̂2∂f̂x]=M1+M2A,
where matrices M1,M2, and A are given by(3.33)M1=[λ1m̂10cz-c1m̂10λ2m̂2cz-c2m̂2],(3.34)M2=[-λzm̂1s1-szm̂1-λzm̂2s2-szm̂2],(3.35)A=[∂θz∂θ1∂θz∂θ2∂θz∂f̂x∂f̂y∂θ1∂f̂y∂θ2∂f̂y∂f̂x],
where s1=sinθ1, s2=sinθ2, sz=sinθz, c1=cosθ1, c2=cosθ2, cz=cosθz, λ1=f̂xs1+f̂yc1, λ2=f̂xs2+f̂yc2 and λz=f̂xsz+f̂ycz. Differentiating (3.29) with respect to θ1,θ2, and f̂x yields(3.36)[∂R∂θ1∂R∂θ2∂R∂f̂x]=B+CA,
where matrices B and C are given by(3.37)B=[∂R̂∂θ1∂R̂∂θ2∂R̂∂f̂x],C=[∂R̂∂θz∂R̂∂f̂y].
Note that all entries of the matrices B and C can be obtained directly from (3.29) along with the change of variables θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix B).To compute all entries of the matrixA, we differentiate (3.27) and (3.28) with respect to θ1,θ2, and f̂x, and this results in(3.38)DA=F,
where matrices D and F are given by(3.39)D=[∂Γo∂θz∂Γo∂f̂y∂Γv∂θz∂Γv∂f̂y],F=-[∂Γo∂θ1∂Γo∂θ2∂Γo∂f̂x∂Γv∂θ1∂Γv∂θ2∂Γv∂f̂x].
Similarly, all entries of the matrices D and F can be obtained directly from (3.27) and (3.28) along with the change of variables θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix B). Once the matrix A is solved from (3.38), it is inserted into (3.32) and (3.36) to obtain all entries of the matrix gp. Due to the complexity of all functions resulting from variable transformations, the matrices B,C,D, and F are evaluated numerically by standard Gaussian quadrature.


## 3.2.3. Member Containing Inflection Point at Member End

Finally, consider a member containing an inflection point at one of its ends, or, equivalently, the bending moment possesses the same sign throughout the member and vanishes only at one of its ends. This particular case arises when one of the applied end moments{m̂1,m̂2} vanishes. The deformed configuration of the member possesses a single-curvature shape, and, in addition, the moment-dependent function ϑ becomes a constant function with its value equal to either 1 or −1 depending on the direction of the nonzero applied moments. Without loss of generality, the development presented below focuses only on the member containing an inflection point at its right end. While results for the member containing an inflection point at its left end are also needed, the treatment of such case follows the same procedure and is not presented here for brevity.Now, let us restrict our attention to a member subjected only to a nonzerom̂1. It should be noted that this particular member can be treated as a special case of a double-curvature member discussed in the Section 3.2.2 by simply taking θz=θ2 and ξz=1. As a consequence, basic equations and procedures adopted in the previous case can, after the proper specialization, be applied to this particular case. By replacing θz=θ2 into (3.25), it leads to(3.40)-ψ∫θ1θ2Fz1(θ,θ2;f̂x,f̂y)dθ=1,(3.41)-ψ∫θ1θ2sin⁡⁡θFz1(θ,θ2;f̂x,f̂y)dθ=0,(3.42)-ψ∫θ1θ2(cos⁡θ-1)Fz1(θ,θ2;f̂x,f̂y)dθ=û2,
where a function Fz1 is given by (3.43)Fz1(θ,θ2;f̂x,f̂y)=12f̂x(cos⁡⁡θ2-cos⁡⁡θ)-2f̂y(sin⁡⁡θ2-sin⁡⁡θ).By introducing the same type of variable transformations as employed in the previous case, (3.40)–(3.42) become(3.44)Γ̅o(θ1,θ2,f̂x,f̂y)=∫ϕ1π/2fo(ϕ;p̅)dϕ-f̂s=0,(3.45)Γ̅v(θ1,θ2,f̂x,f̂y)=∫ϕ1π/2fv(ϕ,θo;p̅)dϕ=0,(3.46)R̅(θ1,û2,f̂x)=R̅̂(θ1,θ2,û2,f̂x,f̂y)=û2-1f̂s∫ϕ1π/2fu(ϕ,θo;p̅)dϕ.
Since the right-end moment is prescribed (i.e., m̂2=0), the rotation at the right end θ2 is no longer an independent quantity but can be obtained in terms of {θ1,f̂x}.Let us redefinef̅=[f̅pf̅r]T where f̅p={f̂x2,m̂1,ℛ̅} and f̅r={f̂x1,f̂y1,f̂y2} and redefine u̅={û2,θ1,f̂x}T. Consistent with these new definitions, the reduced gradient matrix takes the form g̅=[g̅pTg̅rT]T where the submatrices g̅p and g̅r are given by(3.47)g̅p=[∂f̂x2∂û2∂f̂x2∂θ1∂f̂x2∂f̂x∂m̂1∂û2∂m̂1∂θ1∂m̂1∂f̂x∂R̅∂û2∂R̅∂θ1∂R̅∂f̂x],g̅r=1d̂[-g̅11d̂-g̅12d̂-g̅13d̂-s̅̂g̅22g̅23s̅̂-g̅22-g̅23],
where g̅ij denotes any entry of the submatrix g̅p and s̅̂=m̂1/d̂. As clearly indicated in the previous section, the entries g̅11,g̅12, and g̅21 vanish and g̅13=g̅31=1. The remaining entries are contained in a submatrix g̅p given by(3.48)g̅p=[∂m̂1∂θ1∂m̂1∂f̂x∂R̅∂θ1∂R̅∂f̂x].
By imposing the remaining moment boundary conditions at the left end of the member, we obtain(3.49)m̂12+2f̂x(cos⁡⁡θ1-cos⁡⁡θ2)-2f̂y(sin⁡⁡θ1-sin⁡⁡θ2)=0.Taking derivative of (3.49) with respect to θ1 and f̂x leads to(3.50)[∂m̂1∂θ1∂m̂1∂f̂x]=M̅1+M̅2A̅,
where matrices M̅1, M̅2, and A̅ are given by(3.51)M̅1=[λ1m̂1c2-c1m̂1],M̅2=[-λ2m̂1s1-s2m̂1],A̅=[∂θ2∂θ1∂θ2∂f̂x∂f̂y∂θ1∂f̂y∂f̂x].
By differentiating (3.46) with respect to θ1 and f̂x, it results in(3.52)[∂R̅∂θ1∂R̅∂f̂x]=B̅+C̅A̅,
where matrices B̅ and C̅ are(3.53)B̅=[∂R̅̂∂θ1∂R̅̂∂f̂x],C̅=[∂R̅̂∂θ2∂R̅̂∂f̂y].
Note that all entries of the matrices B̅ and C̅ can be obtained directly from (3.46) along with the transformations θ̅=π+(θ-θo) and p̅sinϕ=sin(θ̅/2) (see explicit results in Appendix C). To compute all entries of the matrix A̅, (3.44) and (3.45) are differentiated with respect to θ1 and f̂x and this yields(3.54)D̅A̅=F̅,
where matrices D̅ and F̅  are given by(3.55)D̅=[∂Γ̅o∂θ2∂Γ̅o∂f̂y∂Γ̅v∂θ2∂Γ̅v∂f̂y],F̅=-[∂Γ̅o∂θ1∂Γ̅o∂f̂x∂Γ̅v∂θ1∂Γ̅v∂f̂x].
Explicit expressions for all entries of the matrices D̅ and F̅ are reported in Appendix C. Once the matrix A̅ is solved from (3.54), it is substituted into (3.50) and (3.52) to obtain all entries of the matrix g̅p. Again, due to the complexity of all functions resulting from variable transformations, the matrices B̅,C̅,D̅, and F̅ are evaluated numerically using Gaussian quadrature.


## 3.3. Local Element Tangent Stiffness Matrix

Consider now a member with general boundary conditions as shown schematically in Figure3. Let {x,y} be a local coordinate system of the undeformed member and {x*,y*} be the coordinate system of the deformed member defined such that a chord connecting its end points always lies on the x* axis. With this specific choice of {x*,y*}, behavior of the member observed from this coordinate system is identical to that of a simply supported member.Figure 3
Schematic of undeformed and deformed configurations of a member subjected to general boundary conditions.The normalized end loads and normalized end displacements and rotations are denoted by{f̂x1,f̂y1,m̂1,f̂x2,f̂y2,m̂2} and {û1,v̂1,θ1,û2,v̂2,θ2} in the {x,y} coordinate system and by {f̂x1*,f̂y1*,m̂1*,f̂x2*,f̂y2*,m̂2*} and {û2*,θ1*,θ2*} in the {x*,y*} coordinate system. From geometric consideration of the deformed configuration, {û2*,θ1*,θ2*} can readily be expressed in terms of {û1,v̂1,θ1,û2,v̂2,θ2} by(3.56)θ1=θ1*+ϕ,θ2=θ2*+ϕ,û2*=(1+û2-û1)cos⁡⁡ϕ+(v̂2-v̂1)sin⁡⁡ϕ-1,
where ϕ is a chord rotation governed by(3.57)(1+û2-û1)sin⁡⁡ϕ-(v̂2-v̂1)cos⁡⁡ϕ=0.Letf̂x* be the internal force in x* direction and R* be the residual defined in the {x*,y*} system in the same way as given by (3.10). In the {x,y} coordinate system, we choose {f̂x,ℛ} such that(3.58)f̂x=f̂x*,R=R*.By applying the coordinate transformation, a vectorf={f̂x1,f̂y1,m̂1,f̂x2,f̂y2,m̂2,ℛ} is related to a vector f*={f̂x2*,m̂1*,m̂2*,ℛ*,f̂x1*,f̂y1*,f̂y2*} by(3.59)f=R(ϕ)f*,
where R(ϕ) is a transformation matrix of dimensions 7×7 given by(3.60)R(ϕ)=[0000cϕ-sϕ00000sϕcϕ00100000cϕ00000-sϕsϕ00000cϕ00100000001000],
in which sϕ=sinϕ and cϕ=cosϕ. By defining u={û1,v̂1,θ1,û2,v̂2,θ2,f̂x} and u*={û2*,θ1*,θ2*,f̂x*} and then recalling (3.56)–(3.58), we obtain the relation u*=u*(u). From the fact that behavior of the member in the {x*,y*} system is identical to that for a simply supported member, f* and u* are therefore related by f*=f*(u*). Combining (3.59), u*=u*(u) and f*=f*(u*) leads to the relation f=f(u)=R(ϕ)f*(u*(u)), and, from Taylor series expansion, this nonlinear function possesses the best linear approximation in the neighborhood of a vector u0(3.61)f(u)=f(u0)+kl(u0)(u-u0),
where kl is a local element tangent stiffness matrix given by(3.62)kl=∂R∂ϕf*∂ϕ∂u+Rg∂u*∂u,
in which the relation g=∂f*/∂u* is utilized.For a special case in which the member contains an inflection point at its right end, the end momentm̂2 vanishes and the corresponding end rotation θ2 is eliminated from the set of unknowns. Let us first define following reduced vectors f̅={f̂x1,f̂y1,m̂1,f̂x2,f̂y2,ℛ}, f̅*={f̂x2*,m̂1*,ℛ*,f̂x1*,f̂y1*,f̂y2*}, u̅={û1,v̂1,θ1,û2,v̂2,f̂x}, and u̅*={û2*,θ1*,f̂x*}. By applying the coordinate transformation, the relationship between the reduced vectors f̅ and f̅* is given by(3.63)f̅=R̅(ϕ)f̅*,
where R̅(ϕ) is a reduced transformation matrix of dimensions 6×6 given by(3.64)R̅(ϕ)=[000cϕ-sϕ0000sϕcϕ0010000cϕ0000-sϕsϕ0000cϕ001000].
The relation u̅*=u̅*(u̅) is then obtained by recalling (3.56)–(3.58) and f̅*=f̅*(u̅*) results from the fact that the behavior of the member in the {x*,y*} system is identical to that for a simply supported member. Combining (3.63), u̅*=u̅*(u̅) and f̅*=f̅*(u̅*) yields f̅=R̅(ϕ)f̅*(u̅*(u̅)), and, from Taylor series expansion, the best linear approximation of this nonlinear function in the neighborhood of a vector u̅0 takes the form(3.65)f̅(u̅)=f̅(u̅0)+k̅l(u̅0)(u̅-u̅0),
where k̅l is a reduced local element tangent stiffness matrix given by(3.66)k̅l=∂R̅∂ϕf̅*∂ϕ∂u̅+R̅g̅∂u̅*∂u̅,
in which the relation g̅=∂f̅*/∂u̅* is utilized. Note that the local element tangent stiffness matrix for this particular case is of dimensions 6×6.


## 3.4. Global Element Tangent Stiffness Matrix

Let the orientation of a member in the undeformed state be denoted by an angleβ measured from the global X-axis to the local x-axis (defined in the Section 3.3). The element tangent stiffness matrix kg referring to the global coordinate system {X,Y} can be obtained using the following standard coordinate transformation(3.67)kg=QTklQ,
where Q is a transformation matrix of dimension 7×7 given by(3.68)Q=[cβsβ00000-sβcβ000000010000000cβsβ00000-sβcβ0000000100000001],
in which sβ=sinβ and cβ=cosβ. Similarly, for a member containing an inflection point at the right end, we obtain(3.69)k̅g=Q̅Tk̅lQ̅,
where k̅g is the reduced, global element tangent stiffness matrix and Q̅ is a reduced transformation matrix given by(3.70)Q̅=[cβsβ0000-sβcβ0000001000000cβsβ0000-sβcβ0000001].


## 4. Structure Stiffness Equations

The best linear approximation of nonlinear algebraic equations governing the entire structure can readily be obtained by a direct assembly of the linear approximation of all members. This strategy employs two key ingredients, that is, the compatibility of the displacement and rotation at nodes and ends of members and the equilibrium between external loads and member end forces at nodes. The resulting linear approximation is given by(4.1)P=Po+Kt(U-Uo),
where P is a vector of nodal loads and zero residuals of all members, U is a vector of nodal displacements, nodal rotations and the internal axial force of all members, {Po,Uo} are vectors {P,U} at the reference state, and  Kt is the structure tangent stiffness matrix. Note that the matrix Kt can readily be obtained from a direct assembly of the global element tangent stiffness matrices kg or k̅g of all members. The linear approximation (4.1) is then used in the Newton-Raphson iteration to search for a solution of a system of nonlinear equations governing the entire structure.


## 5. Verifications and Results

As a means to verify both the formulation and numerical implementation and also to demonstrate the capability and versatility of the proposed technique, extensive numerical experiments are performed for various structures. In the verification procedure, a set of simple boundary value problems is considered first and obtained results are compared with existing analytical solutions, and, subsequently, more complex structures are analyzed and results are verified by those obtained from a reliable commercial FEM package, ANSYS.


### 5.1. Cantilever Beam Subjected to Concentrated Moments

Consider a cantilever beam of lengthL and flexural rigidity EI and subjected to two concentrated moments −M (negative sign simply indicating that the applied moment is in clockwise direction) and 1.5 M where the former is applied at the tip and the latter is applied at the mid-span as shown in Figure 4(a). It is clear from equilibrium that the bending moment within the left half of the beam is equal to 0.5 M whereas, in its right half, the bending moment is equal to –M. Three uniform meshes (consisting of 2, 4, and 8 identical members) employed in the analysis are illustrated in Figure 4(b). Responses of the beam are obtained for several levels of the applied moment, that is, m̂=ML/EI∈{1,2,3,4,5}.(a) Schematic of cantilever beam subjected to two concentrated moments and (b) three meshes adopted in the analysis.
(a)(b)Since the bending moment for the left and right halves of the beam is constant and the internal resultant forces identically vanish, the closed form solution for the rotation and the displacement can readily be obtained from the direct integration of the governing equations (3.2). The explicit solutions are given by(5.1)θ={0.5m̂ξ,0≤ξ≤0.5,m̂(0.75-ξ),0.5≤ξ≤1,û={2sin⁡⁡(0.5m̂ξ)m̂-ξ,0≤ξ≤0.5,sin⁡⁡(0.25m̂)+sin⁡⁡((0.75-ξ)m̂)m̂+ξ-1,0.5≤ξ≤1,v̂={2-2cos⁡⁡(0.5m̂ξ)m̂,0≤ξ≤0.5,2-3cos⁡⁡(0.25m̂)+cos⁡⁡((0.75-ξ)m̂)m̂,0.5≤ξ≤1.
The deflected shapes of the beam for different values of m̂ are shown in Figure 5. Numerical results obtained for all three meshes are reported only at nodal points and compared with the analytical solution given by (5.1). As evident from this set of results, numerical solutions exhibit excellent agreement with the benchmark solution. It is important to emphasize that the accuracy of numerical solutions obtained from the proposed technique is independent of the level of mesh refinement; use of three meshes in the analysis is mainly to verify the implementation of the direct stiffness strategy for structures consisting of multiple members. While results are reported only at nodes, responses at any point within the member can readily be obtained once the nodal quantities are determined.Figure 5
Deformed shapes of cantilever beam subjected to two concentrated moments.


### 5.2. Frame Subjected to Concentrated Moments

Next, consider a more complex boundary value problem associated with a frame consisting of a column and two overhanging beams. The column and the two beams are of the same lengthL and flexural rigidity EI, and the frame is subjected to three concentrated moments {M1,M2,M3} as shown schematically in Figure 6(a). In the analysis, we choose {M1,M2,M3} such that M1=M2=αEI/L and M3=-2.5αEI/L, and three meshes (consisting of 3, 6, and 12 members) are adopted as shown in Figure 6(b). Note again that use of different meshes in the experiments is merely to demonstrate capability of the technique to model structures consisting of multiple members.(a) Schematic of frame subjected to three moments and (b) three meshes adopted in the analysis.
(a)(b)This particular loading condition yields a constant bending moment within the two beams and the column and zero resultant forces over the entire structure. Similar to the previous case, the analytical solution for the rotation and displacement at any point within the structure can readily be obtained by directly solving the governing equations (3.2) for each beam and column along with the use of boundary conditions at the fixed base and the continuity conditions at the junction. The deflected shapes of the rigid frame for various values of loading parameter α are reported in Figure 7. Again, numerical results for the displacement at the nodal points obtained from all three meshes coincide with the analytical solutions, and, in addition, no dependence on the level of mesh refinement is observed. These experiments again confirm the validity of the current implementation.Figure 7
Deformed shapes of frame subjected to three concentrated moments.


### 5.3. Opened Square Frame Subjected to Pair of Opposite Forces

Consider an opened square frame subjected to a pair of opposite vertical loads as shown in Figure8(a). The frame consists of a horizontal member, two vertical members of the same length L and two overhanging members of length 0.5 L, and all members have the same flexural rigidity EI. As clearly demonstrated by the two previous problems, the current technique displays an attractive feature, namely, the independence of a level of mesh refinement. Thus, without loss of accuracy of numerical solutions, it is common to discretize the structure using the minimum number of elements to reduce the computational cost. For this particular case, a mesh consisting of 3 horizontal members and 2 vertical members is adopted as shown in Figure 8(b).(a) Schematic of opened square frame subjected to pair of opposite vertical forces and (b) mesh used in analysis.
(a)(b)The deflected shapes of the structure are reported in Figure9 for many values of normalized vertical loads f*=FL2/EI∈{0.2,0.6,1.0,1.4,1.8}. From this set of results, the numerical solutions exhibit excellent agreement with the benchmark solutions obtained from a reliable commercial FEM package, ANSYS; in fact, the computed results from the proposed technique are nearly indistinguishable from those obtained from ANSYS. It is also worth pointing out that in the construction of a (converged) benchmark solution by ANSYS, a series of meshes were considered and a reasonably fine mesh was needed to achieve such highly accurate results.Figure 9
Deformed shapes of opened square frame subjected to pair of opposite vertical forces.


### 5.4. Diamond Box Frame Subjected to Vertical Force

As a final example, consider a diamond box frame subjected to a vertical downward forceF as shown schematically in Figure 10(a). The frame consists of four 45-degree bevel elements of the same length 2L, and all elements have the same flexural rigidity EI. In the analysis, the structure is discretized into four elements (see Figure 10(b)) and various values of a normalized load (i.e., f̂y=FL2/EI∈{5,10,15,20}) are treated.(a) Schematic of diamond box frame subjected to vertical force at upper tip and (b) mesh adopted in analysis.
(a)(b)Deflected shapes of a structure for differentf̂y are reported in Figure 11. It can be concluded again from these results that the proposed technique yields highly accurate solutions for any level of applied loads treated. In particular, the computed displacements and rotations coincide with those obtained from ANSYS. In addition, the nonlinear relationship between the vertical applied load and the vertical displacement at the upper tip is also investigated and results are reported in Figure 12. Besides the obvious monotonic increase of the applied load versus the displacement, it is observed that the stiffness of the structure (represented by the slope of the load-displacement curve) gradually decreases, for small f̂y and then starts to increase rapidly for large f̂y. This is due to that for small f̂y, change of structure configuration is not significant and the axial force within all members is still in compression, and this results in a reduction of the member stiffness due to the influence of geometric nonlinearity. As f̂y is sufficiently large, the configuration of the structure changes considerably from the initial state (see Figure 11) and the axial force within the members switches from compression to tension and thus increasing the member stiffness.Figure 11
Deformed shapes of diamond box frame subjected to vertical force at upper tip.Figure 12
Relation between downward vertical displacement of upper tip and normalized vertical force.


## 5.1. Cantilever Beam Subjected to Concentrated Moments

Consider a cantilever beam of lengthL and flexural rigidity EI and subjected to two concentrated moments −M (negative sign simply indicating that the applied moment is in clockwise direction) and 1.5 M where the former is applied at the tip and the latter is applied at the mid-span as shown in Figure 4(a). It is clear from equilibrium that the bending moment within the left half of the beam is equal to 0.5 M whereas, in its right half, the bending moment is equal to –M. Three uniform meshes (consisting of 2, 4, and 8 identical members) employed in the analysis are illustrated in Figure 4(b). Responses of the beam are obtained for several levels of the applied moment, that is, m̂=ML/EI∈{1,2,3,4,5}.(a) Schematic of cantilever beam subjected to two concentrated moments and (b) three meshes adopted in the analysis.
(a)(b)Since the bending moment for the left and right halves of the beam is constant and the internal resultant forces identically vanish, the closed form solution for the rotation and the displacement can readily be obtained from the direct integration of the governing equations (3.2). The explicit solutions are given by(5.1)θ={0.5m̂ξ,0≤ξ≤0.5,m̂(0.75-ξ),0.5≤ξ≤1,û={2sin⁡⁡(0.5m̂ξ)m̂-ξ,0≤ξ≤0.5,sin⁡⁡(0.25m̂)+sin⁡⁡((0.75-ξ)m̂)m̂+ξ-1,0.5≤ξ≤1,v̂={2-2cos⁡⁡(0.5m̂ξ)m̂,0≤ξ≤0.5,2-3cos⁡⁡(0.25m̂)+cos⁡⁡((0.75-ξ)m̂)m̂,0.5≤ξ≤1.
The deflected shapes of the beam for different values of m̂ are shown in Figure 5. Numerical results obtained for all three meshes are reported only at nodal points and compared with the analytical solution given by (5.1). As evident from this set of results, numerical solutions exhibit excellent agreement with the benchmark solution. It is important to emphasize that the accuracy of numerical solutions obtained from the proposed technique is independent of the level of mesh refinement; use of three meshes in the analysis is mainly to verify the implementation of the direct stiffness strategy for structures consisting of multiple members. While results are reported only at nodes, responses at any point within the member can readily be obtained once the nodal quantities are determined.Figure 5
Deformed shapes of cantilever beam subjected to two concentrated moments.


## 5.2. Frame Subjected to Concentrated Moments

Next, consider a more complex boundary value problem associated with a frame consisting of a column and two overhanging beams. The column and the two beams are of the same lengthL and flexural rigidity EI, and the frame is subjected to three concentrated moments {M1,M2,M3} as shown schematically in Figure 6(a). In the analysis, we choose {M1,M2,M3} such that M1=M2=αEI/L and M3=-2.5αEI/L, and three meshes (consisting of 3, 6, and 12 members) are adopted as shown in Figure 6(b). Note again that use of different meshes in the experiments is merely to demonstrate capability of the technique to model structures consisting of multiple members.(a) Schematic of frame subjected to three moments and (b) three meshes adopted in the analysis.
(a)(b)This particular loading condition yields a constant bending moment within the two beams and the column and zero resultant forces over the entire structure. Similar to the previous case, the analytical solution for the rotation and displacement at any point within the structure can readily be obtained by directly solving the governing equations (3.2) for each beam and column along with the use of boundary conditions at the fixed base and the continuity conditions at the junction. The deflected shapes of the rigid frame for various values of loading parameter α are reported in Figure 7. Again, numerical results for the displacement at the nodal points obtained from all three meshes coincide with the analytical solutions, and, in addition, no dependence on the level of mesh refinement is observed. These experiments again confirm the validity of the current implementation.Figure 7
Deformed shapes of frame subjected to three concentrated moments.


## 5.3. Opened Square Frame Subjected to Pair of Opposite Forces

Consider an opened square frame subjected to a pair of opposite vertical loads as shown in Figure8(a). The frame consists of a horizontal member, two vertical members of the same length L and two overhanging members of length 0.5 L, and all members have the same flexural rigidity EI. As clearly demonstrated by the two previous problems, the current technique displays an attractive feature, namely, the independence of a level of mesh refinement. Thus, without loss of accuracy of numerical solutions, it is common to discretize the structure using the minimum number of elements to reduce the computational cost. For this particular case, a mesh consisting of 3 horizontal members and 2 vertical members is adopted as shown in Figure 8(b).(a) Schematic of opened square frame subjected to pair of opposite vertical forces and (b) mesh used in analysis.
(a)(b)The deflected shapes of the structure are reported in Figure9 for many values of normalized vertical loads f*=FL2/EI∈{0.2,0.6,1.0,1.4,1.8}. From this set of results, the numerical solutions exhibit excellent agreement with the benchmark solutions obtained from a reliable commercial FEM package, ANSYS; in fact, the computed results from the proposed technique are nearly indistinguishable from those obtained from ANSYS. It is also worth pointing out that in the construction of a (converged) benchmark solution by ANSYS, a series of meshes were considered and a reasonably fine mesh was needed to achieve such highly accurate results.Figure 9
Deformed shapes of opened square frame subjected to pair of opposite vertical forces.


## 5.4. Diamond Box Frame Subjected to Vertical Force

As a final example, consider a diamond box frame subjected to a vertical downward forceF as shown schematically in Figure 10(a). The frame consists of four 45-degree bevel elements of the same length 2L, and all elements have the same flexural rigidity EI. In the analysis, the structure is discretized into four elements (see Figure 10(b)) and various values of a normalized load (i.e., f̂y=FL2/EI∈{5,10,15,20}) are treated.(a) Schematic of diamond box frame subjected to vertical force at upper tip and (b) mesh adopted in analysis.
(a)(b)Deflected shapes of a structure for differentf̂y are reported in Figure 11. It can be concluded again from these results that the proposed technique yields highly accurate solutions for any level of applied loads treated. In particular, the computed displacements and rotations coincide with those obtained from ANSYS. In addition, the nonlinear relationship between the vertical applied load and the vertical displacement at the upper tip is also investigated and results are reported in Figure 12. Besides the obvious monotonic increase of the applied load versus the displacement, it is observed that the stiffness of the structure (represented by the slope of the load-displacement curve) gradually decreases, for small f̂y and then starts to increase rapidly for large f̂y. This is due to that for small f̂y, change of structure configuration is not significant and the axial force within all members is still in compression, and this results in a reduction of the member stiffness due to the influence of geometric nonlinearity. As f̂y is sufficiently large, the configuration of the structure changes considerably from the initial state (see Figure 11) and the axial force within the members switches from compression to tension and thus increasing the member stiffness.Figure 11
Deformed shapes of diamond box frame subjected to vertical force at upper tip.Figure 12
Relation between downward vertical displacement of upper tip and normalized vertical force.


## 6. Conclusion

A simple, systematic method has been developed for analysis of flexure-dominating skeleton structures undergoing large displacement and rotation. The technique is based upon the direct stiffness method along with the use of exact element tangent stiffness matrices. These matrices have been derived at a member level using a classical elastica approach and the constraint condition posed by member inextensibility has been directly incorporated in such development. This therefore results in the element tangent stiffness matrix of dimension7×7. It is also worth noting that the resulting tangent stiffness matrix possesses two positive features: (i) it is exact in the sense that it involves no approximation of both the solution form and the governing equations, and (ii) all entries of the matrix are given in an explicit form concerning the elliptic or other similar integrals. The former feature enhances the rate of convergence of a nonlinear solver, and, when properly incorporated with the evaluation of exact residuals, it can in principle yield numerical solutions of the same quality as an analytical solution. The latter feature is well-suited for numerical evaluation of the tangent stiffness matrix by a standard Gaussian quadrature.As evident from various numerical experiments, the current technique offers two crucial benefits. Firstly, the method provides a simple and systematic means to model large structures of various geometries (consisting of multiple members with different orientations) by using exact kinematics (i.e., exact curvature-displacement relationship), and, secondly, it provides “exact” numerical solutions (within round off errors and errors caused by a nonlinear solver and numerical quadrature) that are independent of mesh refinements. One practical contribution of the current investigation is that it provides an accurate computational tool well suited for analysis of structures undergoing large displacement and rotation, for example, very flexible structures, moment-resisting cables, slender drill string rods, and so forth. According to its high accuracy, the proposed technique can also be employed to generate benchmark solutions for a comparison purpose. As a final remark, the proposed technique can further be generalized to treat two important classes of problems, one associated with the treatment of material nonlinearity and the other corresponding to nonlinear structural dynamics. It is apparent that, under severe time-dependent loading conditions (e.g., earthquake and blasting loads), not only the inertia effect becomes significant but also structures generally undergo both large displacement and large deformation prior to collapse.



---

*Source: 592958-2011-12-19.xml*