Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2011 / Article

Research Article | Open Access

Volume 2011 |Article ID 592958 | https://doi.org/10.1155/2011/592958

Jaroon Rungamornrat, Peerasak Tangnovarad, "Analysis of Linearly Elastic Inextensible Frames Undergoing Large Displacement and Rotation", Mathematical Problems in Engineering, vol. 2011, Article ID 592958, 37 pages, 2011. https://doi.org/10.1155/2011/592958

Analysis of Linearly Elastic Inextensible Frames Undergoing Large Displacement and Rotation

Academic Editor: Mohammad Younis
Received18 Apr 2011
Accepted09 Aug 2011
Published19 Dec 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.

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., [37]). 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 [1021] 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 of arbitrary 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., [2832]) 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 length 𝐿 and moment of inertia 𝐼 and made of an elastic material with Young’s modulus 𝐸. An undeformed configuration of this member occupies a straight line 𝑦=0,0𝑥𝐿, 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 (𝑆,0)” and “point (𝑆,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 (𝑆,0) in the undeformed state. During a particular motion, the cross-section at (𝑆,0) in the undeformed state displaces to a point (𝑆+𝑢,𝑣) in the deformed state where 𝑢=𝑢(𝑆) and 𝑣=𝑣(𝑆) denote the 𝑥-component and the 𝑦-component of the displacement at point (𝑆,0), respectively. Let 𝑓𝑥=𝑓𝑥(𝑆),𝑓𝑦=𝑓𝑦(𝑆), and 𝑚=𝑚(𝑆) denote a resultant internal force in 𝑥-direction, a resultant internal force in 𝑦-direction, and a resultant bending moment at the cross-section (𝑆,0), respectively.

By enforcing static equilibrium of an infinitesimal element 𝑑𝑆 in the deformed state (see its free body diagram in Figure 1(b)), it leads to three differential equations𝑑𝑓𝑥𝑑𝑆=0,(2.1)𝑑𝑓𝑦𝑑𝑆=0,(2.2)𝑑𝑚𝑑𝑆=𝑓𝑥sin𝜃+𝑓𝑦cos𝜃,(2.3) where 𝜃 is the rotation of any cross-section. Clearly indicated by (2.1) and (2.2), the internal forces 𝑓𝑥 and 𝑓𝑦 must be constant for the entire member. It is evident from geometry of the element 𝑑𝑆 in the deformed state that the rotation 𝜃 can be related to the two components of the displacement 𝑢 and 𝑣 bysin𝜃=𝑑𝑣𝑑𝑆,cos𝜃=1+𝑑𝑢𝑑𝑆.(2.4) From the kinematics assumption (ii), the curvature 𝜅 and the rotation 𝜃 are related through 𝜅=𝑑𝜃/𝑑𝑆, and, by using assumptions (i) and (vi), we then obtain the linear moment-curvature relationship: 𝑚=EI𝜅. By using these two relations, (2.3) becomes𝑑2𝜃𝑑𝜉2=𝑓𝑥𝑓sin𝜃+𝑦cos𝜃,(2.5) where nondimensional parameters appearing in above equation are defined by 𝜉=𝑆/𝐿, 𝑓𝑥=𝑓𝑥𝐿2/EI, and 𝑓𝑦=𝑓𝑦𝐿2/EI. By performing a direct integration of (2.5), it leads to 𝑑𝜃𝑑𝜉2𝑓=𝐶2𝑥𝑓cos𝜃+2𝑦sin𝜃,(2.6) where 𝐶 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 𝑑𝜃/𝑑𝜉 is uniquely dictated by that of the bending moment 𝑚. As a consequence, 𝑑𝜉/𝑑𝜃 can readily be solved from (2.6) to obtain𝑑𝜉=𝜗𝑑𝜃𝑚𝑓𝐶2𝑥𝑓cos𝜃+2𝑦sin𝜃,(2.7) where 𝜗(𝑚) is a moment-dependent function defined such that 𝜗(𝑚)=1 if 𝑚=𝑚𝐿/EI>0 and 𝜗(𝑚)=1 if 𝑚<0. By combining (2.4) and (2.7), it leads to following two relations:𝑑̂𝑣=𝜗𝑑𝜃𝑚sin𝜃𝑓𝐶2𝑥𝑓cos𝜃+2𝑦,sin𝜃𝑑̂𝑢=𝜗𝑑𝜃𝑚(cos𝜃1)𝑓𝐶2𝑥𝑓cos𝜃+2𝑦,sin𝜃(2.8) where ̂𝑢=𝑢/𝐿 and ̂𝑣=𝑣/𝐿. 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 (𝑆=0) and a roller support at its right end (𝑆=𝐿)) subjected to end loads {𝑚1,𝑚2,𝑓𝑥2} as shown in Figure 2. By imposing the moment boundary condition at the right end, that is, 𝑑𝜃/𝑑𝜉(𝜃2)=𝑚2=𝑚2𝐿/EI, we obtain the constant 𝐶 as𝐶=𝑚22𝑓+2𝑥cos𝜃2𝑓2𝑦sin𝜃2.(3.1) By inserting the constant 𝐶 into the relations (2.7)-(2.8), it yields𝑑𝜉𝐹𝑑𝜃=𝜗𝑚𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2,𝑑̂𝑣𝑑𝜃=𝜗𝑚sin𝜃𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2,𝑑̂𝑢𝑑𝜃=𝜗𝑚(cos𝜃1)𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2,(3.2) where the function 𝐹 is defined by𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2=1𝑚22𝑓+2𝑥cos𝜃2𝑓cos𝜃2𝑦sin𝜃2sin𝜃.(3.3) By imposing the remaining two natural boundary conditions at the left and right ends of a member, we obtain𝑚22𝑚21𝑓+2𝑥cos𝜃2cos𝜃1𝑓2𝑦sin𝜃2sin𝜃1𝑓=0,(3.4)𝑥2𝑓𝑥=0,(3.5) where 𝑚1=𝑚1𝐿/EI and 𝑓𝑥2=𝑓𝑥2𝐿2/EI. Normalized support reactions {𝑓𝑥1,𝑓𝑦1,𝑓𝑦2} can readily be computed from equilibrium of the entire member, and final results are given by𝑓𝑥1𝑓=𝑥,𝑓𝑦1=𝑚1+𝑚2𝑑,𝑓𝑦2=𝑚1+𝑚2𝑑,(3.6) where 𝑑=1+̂𝑢2 and ̂𝑢2=𝑢2/𝐿. Next, by integrating (3.2) from 𝜃=𝜃1 to 𝜃=𝜃2, it leads to a system of three nonlinear algebraic equations 𝜃2𝜃1𝜗𝐹𝑚𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑑𝜃=1,(3.7)𝜃2𝜃1𝜗𝑚sin𝜃𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑑𝜃=0,(3.8)𝜃2𝜃1𝜗𝑚cos𝜃𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑑𝑑𝜃=.(3.9) For a given set of end loads {𝑓𝑥2,𝑚1,𝑚2}, the end displacement and rotations {̂𝑢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 {𝑓𝑥,𝑓𝑦}. On the other hand, the problem finding the end loads {𝑓𝑥2,𝑚1,𝑚2} for a particular set of prescribed displacement and rotations {̂𝑢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, {̂𝑢2,𝜃1,𝜃2} cannot be specified arbitrarily. However, if {𝑓𝑥,𝜃1,𝜃2} are prescribed instead, the end loads {𝑓𝑥2,𝑚1,𝑚2} can be obtained by solving (3.4), (3.5), (3.7), and (3.8) simultaneously, and the end displacement ̂𝑢2 can subsequently be computed from (3.9). However, lack of the displacement component ̂𝑢2 renders a set {𝑓𝑥,𝜃1,𝜃2} not well suited for the development of a solution procedure by a direct stiffness method.

In the present study, {̂𝑢2,𝜃1,𝜃2,𝑓𝑥} is chosen as a set of primary unknowns. It is worth noting that to allow ̂𝑢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𝓡𝑑𝜃2𝜃1𝜗𝑚cos𝜃𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑑𝜃.(3.10) Supplemented by (3.10), for any given set {̂𝑢2,𝜃1,𝜃2,𝑓𝑥}, the quantities {𝑓𝑥2,𝑚1,𝑚2,𝓡} can be obtained from (3.5), (3.7), (3.8), and (3.10) with use of (3.4) to get rid of 𝑓𝑦. It should be emphasized that {̂𝑢2,𝜃1,𝜃2,𝑓𝑥} and the corresponding {𝑓𝑥2,𝑚1,𝑚2,𝓡} are solutions of the boundary value problem if and only if the residual vanishes.

Let 𝐟 be a vector defined by 𝐟=[𝐟𝑝𝐟𝑟]𝑇 where 𝐟𝑝𝑓={𝑥2,𝑚1,𝑚2,𝓡} and 𝐟𝑟𝑓={𝑥1,𝑓𝑦1,𝑓𝑦2} and let 𝐮 be a vector defined by 𝐮={̂𝑢2,𝜃1,𝜃2,𝑓𝑥}𝑇. From the relations (3.4)–(3.8) and (3.10), it can readily be verified that 𝐟=𝐟(𝐮), and, from Taylor series expansion, the best linear approximation of this nonlinear function 𝐟 in the neighborhood of a vector 𝐮0 takes a form𝐟𝐮(𝐮)=𝐟0𝐮+𝐠0𝐮𝐮0,(3.11) where a gradient matrix 𝐠 is given by𝐠=𝜕𝐟=𝐠𝜕𝐮𝑝𝐠𝑟.(3.12) The submatrix 𝐠𝑝 is expressed explicitly in terms of differentiations as𝐠𝑝=𝜕𝑓𝑥2𝜕̂𝑢2𝜕𝑓𝑥2𝜕𝜃1𝜕𝑓𝑥2𝜕𝜃2𝜕𝑓𝑥2𝜕𝑓𝑥𝜕𝑚1𝜕̂𝑢2𝜕𝑚1𝜕𝜃1𝜕𝑚1𝜕𝜃2𝜕𝑚1𝜕𝑓𝑥𝜕𝑚2𝜕̂𝑢2𝜕𝑚2𝜕𝜃1𝜕𝑚2𝜕𝜃2𝜕𝑚2𝜕𝑓𝑥𝜕𝓡𝜕̂𝑢2𝜕𝓡𝜕𝜃1𝜕𝓡𝜕𝜃2𝜕𝓡𝜕𝑓𝑥.(3.13) By defining 𝑔𝑖𝑗 as an entry located at the 𝑖th row and 𝑗th column of the submatrix 𝐠𝑝, the submatrix 𝐠𝑟 can readily be obtained in terms of 𝑔𝑖𝑗 as𝐠𝑟=1𝑑𝑔11𝑑𝑔12𝑑𝑔13𝑑𝑔14𝑑̂𝑠𝑔22+𝑔32𝑔23+𝑔33𝑔24+𝑔34̂𝑠𝑔22𝑔32𝑔23𝑔33𝑔24𝑔34,(3.14) where ̂𝑠=(𝑚1+𝑚2𝑑)/. As clearly indicated by (3.5), (3.7), (3.8), and (3.10), entries 𝑔11,𝑔12,𝑔13,𝑔21, and 𝑔31 vanish whereas 𝑔14=𝑔41=1. The remaining entries of 𝐠𝑝 are contained in a submatrix 𝓰𝑝 given by𝓰𝑝=𝜕𝑚1𝜕𝜃1𝜕𝑚1𝜕𝜃2𝜕𝑚1𝜕𝑓𝑥𝜕𝑚2𝜕𝜃1𝜕𝑚2𝜕𝜃2𝜕𝑚2𝜕𝑓𝑥𝜕𝓡𝜕𝜃1𝜕𝓡𝜕𝜃2𝜕𝓡𝜕𝑓𝑥.(3.15) It should be remarked that determination of all entries of the submatrix 𝓰𝑝 is nontrivial and requires implicit differentiations.

3.2. Derivation of Submatrix 𝑝

The explicit form of the submatrix 𝓰𝑝 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 𝜗(𝑚) and the singularity behavior of the function 𝐹 defined by (3.3).

3.2.1. Member Containing No Inflection Point

The normalized bending moment is strictly positive (i.e., 𝑚(𝜉)>0) or strictly negative (i.e., 𝑚(𝜉)<0) for the entire beam (i.e., 𝜉[0,1]) when the applied end moments {𝑚1,𝑚2} are nonzero and possess different sign (i.e., 𝑚1𝑚2<0). The deformed configuration of the member for this particular case possesses a single curvature, and, in addition, 𝜗(𝑚) becomes a constant function with its value equal to either 1 or −1 depending on the sign of 𝑚; more specifically, 𝜗(𝑚)=1for 𝑚>0 and 𝜗(𝑚)=1 for 𝑚<0. It is worth noting that the function 𝐹, 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 𝑚0 for the entire member. Such desirable feature of 𝐹 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Γ1𝜃1,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2=𝜃2𝜃1𝜗𝐹𝑚𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2Γ𝑑𝜃1=0,(3.16)2𝜃1,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2=𝜃2𝜃1𝜗𝑚sin𝜃𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝓡𝜃𝑑𝜃=0,(3.17)𝓡=1,𝜃2,̂𝑢2;𝑓𝑥,𝑓𝑦,𝑚2=𝑑𝜃2𝜃1𝜗𝑚cos𝜃𝐹𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑑𝜃.(3.18) It is noted that the relation (3.4) implicitly defines 𝑓𝑦 in terms of {𝑚1,𝑚2} and {𝜃1,𝜃2,𝑓𝑥}. As a result, by taking derivative of (3.16)–(3.18) with respect to {𝜃1,𝜃2,𝑓𝑥} along with employing the chain rule of differentiation, it leads to𝐒𝓰𝑝=𝐃,(3.19) where the matrices 𝐒 and 𝐃 are given by𝐒=𝜕Γ1𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑚1𝜕Γ1𝜕𝑚2+𝜕Γ1𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑚20𝜕Γ2𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑚1𝜕Γ2𝜕𝑚2+𝜕Γ2𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑚20𝜕𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑚1𝜕𝜕𝑚2𝜕𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑚21,(3.20)𝐃=𝜕Γ1𝜕𝜃1+𝜕Γ1𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝜃1𝜕Γ1𝜕𝜃2+𝜕Γ1𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝜃2𝜕Γ1𝜕𝑓𝑥+𝜕Γ1𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑓𝑥𝜕Γ2𝜕𝜃1+𝜕Γ2𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝜃1𝜕Γ2𝜕𝜃2+𝜕Γ2𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝜃2𝜕Γ2𝜕𝑓𝑥+𝜕Γ2𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑓𝑥𝜕𝜕𝜃1𝜕𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝜃1𝜕𝜕𝜃2𝜕𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝜃2𝜕𝜕𝑓𝑥𝜕𝜕𝑓𝑦𝜕𝑓𝑦𝜕𝑓𝑥.(3.21) Entries of the matrices 𝐒 and 𝐃 can be expressed in a closed form and the submatrix 𝑝 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 𝜉𝑧(0,1) or, equivalently, the bending moment vanishing at 𝜉𝑧 and 𝑚(𝜉1)𝑚(𝜉2)<0 for 𝜉1[0,𝜉𝑧) and 𝜉2(𝜉𝑧,1]. This particular case arises when the applied end moments {𝑚1,𝑚2} are nonzero and possess the same sign (i.e., 𝑚1𝑚2>0). The resulting deformed configuration of the member possesses a double-curvature shape. As a result, the moment-dependent function 𝜗 is discontinuous at 𝜉𝑧 and takes different values on both sides of the inflection point. For the applied end moments 𝑚1,𝑚2>0, it results in 𝜗=1 for 𝜉[0,𝜉𝑧) and 𝜗=1for 𝜉(𝜉𝑧,1], and for 𝑚1,𝑚2<0, it results in 𝜗=1 for 𝜉[0,𝜉𝑧) and 𝜗=1 for 𝜉(𝜉𝑧,1].

For this special case, the constant 𝐶 appearing in (2.6) is alternatively obtained by imposing a condition at the inflection point, that is, 𝑑𝜃/𝑑𝜉(𝜃𝑧)=0, and this leads to𝑓𝐶=2𝑥cos𝜃𝑧𝑓2𝑦sin𝜃𝑧,(3.22) where 𝜃𝑧 is the rotation at the inflection point. Upon using (3.22), the relations (2.7)-(2.8) become𝑑𝜉𝐹𝑑𝜃=𝜗𝑚𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦,𝑑̂𝑣𝑑𝜃=𝜗𝑚sin𝜃𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦,𝑑̂𝑢𝑑𝜃=𝜗𝑚(cos𝜃1)𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦,(3.23) where the function 𝐹𝑧 is defined by𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦=12𝑓𝑥cos𝜃𝑧𝑓cos𝜃2𝑦sin𝜃𝑧sin𝜃.(3.24) By integrating equations (3.23) over the entire member, we obtain𝜓𝜃𝑧𝜃1𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦𝑑𝜃+𝜃2𝜃𝑧𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦𝜓𝑑𝜃=1,𝜃𝑧𝜃1sin𝜃𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦𝑑𝜃+𝜃2𝜃𝑧sin𝜃𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦𝜓𝑑𝜃=0,𝜃𝑧𝜃1(cos𝜃1)𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦𝑑𝜃+𝜃2𝜃𝑧(cos𝜃1)𝐹𝑧𝜃,𝜃𝑧;𝑓𝑥,𝑓𝑦𝑑𝜃=̂𝑢2,(3.25) where a constant 𝜓 is defined such that 𝜓=1 if 𝑚1,𝑚2>0 and 𝜓=1 if 𝑚1,𝑚2<0.

It is evident from (3.24) that the function 𝐹𝑧 is weakly singular at the inflection point (i.e., at 𝜃=𝜃𝑧); 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 𝑓2𝑠=𝑓2𝑥+𝑓2𝑦,cos𝜃o=𝑓𝑥/𝑓2𝑠, and sin𝜃o𝑓=𝑦/𝑓2𝑠, (ii) applying change of variable 𝜃=𝜋+(𝜃𝜃𝑜) and identity cos𝜃=12sin2(𝜃/2), and (iii) introducing another variable transformation 𝑝sin𝜙=sin(𝜃/2) with 𝑝=sin(𝜃𝑧/2). The function 𝐹𝑧 simply reduces to 𝐹𝑧=12𝑓2𝑠𝜃cos𝑧𝜃𝑜cos𝜃𝜃𝑜=14𝑓2𝑠sin2𝜃𝑧/2sin2=1𝜃/22𝑓𝑠𝑝cos𝜙,(3.26) where 𝜃𝑧=𝜋+(𝜃𝑧𝜃𝑜). Finally, the nonlinear algebraic equations (3.25) can be expressed in terms of integrals over a dummy variable 𝜙 asΓ𝑜𝜃1,𝜃2,𝜃𝑧,𝑓𝑥,𝑓𝑦=𝜙𝜋/21𝑓𝑜𝜙;𝑝𝑑𝜙+𝜙𝜋/22𝑓𝑜𝜙;𝑝𝑓𝑑𝜙𝑠Γ=0,(3.27)𝑣𝜃1,𝜃2,𝜃𝑧,𝑓𝑥,𝑓𝑦=𝜙𝜋/21𝑓𝑣𝜙,𝜃𝑜;𝑝𝑑𝜙+𝜙𝜋/22𝑓𝑣𝜙,𝜃𝑜;𝑝𝜃𝑑𝜙=0,(3.28)1,𝜃2,̂𝑢2,𝑓𝑥=𝜃1,𝜃2,𝜃𝑧,̂𝑢2,𝑓𝑥,𝑓𝑦=̂𝑢21𝑓𝑠𝜙𝜋/21𝑓𝑢𝜙,𝜃𝑜;𝑝𝑑𝜙+𝜙𝜋/22𝑓𝑢𝜙,𝜃𝑜;𝑝,𝑑𝜙(3.29) where 𝑝sin𝜙1=sin(𝜃1/2), 𝑝sin𝜙2=sin(𝜃2/2), 𝜃1=𝜋+(𝜃1𝜃o), 𝜃2=𝜋+(𝜃2𝜃o), and 𝑓𝑜𝜙;𝑝=11𝑝2sin2𝜙,𝑓𝑣𝜙,𝜃𝑜;𝑝=sin𝜃𝑜1𝑝2sin2𝜙2sin𝜃𝑜1𝑝2sin2𝜙2𝜓𝑝cos𝜃𝑜𝑓sin𝜙,𝑢𝜙,𝜃𝑜;𝑝=cos𝜃𝑜11𝑝2sin2𝜙2cos𝜃𝑜1𝑝2sin2𝜙+2𝜓𝑝sin𝜃𝑜sin𝜙.(3.30)

By enforcing the remaining two moment boundary conditions at both ends of the member, it leads to𝑚21𝑓+2𝑥cos𝜃1cos𝜃𝑧𝑓2𝑦sin𝜃1sin𝜃𝑧=0,𝑚22𝑓+2𝑥cos𝜃2cos𝜃𝑧𝑓2𝑦sin𝜃2sin𝜃𝑧=0.(3.31) By taking derivative of (3.31) with respect to 𝜃1, 𝜃2, and 𝑓𝑥, it results in𝜕𝑚1𝜕𝜃1𝜕𝑚1𝜕𝜃2𝜕𝑚1𝜕𝑓𝑥𝜕𝑚2𝜕𝜃1𝜕𝑚2𝜕𝜃2𝜕𝑚2𝜕𝑓𝑥=𝐌1+𝐌2𝐀,(3.32) where matrices 𝐌1, 𝐌2, and 𝐀 are given by𝐌1=𝜆1𝑚10𝑐𝑧𝑐1𝑚10𝜆2𝑚2𝑐𝑧𝑐2𝑚2𝐌,(3.33)2=𝜆𝑧𝑚1𝑠1𝑠𝑧𝑚1𝜆𝑧𝑚2𝑠2𝑠𝑧𝑚2,(3.34)𝐀=𝜕𝜃𝑧𝜕𝜃1𝜕𝜃𝑧𝜕𝜃2𝜕𝜃𝑧𝜕𝑓𝑥𝜕𝑓𝑦𝜕𝜃1𝜕𝑓𝑦𝜕𝜃2𝜕𝑓𝑦𝜕𝑓𝑥,(3.35) where 𝑠1=sin𝜃1, 𝑠2=sin𝜃2, 𝑠𝑧=sin𝜃𝑧, 𝑐1=cos𝜃1, 𝑐2=cos𝜃2, 𝑐𝑧=cos𝜃𝑧, 𝜆1=𝑓𝑥𝑠1+𝑓𝑦𝑐1, 𝜆2=𝑓𝑥𝑠2+𝑓𝑦𝑐2 and 𝜆𝑧=𝑓𝑥𝑠𝑧+𝑓𝑦𝑐𝑧. Differentiating (3.29) with respect to 𝜃1,𝜃2, and 𝑓𝑥 yields𝜕𝜕𝜃1𝜕𝜕𝜃2𝜕𝜕𝑓𝑥=𝐁+𝐂𝐀,(3.36) where matrices 𝐁 and 𝐂 are given by𝜕𝐁=𝜕𝜃1𝜕𝜕𝜃2𝜕𝜕𝑓𝑥𝜕,𝐂=𝜕𝜃𝑧𝜕𝜕𝑓𝑦.(3.37) Note that all entries of the matrices B and C can be obtained directly from (3.29) along with the change of variables 𝜃=𝜋+(𝜃𝜃𝑜) and 𝑝sin𝜙=sin(𝜃/2) (see explicit results in Appendix B).

To compute all entries of the matrix 𝐀, we differentiate (3.27) and (3.28) with respect to 𝜃1,𝜃2, and 𝑓𝑥, and this results in𝐃𝐀=𝐅,(3.38) where matrices 𝐃 and 𝐅 are given by𝐃=𝜕Γ𝑜𝜕𝜃𝑧𝜕Γ𝑜𝜕𝑓𝑦𝜕Γ𝑣𝜕𝜃𝑧𝜕Γ𝑣𝜕𝑓𝑦,𝐅=𝜕Γ𝑜𝜕𝜃1𝜕Γ𝑜𝜕𝜃2𝜕Γ𝑜𝜕𝑓𝑥𝜕Γ𝑣𝜕𝜃1𝜕Γ𝑣𝜕𝜃2𝜕Γ𝑣𝜕𝑓𝑥.(3.39) Similarly, all entries of the matrices 𝐃 and 𝐅 can be obtained directly from (3.27) and (3.28) along with the change of variables 𝜃=𝜋+(𝜃𝜃𝑜) and 𝑝sin𝜙=sin(𝜃/2) (see explicit results in Appendix B). Once the matrix 𝐀 is solved from (3.38), it is inserted into (3.32) and (3.36) to obtain all entries of the matrix 𝓰𝑝. Due to the complexity of all functions resulting from variable transformations, the matrices 𝐁,𝐂,𝐃, and 𝐅 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 {𝑚1,𝑚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 nonzero 𝑚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 𝜃𝑧=𝜃2 and 𝜉𝑧=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 𝜃𝑧=𝜃2 into (3.25), it leads to𝜓𝜃2𝜃1𝐹𝑧1𝜃,𝜃2;𝑓𝑥,𝑓𝑦𝑑𝜃=1,(3.40)𝜓𝜃2𝜃1sin𝜃𝐹𝑧1𝜃,𝜃2;𝑓𝑥,𝑓𝑦𝑑𝜃=0,(3.41)𝜓𝜃2𝜃1(cos𝜃1)𝐹𝑧1𝜃,𝜃2;𝑓𝑥,𝑓𝑦𝑑𝜃=̂𝑢2,(3.42) where a function 𝐹𝑧1 is given by 𝐹𝑧1𝜃,𝜃2;𝑓𝑥,𝑓𝑦=12𝑓𝑥cos𝜃2𝑓cos𝜃2𝑦sin𝜃2sin𝜃.(3.43)

By introducing the same type of variable transformations as employed in the previous case, (3.40)–(3.42) becomeΓ𝑜𝜃1,𝜃2,𝑓𝑥,𝑓𝑦=𝜙𝜋/21𝑓𝑜𝜙;𝑝𝑓𝑑𝜙𝑠=0,(3.44)Γ𝑣𝜃1,𝜃2,𝑓𝑥,𝑓𝑦=𝜙𝜋/21𝑓𝑣𝜙,𝜃𝑜;𝑝𝑑𝜙=0,(3.45)𝜃1,̂𝑢2,𝑓𝑥=𝜃1,𝜃2,̂𝑢2,𝑓𝑥,𝑓𝑦=̂𝑢21𝑓𝑠𝜙𝜋/21𝑓𝑢𝜙,𝜃𝑜;𝑝𝑑𝜙.(3.46) Since the right-end moment is prescribed (i.e., 𝑚2=0), the rotation at the right end 𝜃2 is no longer an independent quantity but can be obtained in terms of {𝜃1,𝑓𝑥}.

Let us redefine 𝐟=[𝐟𝑝𝐟𝑟]𝑇 where 𝐟𝑝𝑓={𝑥2,𝑚1,} and 𝐟𝑟𝑓={𝑥1,𝑓𝑦1,𝑓𝑦2} and redefine 𝐮={̂𝑢2,𝜃1,𝑓𝑥}𝑇. Consistent with these new definitions, the reduced gradient matrix takes the form 𝐠=[𝐠𝑇𝑝𝐠𝑇𝑟]𝑇 where the submatrices 𝐠𝑝 and 𝐠𝑟 are given by𝐠𝑝=𝜕𝑓𝑥2𝜕̂𝑢2𝜕𝑓𝑥2𝜕𝜃1𝜕𝑓𝑥2𝜕𝑓𝑥𝜕𝑚1𝜕̂𝑢2𝜕𝑚1𝜕𝜃1𝜕𝑚1𝜕𝑓𝑥𝜕𝜕̂𝑢2𝜕𝜕𝜃1𝜕𝜕𝑓𝑥,𝐠𝑟=1𝑑𝑔11𝑑𝑔12𝑑𝑔13𝑑̂𝑠𝑔22𝑔23̂𝑠𝑔22𝑔23,(3.47) where 𝑔𝑖𝑗 denotes any entry of the submatrix 𝐠𝑝 and ̂𝑠=𝑚1/𝑑. As clearly indicated in the previous section, the entries 𝑔11,𝑔12, and 𝑔21 vanish and 𝑔13=𝑔31=1. The remaining entries are contained in a submatrix 𝓰𝑝 given by𝓰𝑝=𝜕𝑚1𝜕𝜃1𝜕𝑚1𝜕𝑓𝑥𝜕𝓡𝜕𝜃1𝜕𝓡𝜕𝑓𝑥.(3.48) By imposing the remaining moment boundary conditions at the left end of the member, we obtain𝑚21𝑓+2𝑥cos𝜃1cos𝜃2𝑓2𝑦sin𝜃1sin𝜃2=0.(3.49)

Taking derivative of (3.49) with respect to 𝜃1 and 𝑓𝑥 leads to𝜕𝑚1𝜕𝜃1𝜕𝑚1𝜕𝑓𝑥=𝐌1+𝐌2𝐀,(3.50) where matrices 𝐌1, 𝐌2, and 𝐀 are given by𝐌1=𝜆1𝑚1𝑐2𝑐1𝑚1,𝐌2=𝜆2𝑚1𝑠1𝑠2𝑚1,𝐀=𝜕𝜃2𝜕𝜃1𝜕𝜃2𝜕𝑓𝑥𝜕𝑓𝑦𝜕𝜃1𝜕𝑓𝑦𝜕𝑓𝑥.(3.51) By differentiating (3.46) with respect to 𝜃1 and 𝑓𝑥, it results in𝜕𝓡𝜕𝜃1𝜕𝓡𝜕𝑓𝑥=𝐁+𝐂𝐀,(3.52) where matrices 𝐁 and 𝐂 are𝜕𝐁=𝓡𝜕𝜃1𝜕𝓡𝜕𝑓𝑥,𝜕𝐂=𝓡𝜕𝜃2𝜕𝓡𝜕𝑓𝑦.(3.53) Note that all entries of the matrices 𝐁 and 𝐂 can be obtained directly from (3.46) along with the transformations 𝜃=𝜋+(𝜃𝜃𝑜) and 𝑝sin𝜙=sin(𝜃/2) (see explicit results in Appendix C). To compute all entries of the matrix 𝐀, (3.44) and (3.45) are differentiated with respect to 𝜃1 and 𝑓𝑥 and this yields𝐃𝐀=𝐅,(3.54) where matrices 𝐃 and 𝐅 are given by𝜕𝐃=Γ𝑜𝜕𝜃2𝜕Γ𝑜𝜕𝑓𝑦𝜕Γ𝑣𝜕𝜃2𝜕Γ𝑣𝜕𝑓𝑦,𝜕𝐅=Γ𝑜𝜕𝜃1𝜕Γ𝑜𝜕𝑓𝑥𝜕Γ𝑣𝜕𝜃1𝜕Γ𝑣𝜕𝑓𝑥.(3.55) Explicit expressions for all entries of the matrices 𝐃 and 𝐅 are reported in Appendix C. Once the matrix 𝐀 is solved from (3.54), it is substituted into (3.50) and (3.52) to obtain all entries of the matrix 𝓰𝑝. Again, due to the complexity of all functions resulting from variable transformations, the matrices 𝐁,𝐂,𝐃, and 𝐅 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 Figure 3. Let {𝑥,𝑦} be a local coordinate system of the undeformed member and {𝑥,𝑦} be the coordinate system of the deformed member defined such that a chord connecting its end points always lies on the 𝑥 axis. With this specific choice of {𝑥,𝑦}, behavior of the member observed from this coordinate system is identical to that of a simply supported member.

The normalized end loads and normalized end displacements and rotations are denoted by {𝑓𝑥1,𝑓𝑦1,𝑚1,𝑓𝑥2,𝑓𝑦2,𝑚2} and {̂𝑢1,̂𝑣1,𝜃1,̂𝑢2,̂𝑣2,𝜃2} in the {𝑥,𝑦} coordinate system and by {𝑓𝑥1,𝑓𝑦1,𝑚1,𝑓𝑥2,𝑓𝑦2,𝑚2} and {̂𝑢2,𝜃1,𝜃2} in the {𝑥,𝑦} coordinate system. From geometric consideration of the deformed configuration, {̂𝑢2,𝜃1,𝜃2} can readily be expressed in terms of {̂𝑢1,̂𝑣1,𝜃1,̂𝑢2,̂𝑣2,𝜃2} by𝜃1=𝜃1𝜃+𝜙,2=𝜃2+𝜙,̂𝑢2=1+̂𝑢2̂𝑢1̂𝑣cos𝜙+2̂𝑣1sin𝜙1,(3.56) where 𝜙 is a chord rotation governed by1+̂𝑢2̂𝑢1̂𝑣sin𝜙2̂𝑣1cos𝜙=0.(3.57)

Let 𝑓𝑥 be the internal force in 𝑥 direction and 𝓡 be the residual defined in the {𝑥,𝑦} system in the same way as given by (3.10). In the {𝑥,𝑦} coordinate system, we choose {𝑓𝑥,} such that𝑓𝑥=𝑓𝑥,=.(3.58)

By applying the coordinate transformation, a vector 𝑓𝐟={𝑥1,𝑓𝑦1,𝑚1,𝑓𝑥2,𝑓𝑦2,𝑚2,} is related to a vector 𝐟𝑓={𝑥2,𝑚1,𝑚2,,𝑓𝑥1,𝑓𝑦1,𝑓𝑦2} by𝐟=𝐑(𝜙)𝐟,(3.59) where 𝐑(𝜙) is a transformation matrix of dimensions 7×7 given by𝐑(𝜙)=0000𝑐𝜙𝑠𝜙00000𝑠𝜙𝑐𝜙0𝑐0100000𝜙00000𝑠𝜙𝑠𝜙00000𝑐𝜙00100000001000,(3.60) in which 𝑠𝜙=sin𝜙 and 𝑐𝜙=cos𝜙. By defining 𝐮={̂𝑢1,̂𝑣1,𝜃1,̂𝑢2,̂𝑣2,𝜃2,𝑓𝑥} and 𝐮={̂𝑢2,𝜃1,𝜃2,𝑓𝑥} and then recalling (3.56)–(3.58), we obtain the relation 𝐮=𝐮(𝐮). From the fact that behavior of the member in the {𝑥,𝑦} system is identical to that for a simply supported member, 𝐟 and 𝐮 are therefore related by 𝐟=𝐟(𝐮). Combining (3.59), 𝐮=𝐮(𝐮) and 𝐟=𝐟(𝐮) leads to the relation 𝐟=𝐟(𝐮)=𝐑(𝜙)𝐟(𝐮(𝐮)), and, from Taylor series expansion, this nonlinear function possesses the best linear approximation in the neighborhood of a vector 𝐮0𝐟𝐮(𝐮)=𝐟0+𝐤𝑙𝐮0𝐮𝐮0,(3.61) where 𝐤𝑙 is a local element tangent stiffness matrix given by𝐤𝑙=𝜕𝐑𝐟𝜕𝜙𝜕𝜙𝜕𝐮+𝐑𝐠𝜕𝐮𝜕𝐮,(3.62) in which the relation 𝐠=𝜕𝐟/𝜕𝐮 is utilized.

For a special case in which the member contains an inflection point at its right end, the end moment 𝑚2 vanishes and the corresponding end rotation 𝜃2 is eliminated from the set of unknowns. Let us first define following reduced vectors 𝑓𝐟={𝑥1,𝑓𝑦1,𝑚1,𝑓𝑥2,𝑓𝑦2,}, 𝐟𝑓={𝑥2,𝑚1,,𝑓𝑥1,𝑓𝑦1,𝑓𝑦2}, 𝐮={̂𝑢1,̂𝑣1,𝜃1,̂𝑢2,̂𝑣2,𝑓𝑥}, and 𝐮={̂𝑢2,𝜃1,𝑓𝑥}. By applying the coordinate transformation, the relationship between the reduced vectors 𝐟 and 𝐟 is given by𝐟=𝐑(𝜙)𝐟,(3.63) where 𝐑(𝜙) is a reduced transformation matrix of dimensions 6×6 given by𝐑(𝜙)=000𝑐𝜙𝑠𝜙0000𝑠𝜙𝑐𝜙0𝑐010000𝜙0000𝑠𝜙𝑠𝜙0000𝑐𝜙001000.(3.64) The relation 𝐮=𝐮(𝐮) is then obtained by recalling (3.56)–(3.58) and 𝐟=𝐟(𝐮) results from the fact that the behavior of the member in the {𝑥,𝑦} system is identical to that for a simply supported member. Combining (3.63), 𝐮=𝐮(𝐮) and 𝐟=𝐟(𝐮) yields 𝐟=𝐑(𝜙)𝐟(𝐮(𝐮)), and, from Taylor series expansion, the best linear approximation of this nonlinear function in the neighborhood of a vector 𝐮0 takes the form𝐟𝐮=𝐟𝐮0+𝐤𝑙𝐮0𝐮𝐮0,(3.65) where 𝐤𝑙 is a reduced local element tangent stiffness matrix given by𝐤𝑙=𝜕𝐑𝜕𝜙𝐟𝜕𝜙𝜕𝐮+𝐑𝐠𝜕𝐮𝜕𝐮,(3.66) in which the relation 𝐠=𝜕𝐟/𝜕𝐮 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 𝑋-axis to the local 𝑥-axis (defined in the Section 3.3). The element tangent stiffness matrix 𝐤𝑔 referring to the global coordinate system {𝑋,𝑌} can be obtained using the following standard coordinate transformation𝐤𝑔=𝐐𝑇𝐤𝑙𝐐,(3.67) where 𝐐 is a transformation matrix of dimension 7×7 given by𝑐𝐐=𝛽𝑠𝛽00000𝑠𝛽𝑐𝛽000000010000000𝑐𝛽𝑠𝛽00000𝑠𝛽𝑐𝛽0000000100000001,(3.68) in which 𝑠𝛽=sin𝛽 and 𝑐𝛽=cos𝛽. Similarly, for a member containing an inflection point at the right end, we obtain𝐤𝑔=𝐐𝑇𝐤𝑙𝐐,(3.69) where 𝐤𝑔 is the reduced, global element tangent stiffness matrix and 𝐐 is a reduced transformation matrix given by𝑐𝐐=𝛽𝑠𝛽0000𝑠𝛽𝑐𝛽0000001000000𝑐𝛽𝑠𝛽0000𝑠𝛽𝑐𝛽0000001.(3.70)

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) where 𝐏 is a vector of nodal loads and zero residuals of all members, 𝐔 is a vector of nodal displacements, nodal rotations and the internal axial force of all members, {𝐏𝑜,𝐔𝑜} are vectors {𝐏,𝐔} at the reference state, and 𝐊𝑡 is the structure tangent stiffness matrix. Note that the matrix 𝐊𝑡 can readily be obtained from a direct assembly of the global element tangent stiffness matrices 𝐤𝑔 or 𝐤𝑔 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 length 𝐿 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𝐿/EI{1,2,3,4,5}.

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𝜃=0.5𝑚𝜉,0𝜉0.5,𝑚(0.75𝜉),0.5𝜉1,̂𝑢=2sin0.5𝑚𝜉𝑚𝜉,0𝜉0.5,sin0.25𝑚+sin(0.75𝜉)𝑚̂𝑚+𝜉1,0.5𝜉1,𝑣=22cos0.5𝑚𝜉(𝑚,0𝜉0.5,23cos0.25𝑚+cos0.75𝜉)𝑚𝑚,0.5𝜉1.(5.1) The deflected shapes of the beam for different values of 𝑚 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.

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 length 𝐿 and flexural rigidity EI, and the frame is subjected to three concentrated moments {𝑀1,𝑀2,𝑀3} as shown schematically in Figure 6(a). In the analysis, we choose {𝑀1,𝑀2,𝑀3} such that 𝑀1=𝑀2=𝛼EI/𝐿 and 𝑀3=2.5𝛼EI/𝐿, 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.

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.

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 Figure 8(a). The frame consists of a horizontal member, two vertical members of the same length 𝐿 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).

The deflected shapes of the structure are reported in Figure 9 for many values of normalized vertical loads 𝑓=𝐹𝐿2/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.

5.4. Diamond Box Frame Subjected to Vertical Force

As a final example, consider a diamond box frame subjected to a vertical downward force 𝐹 as shown schematically in Figure 10(a). The frame consists of four 45-degree bevel elements of the same length 2𝐿, 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., 𝑓𝑦=𝐹𝐿2/EI{5,10,15,20}) are treated.

