Numerical fundamentals and interactive computer graphics system for the nonlinear analysis of planar frames

Recent scientific and computational advances have facilitated the analysis of slender structural systems subject to instability. With the employment of more sophisticated numerical tools and algorithms, it is possible to accurately determine the critical points (limit and bifurcation loads) as well as the post-critical behavior of the structural system. In the computational context, efficient data structures are needed to enable code and graphic interface expansion for the generation of models and visualization of the results obtained. Thus, an interactive object-orientated graphic computational system is presented herein. It has been developed using MATLAB/GUI, with pre-processing, analysis and post-processing capacities for planar structural frames. The nonlinear finite element developed and implemented for the structure modeling is formulated considering second order effects. Therefore, with the computational tool presented, the geometric nonlinear effects and stability of the structural system can be directly addressed, and the visualization of the numerical results are accessed through interactive controls that permit data inclusion and analysis verification. The engineerdesigner can see the structural model discretization, the equilibrium path, its deformation configuration, the force and bending moment diagrams at the moment that he runs the program and in each load step. It is also possible to export the images, videos or tables of the obtained numerical results. The example presented demonstrates the capacity of the developed graphic computational system.


Introduction
The numerical analysis of structural systems is of computational nature due to the great number of operations needed for solution of the problem.This characteristic is extremely accentuated when nonlinear effects are considered in conjunction with the application of the finite element method (FEM).Therefore, the development of computational tools is of great importance for the practice of structural calculations and designs.
As the computational models are refined to consider the various characteristics present in the actual structures, i.e. material plasticity, second order effects, and semi-rigid connections, the numerical solutions reached are closer to the actual behavior of the structure (Silva, 2009).This also implies that the computational implementations have also become more complex, together with a considerable increase in processing time and storage cost.Thus, matters such as the program's paradigm choice (code organization strategy), efficiency, access, and memory usage become relevant and can determine the feasibility of the utilization of a given computational system.
The object-orientated programming paradigm (OOP; Stroustrup, 2013) has become a tendency for developing complex and/or large-scale computational systems, since it offers great advantages for maintenance, productivity and higher-quality software, as well as for code usage and extensions when compared to the conventional approach for programming struc-tures.Thus, OOP has herein been used throughout the complete development of the computational system, from the structural problem solution to the graphic interface for insertion in the model, as well as for determining the analysis options and visualizing the results.
Structural modeling requires the insertion of a great number of parameters, such as nodal coordinates, external loads, cross section properties, materials, among others.In addition, there are various options for the analysis and visualization processes.Therefore, the use of iterative graphic interfaces is of great assistance to the analyst/engineer, since it permits better visualization, control and understanding of the process when compared to the use of standard data files.
This article presents an interactive object-orientated graphic computational system, developed using MATLAB/GUI (Kwon and Bang, 1997;MATLAB, 2015), for the geometric nonlinear analysis of planar structural frames.This computational tool is enabled with pre-processing, analysis and post-processing capacities, and the nonlinear finite element developed and implemented for the structure modeling is formulated considering large displacements (second order effects).Therefore, the geometric nonlinear effects and stability of the structural system can be directly addressed, and the visualization of the numerical results are accessed through interactive controls that permit data inclusion and analysis verification.The analyst/engineer can see the structural model discretization, the equilibrium path, its deformation configuration, the force and bending moment diagrams at the moment that he runs the program and in each load step.It is also possible to export the images, videos or tables of the obtained numerical results.
This article is organized as follows: Sections 2 and 3 show the finite element formulations according to the total Lagrangian and corotational reference systems, which deal with the non-linearity associated with the second order effects.Also presented are the internal force vector and the stiffness matrix of these formulations.Section 4 discusses the nonlinear solver strategy adopted, and Section 5 presents the interactive graphic computational system.The Section 6 brings a numerical application of the graphical tool and finally, some conclusions about the computational system are presented.

