Abstract
This paper presents an efficient and accurate numerical technique for analysis of two-dimensional frames accounted for both geometric nonlinearity and nonlinear elastic material behavior. An adaptive remeshing scheme is utilized to optimally discretize a structure into a set of elements where the total displacement can be decomposed into the rigid body movement and one possessing small rotations. This, therefore, allows the force-deformation relationship for the latter part to be established based on small-rotation-based kinematics. Nonlinear elastic material model is integrated into such relation via the prescribed nonlinear moment-curvature relationship. The global force-displacement relation for each element can be derived subsequently using corotational formulations. A final system of nonlinear algebraic equations along with its associated gradient matrix for the whole structure is obtained by a standard assembly procedure and then solved numerically by Newton-Raphson algorithm. A selected set of results is then reported to demonstrate and discuss the computational performance including the accuracy and convergence of the proposed technique.
1. Introduction
It has been well recognized that various existing conventional techniques based either on linearized or simplified kinematics are generally incapable of modeling several crucial structural behaviors (e.g., buckling and postbuckling responses, -delta effects, stability characteristics and bifurcations, etc.) or yield inaccurate results when they are applied to the modeling of structures accounted for highly geometric nonlinearities. For this reason, integration of a proper nonlinear kinematics to gain a more physically suitable mathematical model that can simulate practical situations to the level of complexity involved becomes essential. Analysis software packages fully equipped with sophisticated mathematical models have been emerged continuously with the primary objectives to enhance the capability of the prediction of structural responses under various circumstances, for example, presence of strong nonlinearities, various types of practical excitations, and intricate interactions with surrounding environments.
Treatment of geometric nonlinearities in structural modeling has a very long history since the age of well-known mathematicians such as Euler, Lagrange, and Kirchoff (see extensive historical review in [1]). In that era, a problem of finding exact deformed shape of beams and columns, commonly known as an elastica problem, became very popular and mathematically challenging. However, due to lack of powerful, computationally aided tools and complexity of governing differential equations posed by highly nonlinear kinematics, solutions of most problems were derived based on analytical techniques and, therefore, limited only to simple structural configurations and loading conditions.
As a result of the invention of computational devices and breakthrough in the field of numerical analysis, various powerful and robust numerical techniques have been continuously proposed, for past several decades, and employed to solve more practical and large-scale elastica problems. Some of those investigations relevant to the present study are briefly summarized below not only to indicate the historical background and the recent advances of computational techniques in this specific area but also to identify the gap of knowledge and the novel aspect of the current investigation. Lewis and Monasa [2] investigated the deformed shape of a cantilever beam which is made of nonlinear elastic materials of the Ludwik type (i.e., the uniaxial stress-strain relationship follows a pure power law) and subjected to a concentrated force by using the fourth-order Runge-Kutta method. Later, Monasa and Lewis [3] generalized their previous work to study a Ludwik-type cantilever beam under multiple concentrated forces. Prathep and Varadan [4] developed an iterative numerical scheme to investigate the inelastic postbuckling behavior of a cantilever column. In their work, the material nonlinearity was modeled by using Ramberg-Osgood stress-strain relation. Kounadis and Mallis [5] studied the influence of the slenderness ratio and a material parameter contained in the nonlinear elastic constitutive model on the buckling behavior of a bar. In 1989, a large deflection of a linearly elastic cantilever beam under a rotational, arbitrarily distributed load was studied by B. N. Rao and G. V. Rao [6] using the fourth-order Runge-Kutta technique. Wang [7] proposed a numerical procedure based on a perturbation technique to investigate the postbuckling behavior of a prismatic, roller-clamped, linearly elastic column under the applied end load. Later, Lee and Oh [8] used Runge-Kutta and Regula-Falsi methods to solve the elastica and buckling of a simply supported tapered column subjected to the end load. The postbuckling behavior of a linearly elastic cantilever column under the combined action of a uniformly distributed load and a concentrated force at the tip was also studied by Lee [9]. In the analysis, a technique based on Butcher’s fifth-order Runge-Kutta method was employed to obtain numerical results. Later, Lee [10] employed the shooting method and bisection algorithm to determine the large deflection of a nonlinear elastic cantilever beam under a uniformly distributed load and a concentrated force at its free end. Similar to the work of Lewis and Monasa [2] and Monasa and Lewis [3], the material model of Ludwik-type was assumed in the analysis. Dado et al. [11] investigated the computational performance of various techniques such as a semianalytical method, a numerical integration scheme, and a finite element method in the postbuckling analysis of a column consisting of two segments linked by a rotational spring. In 2005, Dado and Al-sadder [12] studied the large deflection of prismatic and nonprismatic cantilever beams under various types of applied loads. In their work, a technique based on the minimization of the residual error of the governing differential equation along with the representation of the beam rotation in terms of polynomial functions was employed. Vaz and Patel [13] applied both the shooting and Runge-Kutta methods to study the postbuckling behavior of a simply supported column. The material nonlinearity was handled by assuming a bilinear moment-curvature relationship. The influence of the rotational spring at the base of a cantilever column and the tip follower force on its deformed shape was also investigated by Shvartsman [14] by using a technique requiring no iterative procedure. Lacarbonara [15] investigated the postbuckling behavior of a nonprismatic nonlinearly elastic rod by using the higher-order perturbation technique. Later, Wang et al. [16] proposed an efficient and accurate technique, known as Homotopy analysis method (HAM), to solve a large deflection of a cantilever beam subjected to a concentrated force at the free end. Banerjee et al. [17] employed the shooting method along with the adomian decomposition to study a linear elastic cantilever beam under arbitrary loading conditions and containing an interior inflection point. The influence of discontinuity conditions and initial displacements on responses of frames undergoing large displacement and rotation was studied by Hu et al. [18]. In their analysis, a technique called the differential quadrature element method (DQEM) was utilized to construct numerical solutions. Lanc et al. [19] performed a large-displacement analysis of beam-type structures using the one-dimensional nonlinear finite element technique. In their work, the material nonlinearity was treated by using an elastic-plastic constitutive model. Later, Shatarat et al. [20] investigated a nonlinear elastic, rhombus frame subjected to a pair of opposite forces acting along its diagonal. In their study, the relationship between the displacement and applied forces at the corner was obtained by using a semianalytical technique. Chen [21] proposed a technique based on the moment integral treatment to determine the large deflection of a cantilever beam and results were found to be in good agreement with those obtained by the elliptic integral method and the commercial finite element package ANSYS©.
Recently, Rungamornrat and Tangnovarad [22] proposed a semianalytical technique based on a direct stiffness procedure to analyze two-dimensional, linearly elastic, inextensible frames of arbitrary geometry and undergoing large displacement and rotation. In their formulation, the element tangent stiffness matrix and the residual load vector were derived in an exact form using the conventional elliptic integral technique. Later, Douanevanh et al. [23] generalized the work of Rungamornrat and Tangnovarad [22] to remove the inextensibility assumption. It is important to remark that while their technique yielded highly accurate results comparable to analytical solutions without the requirement of mesh refinement, it still requires solving a system of exact nonlinear algebraic equations in the calculation of the element tangent stiffness matrix and evaluation of the load residuals. In addition, that system of equations is strongly dependent on the presence of inflection points within elements. More recently, Sinsamutpadung et al. [24] developed a simplified numerical technique to solve the large displacement and rotation of two-dimensional frames. Their proposed strategy was based primarily on the displacement decomposition and a linearized kinematics for the displacement measured in the corotational coordinate system of an element. While their technique was found promising and robust, the formulation is still limited to frames made of linearly elastic materials.
For various structures found in practices, as their displacement induced by excitations becomes sufficiently large (e.g., structures in the vicinity of collapses), it generally introduces large deformations or strains in certain regions and, as a result, linear constitutive models cannot be utilized to accurately represent the material behavior in this regime. In this paper, the technique proposed by Sinsamutpadung et al. [24] is enhanced by incorporating the material nonlinearity in the analysis of frames undergoing large displacement and rotation. The following sections present basic equations and assumptions, element force-displacement relations in different reference coordinate systems, load-displacement relationship for a structure, solution procedures, a selected set of numerical results, and conclusions and remarks.
2. Basic Equations for Small-Rotation-Based Kinematics
Consider a straight and prismatic member of length , subjected to applied loads at both ends as shown schematically in Figure 1(a). The member is made of a homogeneous, hypothetical material with the moment-curvature and axial-force-strain relationships at any cross section prescribed a priori. An undeformed state of this member, represented by a line connecting the centroid of all cross sections and its subsequent deformed state resulting from the action of applied loads, is also indicated in Figure 1(a).
(a)
(b)
By assuming that the member undergoes small displacement and rotation relative to its undeformed state, the static equilibrium equations can be formulated based on the known undeformed geometry (see Figure 1(b) for a free-body diagram of the undeformed element ) and it finally leads to where , , and are the axial force, shear force, and bending moment at any cross section, respectively. Components of the displacement in the - and -directions, denoted by and , are related to the axial strain and curvature via the small-rotation-based kinematics: It should be remarked that the relation (2) results directly from the approximations and when the rotation is sufficiently small. The axial strain and curvature of any cross section are assumed to be related to the axial force and the bending moment via the fully uncoupled relations: where denotes the cross sectional area of a member, is Young’s modulus, and is any prescribed, sufficiently smooth function. It is worth pointing out that a form of the function generally depends on the stress-strain relationship and shape of the cross section. In the present study, the axial force is assumed to be related linearly to the axial strain and its influence on the moment-curvature relationship is neglected. It should be noted that the use of this simplified constitutive relation for the cross section does not aim to represent the real material behavior but is only for convenience purpose to demonstrate the computational performance of the proposed technique when dealing with nonlinear elasticity. The technique can be further extended to incorporate more complex material models by only changing the constitutive relation (3).
3. Element Force-Displacement Relationship
In this section, the force-displacement relationship and its best linear approximation for a generic member undergoing large displacement and rotation are developed. By following a well-known corotational technique, the total displacement of a member (measured from the undeformed state to its deformed state ) can be decomposed into a deformation-free displacement (measured from to a state which results directly from the rigid translation and rotation) and a displacement measured from to as illustrated in Figure 2. Three different reference coordinate systems, that is, a global system , a local system , and a corotational system are introduced for further reference. It should be remarked that is a moving coordinate system with its origin always located at the end of a deformed member and the -axis coincident with a line connecting the two end points (see Figure 2).
A vector containing the total displacements and rotations at both ends of the member in the global coordinate system is defined by , where the subscript indicates the member end, the symbols and denote the displacement along the -and -axes, and denotes the end rotation. Similarly, a vector of the total end displacements and rotations of the member in the local coordinate system can be defined by . From a standard law of coordinate transformations, can readily be related to by where denotes a conventional transformation matrix which is a function of a member orientation angle . From the geometric consideration, a vector containing nonzero components of the end displacement and rotations of the member measured relative to the moving coordinate system , denoted by where , , and are clearly defined in Figure 2, can be expressed, in terms of and the angle (i.e., an angle between the -axis and the -axis), as where is the initial length of the member and the angle can be obtained in terms of (i.e., ) from the following relation: Vectors containing the end forces and end moments of each member in the global and local coordinate systems, denoted, respectively, by and , where all entries are ordered in a fashion consistent with and , can also be related, again, via the standard law of coordinate transformations: Similarly, a member force vector measured in the corotational coordinate system, denoted by , can be related to by where is a standard transformation matrix which is clearly dependent on the angle .
3.1. Force-Displacement Relation in Corotational Coordinate System
To establish the relationship between the force vector and the vector , it is assumed that the member length resulting from the discretization is sufficiently small such that the rotation at any point in the current state relative to the state is infinitesimal. As a result, the governing equations (1) and (2) apply. It is worth noting that the movement of a member perceived by an observer in the corotational system is identical to that of a member restrained by a pinned support at the origin and the roller support at the other end and subjected to applied end loads as shown in Figure 3.
By forming equilibrium equations for the entire member in its deformed state, all support reactive forces can readily be obtained in terms of the applied loads and the elongation by By employing the first equation of (1), (2), and (3) along with the end boundary conditions, the applied force can be related linearly to by To establish the relation between the applied end moments and the end rotations , a technique based on the generalized moment area or curvature area equations is employed as briefly summarized below (also see details in the work of Pinyochotiwong et al. [25]). First, by using the equilibrium equations (1), the bending moment at any cross section can be expressed in terms of the applied moments as where and are linear functions of defined by The curvature of any cross section can then be obtained by substituting (11) into the second equation of (3). Finally, the end rotations are obtained from the following two curvature area equations: Equations (13) implicitly define the end moments in terms of the end rotations ; that is, and . From Taylor series expansion of (13), the best linearized relationship between and is given by where and are end moments and end rotations at a reference equilibrium state and is a flexural tangent flexibility matrix of a member with its entries explicitly defined by Relations (9), (10), and (13) implicitly define the force vector as a nonlinear function of the displacement vector ; that is, . The best linear approximation of this nonlinear function is given by where and are displacement and force vectors at a reference equilibrium state and is the gradient matrix evaluated at . Entries of can be computed as follows: is obtained directly from (10), is clearly the inverse of the matrix , , and can readily be obtained in terms of other computed entries via the relation (9).
3.2. Force-Displacement Relation in Local Coordinate System
By combining (5), (6), (8), and the implicit function constructed in Section 3.1, it leads to the member force-displacement relationship in the local coordinate system The best linear approximation of (17) can be obtained, again, by Taylor series expansion with the final result: where and are local displacement and force vectors at a reference equilibrium state and is a local element tangent stiffness matrix evaluated at . The form of can be obtained from the chain rule for differentiations and the final expression can be written in a concise form as
where the involved submatrices , , and can readily be computed using the relations (4)-(5) and the fact that the transformation matrix is given explicitly in terms of elementary functions of (see explicit expression of these submatrices in the work of Sinsamutpadung et al. [24]). It can also be verified that is essentially symmetric.
3.3. Force-Displacement Relation in Global Coordinate System
By combining (4), (7), and (17), it yields the member force-displacement relationship in the global coordinate system: Clearly, this relationship allows the end forces and end moments to be determined for any prescribed end displacements and rotations . Again, by performing Taylor series expansion about any admissible vector , a linearized relationship between and in the neighborhood of can be obtained as where and is a global element tangent stiffness matrix evaluated at . From the chain rule for differentiations and the fact that is a constant, it can be verified that
4. Structure Load-Displacement Relation
To establish the structure load-displacement relationship, the given structure is first discretized into straight and prismatic members. By defining and as a vector containing all nodal degrees of freedom (i.e., nodal displacements and rotations) and a vector of applied nodal loads, respectively, the force vector can then be obtained by assembling all global member force vectors and the displacement vector can directly be related to all global member displacement vectors ; that is, By using the fact that is a function of for all along with the relations (23)-(24), the force vector is therefore a function of the displacement vector ; that is, It should be remarked that, for any prescribed displacement vector , the force vector can be determined by applying the following equations consecutively: (24), (4), (5), the relation developed in Section 3.1, (8), (7), and (23). The best linear approximation of the relation (25) in the neighborhood of can also be obtained as where is the structure tangent stiffness matrix evaluated at . By using the relations (23) and (24) and the fact that is available for all members through the relations (21) and (22), can therefore be obtained from a standard assembly procedure of all .
5. Solution Procedure
It should be remarked that the nonlinearity of the relationship between the force vector and the displacement vector stems directly from two main sources, the material nonlinearity present in the moment-curvature model (3) and the nonlinearity associated with the superposed large rigid body rotation in the total displacement decomposition scheme.
Although the procedure for determining the force vector for a given displacement vector is clearly stated in the previous section, it still requires solving the force vector for a given displacement vector in the corotational coordinate system for all members. This nontrivial task can be achieved as follows: (i) the force is computed directly from (10); (ii) the end moments are obtained for any prescribed by solving nonlinear equations (13) using Newton-Raphson scheme; and (iii) the reactive forces are obtained from (9). It is noted that the linearized relation (14) along with (15) is useful for the iterative procedure carried out in step (ii) and the algorithm indicated above is restricted only to the constitutive relation (3).
A more challenging and more practical task is associated with the determination of the deformed state of a structure under a specified set of applied loads. Since the relationship cannot readily be inversed, the solution for the displacement vector for a specified force vector in each load step is obtained numerically via Newton-Raphson iteration. Details of the iterative procedure are given below.(a)Choose an initial guess of the displacement vector, ().(b)Obtain , , and from (24), (4), and (5), respectively, for all .(c)Obtain for each using procedure indicated above and then compute for all using procedure described in Section 3.1.(d)Obtain from (8) and from (19) for all .(e)Obtain from (7) and from (22) for all .(f)Obtain and from direct assembly procedure.(g)Check convergence condition where tol is a specified tolerance. If such condition is satisfied, iteration is terminated and becomes an approximate solution; otherwise, proceed to the next step.(h)Update the displacement vector from and then return to step (b).
It is remarked that, to accelerate the iterative procedure for the current load step, the information from the previous load step is utilized to provide the suitable initial guess. In addition, the nonuniform step size is employed to obtain results for the given loading history and this proves to enhance the computational efficiency of the proposed technique especially when the structure undergoes very large displacement and rotation, and the bending moment within members is far beyond the linear regime.
6. Adaptive Remeshing Algorithm
To ensure that the small-rotation-based kinematics employed in the formulation is suitable for modeling the relative displacement of any member in the corotational coordinate system , meshes used in the analysis must be sufficiently refined. In general, segments or regions whose curvature is significantly large must be discretized into several members whereas those possessing small curvature or being nearly straight require less number of members in the discretization. To achieve this crucial task efficiently and conveniently, a simple automatic remeshing scheme is integrated into the implemented numerical procedure as briefly described below.
For the current load step, the analysis is first carried out using the mesh from the previous load step (the initial mesh is required and must be supplied for the first load step). Then, the maximum relative rotation of all members (i.e., maximum rotation measured in system) is calculated by using the curvature area equations and then is checked with the specified rotation limit. All members whose maximum relative rotation exceeds the given rotation limit are subdivided into two members by adding a new node at their center as shown in Figure 4. Once the new mesh is generated, the analysis is carried out again for this load step to obtain the new equilibrium state (following the procedure described in Section 5). The process of checking the requirement for member subdivision and performing the analysis of the current load step with the new mesh is repeated until the maximum relative rotation of all members in the current mesh is less than the rotation tolerance. It should be remarked that the information of all new nodes and members can be obtained by using results of old members before subdivision. This implemented adaptivity not only reduces effort for mesh generation but also provides a means to obtain a suitable mesh where the refinement is carried out only for regions needed.
7. Numerical Results
To examine both the accuracy and convergence behavior of numerical solutions generated by the proposed technique, various structures under different loading conditions are fully investigated. A simple problem whose solution exists analytically is considered first not only to verify the formulation and implementation but also to demonstrate the relationship between solution error and the level of mesh refinement. The latter should provide the guideline for setting the rotation limit in the adaptive remeshing scheme to control the solution error. In addition, for problems whose analytical solution does not exist, benchmark results obtained from a reliable finite element package are utilized in the verification procedure. Note in particular that all the benchmark solutions used in the verification procedure are obtained based on the same constitutive relation (3) while the geometric nonlinearity is handled by exact kinematics for large displacement and rotation (i.e., the approximations and are not removed). Once the implemented algorithm is tested, it is then applied to analyze more complex and practical structures and a selected set of results is reported and discussed. In the numerical study, the following nonlinear moment-curvature relationship is utilized: where are model constants. It should be noted that and denote the bending moment and curvature at the boundary of a linear regime and are parameters used to describe the nonlinearity of the model. Since the relation (27) contains both linear and nonlinear regimes with the nonlinearity being controlled by , it can be used to model a wide range of moment-curvature relationship of the beam cross section. When are taken to be , , and , the model recovers the linear, bilinear hardening, and bilinear without hardening relationship, respectively (see Figure 5(c)). Graphs of relationship of the family (27) for various values of are shown in Figure 5; clearly, as and increase, the slope of the nonlinear part is descending.
(a)
(b)
(c)
7.1. Cantilever Beam Subjected to End Moment
Consider a cantilever beam of length and cross-sectional area . The beam is made of a homogeneous material with , , , and and it is subjected to a concentrated, counterclockwise moment at the tip as shown in Figure 6(a). First, the analysis is carried out without the automatic remeshing algorithm and five meshes containing 1, 2, 4, 8, and 16 members shown in Figure 6(b) are utilized.
(a)
(b)
To demonstrate the capability of the current technique to analyze structures undergoing large displacement and rotation, values of the applied moment are chosen to be sufficiently large in the numerical experiments. The computed vertical and horizontal displacements at the tip of the beam from all five meshes are normalized by the length and then reported in Table 1 for along with the maximum relative rotation (measured in the corotational coordinate system). The corresponding deformed shapes are shown in Figure 7. It is evident that the numerical results converge to the exact solution as the mesh is refined. This is due to the fact that as the number of members in the mesh increases, the rotation in each member measured in the corotational coordinate system becomes smaller and this, as a result, improves the suitability of using small-rotation-based governing equations (1) and (2). Figure 8 shows quantitatively the relationship between the percent error of both the horizontal and vertical displacements at the tip of the beam and the maximum rotation of all members measured in the corotational coordinate system as the mesh is refined. Since the maximum rotation is a one-to-one correspondence to the number or length of discretized members, results indicated in Figure 8 can be used to obtain the rotation limit to achieve a required level of accuracy of an approximate solution.
To ensure that the maximum rotation measured in the corotational coordinate system of each member is sufficiently small (i.e., not exceeding a specified rotation limit), the automatic remeshing scheme described in Section 6 is integrated in the analysis procedure. Comparison of results with the exact solution for the same applied load level (i.e., ), when such adaptive strategy is utilized, is shown in Table 2 for different values of the rotation limit. According to this set of results, it can be seen that as the rotation limit decreases, the number of members, together with the level of accuracy, increases. Deformed shapes of the beam for and different rotation limits are also reported in Figure 9. The number of members in the final mesh for the given rotation limit and load level is also provided. It is apparent that numerical results, for this particular case, are almost indistinguishable from the exact solution although the large rotation limit (i.e., 20 degree) is used in the adaptive scheme. Note in addition to the fact that the number of members in the final mesh resulting from the adaptive scheme also depends on the curvature of the member; for this particular problem, as the applied end moment increases, the curvature of the beam increases and, as a result, the mesh must be more refined to ensure that the relative rotation does not exceed the rotation limit.
7.2. Diamond-Shaped Frame Subjected to Vertical Load
Consider a diamond-shaped frame with a fixed support at the bottom corner and subjected to a vertical load at the top corner as shown schematically in Figure 10(a). All four sides of the frame have the same length and the same properties, that is, with , , , and . Note in particular that the moment-curvature relationship chosen in the numerical study is bilinear. In the analysis, three different load levels associated with are considered, the initial mesh consisting of four members as shown in Figure 10(b) is utilized, and the rotation limit is set equal to 2 degrees to ensure the high level of accuracy.
(a)
(b)
The deformed shapes obtained from the current technique are reported in Figure 11 along with the benchmark solution generated by a reliable finite element analysis package. The final mesh for each load level resulting from the adaptive scheme is shown in Figure 12. It is evident that as the load increases, the mesh is more refined. In particular, most of members are located in the corner region since the corresponding curvature is relatively large. In addition, results shown in Figure 11 indicate that numerical solutions obtained from the proposed technique are highly accurate and almost identical to the benchmark solutions.
(a)
(b)
(c)
Note also that the deformed shape of the frame can still be accurately captured if the rotation limit changes from 2 to 20 degrees as shown in Figure 13. In particular, the percent difference of the vertical displacement at the top corner, the horizontal displacement at point (indicated in Figure 10), and the vertical displacement at point obtained by using two different rotational limits (i.e., 2 and 20 degrees) are 0.323%, 1.235%, and 0.322%, respectively.
7.3. Cantilever Beam Subjected to End Loads
To demonstrate the influence of large displacement and rotation on values of internal forces such as the shear force, axial force, and bending moment at any cross section, let us consider a cantilever beam of length and cross-sectional area , made of a material with , , , and . The beam is subjected to a horizontal force , a vertical force , and a moment at the tip as shown in Figure 14(a). In the analysis, an automatic remeshing scheme is used with the rotation limit equal to 2 degrees and an initial mesh containing only one member (see Figure 14(b)).
(a)
(b)
In the numerical study, the end loads are applied proportionally such that , where is a load parameter. The deformed shapes obtained from the analysis for , and 50 are reported in Figure 15. It is remarked that three different values of the loading parameter are chosen to represent the beam undergoing small, intermediate, and large deflections, respectively. The axial force, the shear force, and the bending moment at any cross section of the beam for above proportional loads are obtained and reported along with those generated by the linear elastic analysis in Figures 16, 17, and 18. It is evident that results obtained from the proposed technique show significant deviation from the linear solution and such difference increases as the magnitude of applied loads becomes larger. This is mainly due to that the deformed configuration of the member at larger loads is significantly different from the undeformed configuration and this therefore boosts the influence of the axial-shear-bending interaction. In particular, the axial force and shear force obtained from the linear analysis are constant throughout the member (due to that the equilibrium is enforced in the undeformed state) whereas those obtained from the large displacement and rotation analysis vary nonlinearly along the member due to the significant change of the member axis. This similar trend was also observed in the case of the bending moment. However, the deviation from the linear elastic solution for the bending moment is less significant since the influence of the -delta effect is reduced by the considerable reduction of the moment arm of the force when the beam undergoes very large displacement and rotation. Clearly, the increase in moment due to the -delta effect is less than the decrease in moment due to the reduction of the moment arm for and vice versa for . It can be noted also that as increases, the bending stiffness of the beam is further reduced since the region where the bending moment is larger than almost occupies the entire beam.
7.4. Portal Frame Subjected to Horizontal and Vertical Forces
As a final example, let us consider a problem of finding the deformed shape of a portal frame clamped at its bases and subjected to horizontal and vertical loads as shown in Figure 19(a). The given frame is made of a material with Young’s modulus and moment-curvature relationship described by (27) and all members have the same cross-sectional area . An initial mesh containing only three members, indicated in Figure 19(b), is adopted in the analysis and the rotation limit equal to 2 degrees is utilized in the automatic remeshing scheme. To demonstrate the influence of material nonlinearity on the deformed shape of the frame, various values of material parameters are considered along with and .
(a)
(b)
The normalized computed deformed shapes of the frame are reported in Figures 20 and 21 for , , and various material parameters . The load level considered in the analysis is chosen to ensure that the total displacement and rotation of the frame are sufficiently large to demonstrate the capability of the proposed technique to capture highly geometric and material nonlinearities. It is evident from this set of results that the frame becomes softer as increases for a fixed and, similarly, as increases for a fixed . This is clearly consistent with the moment-curvature relationship; that is, the larger the value of and , the higher the material nonlinearity and lesser material stiffness.
8. Conclusions and Remarks
The numerical technique which is capable of modeling two-dimensional, nonlinearly elastic, extensible frames undergoing large displacement and rotation has been developed. The force-deformation relationship for all members in the corotational coordinate system has been established based on the small-rotation-based kinematics. The large displacement and rotation feature has been taken into account in the analysis by superposing the displacement measured relative to the corotational coordinate system by the large rigid body displacement. The material nonlinearity has been integrated in the modeling via the prescribed arbitrary nonlinear moment-curvature relationship. The final system of nonlinear algebraic equations has been solved by Newton-Raphson iterative scheme. To ensure the validity of using the small-rotation-based governing equations in the corotational system, an adaptive remeshing algorithm has been implemented to automatically obtain an optimal mesh for a specified rotation limit. The benefit of this proposed technique is to avoid the direct treatment of highly nonlinear differential equations in the development of force-displacement relationship and the corresponding element tangent stiffness matrix.
From extensive numerical experiments, it has been found that the proposed technique yields very accurate numerical results compared with the benchmark solutions provided that the mesh used in the analysis is sufficiently refined. The adaptive remeshing scheme is proved to provide a convenient means to obtain an optimal mesh used to acquire accurate results; only meshes of regions or segments whose curvature is sufficiently large are refined. It is worth noting that the present study is still restricted to a simple model for the treatment of material nonlinearity. The moment-curvature relationship utilized is free of the axial interaction and applies only to history-independent materials. The former can potentially lead to an inaccurate prediction when applied to axially dominated structures whereas the latter renders the technique incapable of modeling inelastic structural responses. One possible extension of the current work to better model material nonlinearity is to incorporate a fiber model to directly construct the cross-sectional behavior from a basic constitutive model in the level of material points. This additionally eases the treatment of the axial-bending coupling behavior.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
Authors gratefully acknowledge financial supports provided by Stimulus Package 2 (SP2) of Ministry of Education under the theme of Green Engineering for Green Society.