Deflected shapes of a structure for different 𝑓𝑦 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 𝑓𝑦 and then starts to increase rapidly for large 𝑓𝑦. This is due to that for small 𝑓𝑦, 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 𝑓𝑦 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.

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 dimension 7×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.

Appendices

In this section, we demonstrate the derivation of gradient matrices in Section 3.2. The gradient matrix 𝓰𝑝 of a member containing no inflection point (Section 3.2.1) is presented in an explicit form, while, in the case of a member containing an inflection point (Section 3.2.2 and Section 3.2.3), we expand all entries of matrices 𝐁,𝐂,𝐃,𝐅,𝐁,𝐂,𝐃, and 𝐅.

A. Member Containing No Inflection Point

Let us refer to (3.20) and (3.21); the explicit form of matrices 𝐒 and 𝐃 are given by𝜗𝑖𝐒=2𝑖1𝑠2𝑚1𝑠2𝑠1𝜗𝑖2𝑖1𝑠1𝑚2𝑠2𝑠10𝜗𝑖4𝑖2𝑠2𝑚1𝑠2𝑠1𝜗𝑖4𝑖2𝑠1𝑚2𝑠2𝑠10𝜗𝑖5𝑖3𝑠2𝑚1𝑠2𝑠1𝜗𝑖5𝑖3𝑠1𝑚2𝑠2𝑠11,𝜗𝑖𝐃=2𝑖1𝑠2𝜆1𝑠2𝑠1𝜗𝑖2𝑖1𝑠2𝜆2𝑠2𝑠1𝜗𝑖2𝑐1𝑐2+𝑖1𝑠12𝑠2𝑠1𝜗𝑖4𝑖2𝑠2𝜆1𝑠2𝑠1𝜗𝑖4𝑖2𝑠2𝜆2𝑠2𝑠1𝜗𝑖4𝑐1𝑐2+𝑖2𝑠12𝑠2𝑠1𝜗𝑖5𝑖3𝑠2𝜆1𝑠2𝑠1𝜗𝑖5𝑖2𝑠3𝜆2𝑠2𝑠1𝜗𝑖5𝑐1𝑐2+𝑖3𝑠12𝑠2𝑠1+1𝑚11𝑚2𝜗𝑖3𝑠1𝑚1𝑠2𝑚2𝜗𝑖5𝑐1𝑚1𝑐2𝑚2𝜗𝑖6,(A.1) where 𝑠1=sin𝜃1,𝑠2=sin𝜃2, 𝑠12=sin(𝜃1𝜃2), 𝑐1=cos𝜃1, 𝑐2=cos𝜃2, 𝜆1=𝑓𝑥𝑠1+𝑓𝑦𝑐1, 𝜆2=𝑓𝑥𝑠2+𝑓𝑦𝑐2, and integrals {𝑖1,𝑖2,𝑖3,𝑖4,𝑖5,𝑖6} are defined by𝑖1=𝜃2𝜃1𝐹3𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑖𝑑𝜃,2=𝜃2𝜃1𝐹3𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑖sin𝜃𝑑𝜃,3=𝜃2𝜃1𝐹3𝜃,𝜃2;𝑓𝑥,𝑓𝑦,𝑚2𝑖cos𝜃𝑑𝜃,4=𝜃2𝜃1