Total lagragian formulation
In the plane frame formulations, it is assumed that the cross sections of the bar undergo rigid body movement (Reissner, 1972).Thus, the movement of a point on the cross section is composed of a horizontal translation u, a vertical translation v and a rotation θ around the bending axis.In the following formulations, it is considered that the axis of the element coincides with the structural axis.For a generic case of an arbitrarily positioned element in the plane, the transformation of the quantities must be performed, just as in linear analysis.
In addition, (' ) denotes the derivatives in relation to the spatial variable along the axis of the element x 0 .
The displacement vector u of a point at distance y 0 of a neutral line of the cross section of the element is given by (Reissner, 1972): The longitudinal ε xx and transversal ε xy components of the Green-Lagrange deformation tensor are: ε xx = ε -y 0 χ and ε xy = g in which the axial ε, shear g and bending χ deformations, which form the section deformation vector ε = [ε g χ] T , are defined according to: The virtual deformations are found from:

u u
with the matrix B defined as the deformation-displacement matrix that incorporates the nonlinear geometric effects.It is as-sumed that the deformations in the material are small, and as such, the St. Venant material hypothesis is applied.Furthermore, it is assumed that the material behaves elastically in such a way that the normal stress σ xx and the shearing stress can be expressed as: where E and G are the longitudinal and shear modulus of elasticity, respectively.
The resulting forces in a cross section are obtained as follows: where , N, S and M are the normal force, the shearing force, and the bending moment, respectively.In the model adopted, the transversal deformation is assumed to be constant in the section.
To diminish the amount of errors generated with this hypothesis, a Timoshenko correction coefficient κ is introduced in the integration of the transversal stress of the section (Crisfield, 1991).
Applying FEM, the displacements are interpolated by linear functions, through the nodes values, as shown below: Ouro Preto, 72(2), 199-207, apr. jun. | 2019   0 where L 0 is the element's initial length, while the element internal forces vector is obtained according to (Zienkiewicz and Taylor, 1991):

N S M dx
The tangent stiffness matrix of the element is calculated from the variation of the internal force vector with respect to the displacements, and is given by: [ ] where the first one represents the material stiffness matrix regarding the deformation increments undergone by the element; the second one represents the geometric stiff-ness matrix regarding the changes of the forces acting on the element.

Corotational formulation
The movement undergone by an element can be divided into two components: one due to the movements of the rigid body, where the element is translated and rotated as a whole, without its configuration being altered; the other parcel is related to the movements that generate deformations (Crisfield, 1991).The second type of movement is of greater interest, since it generates the forces in the structure.
Thus, to filter the rigid body movement from the total movement of the element, a corotational referential is introduced, as in Figure 1.As this referential follows the element, only the movements that cause deformation are detected.In this corotational system, the non-zero nodal displacements are given by: where u j is the longitudinal displacement, θ i and θ j are the rotational dis-placements; the rigid body rotation is θ r and the length of the deformed configuration L are defined according to: In the corotational formulation, it is important to relate the local system with the global one in order to obtain the internal forces vector F i and the tangent stiffness matrix K in the global system.This is achieved by taking the variation of the corotational displacements with respect to the global ones, as follows: with a, b, c i and c j given by: [ ] [ ] Numerical fundamentals and interactive computer graphics system for the nonlinear analysis of planar frames The internal forces vector F i and the stiffness matrix K in the global system can be related to the same in the local system following the expressions (Crisfield, 1991): As the movements that cause deformations are considered small, it is assumed that the behavior of the element is linear in the local system.Thus, in the local system, the Timoshenko (K T ) and Euler-Bernoulli (K EB ) element stiffness matrices are given by:

Solution of the nonlinear equations
With the FE discretization, a nonlinear algebraic equations system is obtained and can be written as: where u is the system's displacement nodal vector, F r is the reference external forces vector, l is the load parameter, F i is the internal forces vector and r is the residual forces vector.For solution of the nonlinear system, Eq. ( 17), an incremental process applied together with an iterative method, generally the Newton-Raphson one, where given an equilibrium state (u n , l n ) and a predicted solution (Du 0 n+1 , Dl 0 n+1 ) for the new state of equilibrium, a new estimate is obtained from the sub-increments δu k+1 n+1 and δl k+1 n+1 in iteration k according to: The sub-increment δu k+1 n+1 can be calculated using the expression below: ( ) ( )

Iteration strategies
To determine the sub-increment in Eq. 19 and the increments in Eq. 18, it is necessary to determine the subincrement of the load parameter δl k+1 n+1 .The choice of this load parameter is done in such a way that the numerical process is able to pass through load and displacement limit points, as well as bifurcation ones (Pires, 2012).Among the various iteration strategies implemented are: linear, cylindrical and spherical arc-length (Crisfield, 1991); minimum residual displacement norm (Chan, 1988); orthogonal residual procedure (Krenk, 2009).

Load increment strategies
Once the new estimate of the increments is obtained, the system's equilibrium is tested.Various convergence criteria can be used (Crisfield, 1991).
Generally, the force criteria is used, where the norm of r should be less than a previously stipulated tolerance.If the criteria is satisfied, a new equilibrium point is established and a new estimate for the displacement increments Du 0 n+1 becomes necessary.In general, this estimate is obtained as presented below: Various strategies can be used to determine this incremental estimate of the initial load parameter Dl 0 n+1 .Among those implemented are: direct load increment (Krenk, 2009); cylindrical and spherical arc-length (Crisfield, 1991); and constant external work (Silva, 2009).Figure 2 shows, in flowchart form, the nonlinear solving strategy adopted in this study.

Figure 2
The incremental-iterative strategy for solving the nonlinear structural equations.

Interactive graphic computational system
In this study, an interactive computer graphics with capacity for pre-processing, analyzing, and post-processing was developed using MATLAB/GUI (MATLAB, 2015), and employing the paradigm designated object-orientated programming (OOP), in order to help elaborated designs for planar frame structures under large displacements.

Object-orientated programming (OOP)
Object orientation is a programming paradigm based on the composition and interaction of various units of the program called "objects".In this programming paradigm, a set of classes defines the objects in the computational system.Each class determines the behavior and attributes of its objects (nodes, material, sections, elements), as well as their relationship with other objects (Farrell, 2009).Herein, the OOP paradigm was used because of its advantages, such as organization of the generated codes, inheritance, and polymorphism.

Pre-processing
The opening screen of the developed computational system is shown in Figure 3.With the embedded preprocessing capacity, the user can insert all of the relevant parameters for the idealized structural model.As the structure is modified in the program, the structural design is updated in the principal canvas.Thus, the geometry, together with the boundary conditions of the problem, can be visualized and exported.

Figure 3
Initial screen of the interactive graphics system.

Numerical fundamentals and interactive computer graphics system for the nonlinear analysis of planar frames
Figure 4 shows the menus available in the program.In the first group of menus, the user can administrate the database of the analysis; begin a new project; re-open an existing project; or save a current project.Besides these options, the user can obtain assistance about a given task.For more details, see Santana (2015).In the second group of menus, the model parameters are inserted.In other words, the nodes, material, cross section of the bars, elements, support conditions and loads can be introduced.Once a determined menu is selected, a panel referent to the parameters that should be inserted is generated on the left side of the screen, as shown in Figure 5.In the nodes panel, besides the coordinates, the boundary conditions of each node can be inserted, such as loads, supports and prescribed displacements.
Various material and sections types can be chosen in their respective interactive panels, including the concrete described in the NBR 6118 ( 2003) and the steel described in the NBR 8800 ( 2008) standards, as well as common sections encountered in practice.

Structural analysis
In the third group of menus presented in Figure 4, the referent parameters for the structural analysis are generated, as illustrated in Figure 6.In this panel, options such as iteration strategy, load increment strategy, tolerance and number of desired iterations can be specified.This panel also has a "solve" key that initiates the solution process for the structural model.The numerical fundamentals presented in Sections 2, 3 and 4 are the bases of the Structural Analysis option.

Post-processing
With the last group of menus (Figure 4), the results obtained in the analysis can be visualized by means of graphics, animations and tables.In the first sub-menu of this group, it is possible to visualize the deformed configuration; in the second, the equilibrium path; in the third, the internal forces diagram; and in the fourth, the previous results combined, as can be seen in Figure 7.A control panel of the load step is also generated for each of the results visualized.In the fifth submenu of this group, the results can be visualized as tables.These tables can be exported as database files.

Application example
The developed computational system was applied to analyze a pinned supported circular arch with a diameter of L = 100, under a centrally concentrated (Fig. 8a) and eccentric (Fig. 8b) load.This structural system was initially studied by Harrison (1978) and later by Yang and Kuo (1994) and Galvão (2000).The arch was discretized with twenty-six elements, has material with elastic modulus E = 2000, cross section area A = 10 and inertia I = 1.The eccentric load is induced by means of a moment M = 0.001P, applied at the central node.The results obtained were compared with those present in Galvão (2000).The models were analyzed using the Total Lagrangian and Corotational beam formulations according to the Timoshenko theory, and as the nonlinear solution strategy, the cylindrical arc-length was adopted in the iterative process and in the load parameter increment strategy.
Figures 9a and 9b present the circular arch equilibrium paths for the two loading cases (central load, perfect system; and eccentric load, imperfect system), where the variation of the vertical displacement of the arch's central node with respect to the load factor can be observed.See that these two load configurations provoke high nonlinear responses in the arch, characterized by the presence of various load and displacement limit points.However, in the case of the imperfect system, the stiffness of the structure is considerably less, characterizing the unstable behavior of the structure.The loops in the equilibrium paths also demonstrate the existence of various equilibrium configurations for a given value of the load parameter.These curves were obtained by considering load P to be equal to 0.4N to initiate the analysis and the standard Newton-Raphson method was applied.The results obtained with the formulations implemented are in accordance with those of Galvão (2000), who adopted the classical beam theory.Figures 9a and 9b present the circular arch equilibrium paths for the two loading cases (central load, perfect system; and eccentric load, imperfect system), where the variation of the vertical displacement of the arch's central node with respect to the load factor can be observed.See that these two load configurations provoke high nonlinear responses in the arch, characterized by the presence of various load and displacement limit points.However, in the case of the imperfect system, the stiffness of the structure is considerably less, characterizing the unstable behavior of the structure.The loops in the equilibrium paths also demonstrate the existence of various equilibrium configurations for a given value of the load parameter.These curves were obtained by considering load P to be equal to 0.4N to initiate the analysis and the standard Newton-Raphson method was applied.The results obtained with the formulations implemented are in accordance with those of Galvão ( 2000), who adopted the classical beam theory.Finely, as main advantages of the numerical modeling system, there can be highlighted: the ease and straight modeling process (user interactive controls which allows for the inclusion and removal of data); the user different options of nonlinear solvers (analyzes); the resulting visualization, such as the nonlinear equilibrium paths; and it can also serve as an educational tool.The limitations of the tool are: no consideration of the inelastic and semirigid connection effects; it is based on plane truss and beam-column elements and just permits 2D visualization not addressing nonlinear transient structural problems.

Concluding remarks
This work presented an interactive object-orientated graphic computational system for second order analysis of planar structural frames.The computational system developed with MATLAB/GUI proved to be very useful in the modeling process and evaluation of the numerical results.With pre-processing, analysis and post-processing capacities, and nonlinear finite element formulations implemented, this computational tool offers to the analyst/engineer a direct access to the geometric nonlinear effects and stability of the structural system.The visualization of the numerical results is accessed through interactive controls that permit data inclusion and analysis verification.The structural model discretization, equilibrium paths, deformation configurations and resultant force diagrams can easily be seen at the moment that the analyst runs the program and in each load step.
Additionally, the computational program developed here can be used as an educational graphic tool for undergraduate and graduate students from architecture, civil and mechanical engineering courses.
The choice of the OOP paradigm, as previously emphasized, facilitated the organization, maintenance and expansion of the computational system's code.The Total Lagrangian and Corotational beam formulations have shown to be adequate in the analysis of planar frame systems subject to instability.
Other numerical structural analyses using the developed computational system are presented in Santana (2015).Future research will incorporate the issues of ma-terial inelasticity, semi-rigid connections, and their dynamic effects in the structural analysis (Batelo, 2014).

Figure 5
Figure 5Interactive panels of data input.

Figure 6
Figure 6Analysis and results panels.

Figure 7
Figure 7 Visualization of the numerical results.

Figure 8
Figure 8 Pinned circular arch: FE discretization and loading conditions.
Figure 9 Equilibrium paths of the pinned circular arch.
Figure 10 Deformed configurations of the pinned circular arch.
Figure 11 Bending moment diagram of the pinned circular arch.
Received: 31 August 2017 -Resubmitted: 21 August 2018 -Accepted: 2 December 2018.All content of the journal, except where identified, is licensed under a Creative Commons attribution-type BY.