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., [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 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., [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 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 , 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 ” and “point ” are often utilized for convenience in this paper to refer to a cross-section at any state with its centroid located at a point in the undeformed state. During a particular motion, the cross-section at 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 , 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 , respectively.
(a)
(b)
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 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 by 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: . By using these two relations, (2.3) becomes where nondimensional parameters appearing in above equation are defined by , , and . By performing a direct integration of (2.5), it leads to 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 where is a moment-dependent function defined such that if and if . By combining (2.4) and (2.7), it leads to following two relations: 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 and a roller support at its right end ) subjected to end loads as shown in Figure 2. By imposing the moment boundary condition at the right end, that is, , we obtain the constant as By inserting the constant into the relations (2.7)-(2.8), it yields where the function is defined by By imposing the remaining two natural boundary conditions at the left and right ends of a member, we obtain where and . Normalized support reactions can readily be computed from equilibrium of the entire member, and final results are given by where and . Next, by integrating (3.2) from to , it leads to a system of three nonlinear algebraic equations For a given set of end loads , the end displacement and rotations 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 for a particular set of prescribed displacement and rotations 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, cannot be specified arbitrarily. However, if are prescribed instead, the end loads can be obtained by solving (3.4), (3.5), (3.7), and (3.8) simultaneously, and the end displacement can subsequently be computed from (3.9). However, lack of the displacement component renders a set not well suited for the development of a solution procedure by a direct stiffness method.
In the present study, is chosen as a set of primary unknowns. It is worth noting that to allow to be one of independent variables, the strong requirement posed by (3.9) must be relaxed via the introduction of the residual such that Supplemented by (3.10), for any given set , the quantities 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 and the corresponding are solutions of the boundary value problem if and only if the residual vanishes.
Let be a vector defined by where and and let be a vector defined by . 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 takes a form where a gradient matrix is given by The submatrix is expressed explicitly in terms of differentiations as 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 where . As clearly indicated by (3.5), (3.7), (3.8), and (3.10), entries , and vanish whereas . The remaining entries of are contained in a submatrix given by 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., ) or strictly negative (i.e., ) for the entire beam (i.e., ) when the applied end moments are nonzero and possess different sign (i.e., ). 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, for and for . 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 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 It is noted that the relation (3.4) implicitly defines in terms of and . As a result, by taking derivative of (3.16)–(3.18) with respect to along with employing the chain rule of differentiation, it leads to where the matrices and are given by 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 or, equivalently, the bending moment vanishing at and for and . This particular case arises when the applied end moments are nonzero and possess the same sign (i.e., ). 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 , it results in for and for , and for , it results in for and for .
For this special case, the constant appearing in (2.6) is alternatively obtained by imposing a condition at the inflection point, that is, , and this leads to where is the rotation at the inflection point. Upon using (3.22), the relations (2.7)-(2.8) become where the function is defined by By integrating equations (3.23) over the entire member, we obtain where a constant is defined such that if and if .
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 , and , (ii) applying change of variable and identity , and (iii) introducing another variable transformation with . The function simply reduces to where . Finally, the nonlinear algebraic equations (3.25) can be expressed in terms of integrals over a dummy variable as where , , , , and
By enforcing the remaining two moment boundary conditions at both ends of the member, it leads to By taking derivative of (3.31) with respect to , , and , it results in where matrices , , and are given by where , , , , , , , and . Differentiating (3.29) with respect to , and yields where matrices and are given by Note that all entries of the matrices B and C can be obtained directly from (3.29) along with the change of variables and (see explicit results in Appendix B).
To compute all entries of the matrix , we differentiate (3.27) and (3.28) with respect to , and , and this results in where matrices and are given by Similarly, all entries of the matrices and can be obtained directly from (3.27) and (3.28) along with the change of variables and (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 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 . 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 and . 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 into (3.25), it leads to where a function is given by
By introducing the same type of variable transformations as employed in the previous case, (3.40)–(3.42) become Since the right-end moment is prescribed (i.e., ), the rotation at the right end is no longer an independent quantity but can be obtained in terms of .
Let us redefine where and and redefine . Consistent with these new definitions, the reduced gradient matrix takes the form where the submatrices and are given by where denotes any entry of the submatrix and . As clearly indicated in the previous section, the entries , and vanish and . The remaining entries are contained in a submatrix given by By imposing the remaining moment boundary conditions at the left end of the member, we obtain
Taking derivative of (3.49) with respect to and leads to where matrices , , and are given by By differentiating (3.46) with respect to and , it results in where matrices and are Note that all entries of the matrices and can be obtained directly from (3.46) along with the transformations and (see explicit results in Appendix C). To compute all entries of the matrix , (3.44) and (3.45) are differentiated with respect to and and this yields where matrices and are given by 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 and in the coordinate system and by and in the coordinate system. From geometric consideration of the deformed configuration, can readily be expressed in terms of by where is a chord rotation governed by
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
By applying the coordinate transformation, a vector is related to a vector by where is a transformation matrix of dimensions given by in which and . By defining and 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 where is a local element tangent stiffness matrix given by 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 vanishes and the corresponding end rotation is eliminated from the set of unknowns. Let us first define following reduced vectors , , , and . By applying the coordinate transformation, the relationship between the reduced vectors and is given by where is a reduced transformation matrix of dimensions given by 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 takes the form where is a reduced local element tangent stiffness matrix given by in which the relation is utilized. Note that the local element tangent stiffness matrix for this particular case is of dimensions .
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 where is a transformation matrix of dimension given by in which and . Similarly, for a member containing an inflection point at the right end, we obtain where is the reduced, global element tangent stiffness matrix and is a reduced transformation matrix given by
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 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, .
(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 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 as shown schematically in Figure 6(a). In the analysis, we choose such that and , 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)
(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.
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).
(a)
(b)
The deflected shapes of the structure are reported in Figure 9 for many values of normalized vertical loads . 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 , 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., ) are treated.
(a)
(b)
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 . 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 where , , , , , , and integrals are defined by By substituting (A.1) into (3.19), and then solving such system of linear equations, it leads to the explicit forms of the matrix where , and .
B. Member Containing Interior Inflection Point
This section presents explicit results for all entries of matrices , and . By employing chain rule of differentiation and recalling all involved identities and changes of variables, entries of the matrices , and are given by By taking derivatives of functions with respect to , we obtain By taking derivatives of functions with respect to , we obtain All entries of matrices , and can now be obtained via the relation (B.1) along with the results (B.2) and (B.3).
C. Member Containing Inflection Point at Its Right End
This section presents explicit results for all entries of matrices , , , and . By employing chain rule of differentiation and recalling all involved identities and changes of variables, entries of the matrices , , , and are given by By taking derivatives of functions with respect to , we obtain By taking derivatives of functions with respect to , we obtain All entries of matrices , , , and can now be obtained via the relation (C.1) along with the results (C.2) and (C.3).
Acknowledgments
The first author gratefully acknowledges financial supports provided by “Chulalongkorn University for development of new faculty staff” and “Stimulus Package 2 (SP2) of Ministry of Education under the theme of Green Engineering for Green Society”.