INFLUENCE OF MODE INTERACTIONS ON INELASTIC BUCKLING EFFECTIVE STRESSES IN THIN-WALLED COLUMNS USING FINITE STRIP METHOD

D.D. Milašinović1 P. Marić2 Ž. Živanov2 M. Hajduković2 UDK: 624.042.3.072.31 DOI: 10.14415/konferencijaGFS2019.036 Summary: The inelastic buckling, free vibration and damage of prismatic thin-walled columns, considering two sources of non-linearity (geometrical non-linearity due to large deflection and material non-linearity due to inelastic behaviour), are analysed by implementing the semi-analytical finite strip method (FSM). A very basic continuum damage model with one damage parameter is implemented in conjunction with a mathematical material modelling approach named rheological-dynamical analogy (RDA) in order to address stiffness reduction due to inelastic behaviour. The RDA is a type of inelastic analysis, which transforms one category of very complicated material non-linear problems to simpler linear dynamical problems using modal analysis. The effective stress, as defined in the damage mechanics, is used in this paper. An analysis of damaged composite thin-walled wide-flange columns made of a polymer resin, which is reinforced with fibrous material, is carried out. Numerical examples computed by application of extensive hybrid parallelization present the influence of mode interactions on the damage variables and effective stresses in thin-walled columns.


INTRODUCTION
The thin-walled structures are structures which are generally made by joining flat plates at their edges. An important sub-set of these structures, which are the main concern of this paper, are essentially prismatic forms that can have some stiffeners, such as those used in column members, stiffened slabs and box girders. The analysis of the behaviour of these structures is approached using the FSM. The FSM is based on the basis functions (or eigenfunctions), which are derived from the solution of the beam differential equation of transverse vibration, and proved to be an efficient tool for 1. analysing a great deal of structures for which both geometry and material properties can be considered as constant along a main direction. This method was pioneered by Cheung [1], who combined the plane elasticity and the Kirchhoff plate theory. Wang and Dawe [2] have applied the elastic geometrically nonlinear FSM to the large deflection and post-overall-buckling analysis of diaphragm-supported plate structures. Current versions of FSM are capable of modelling the plates with different thickness and material properties along the longitudinal direction of the strip. Naderian at al. proposed an integrated finite strip discretization scheme for the flat shell spline finite strip for modelling of long-span cable-stayed bridges [3]. This approach of FSM is a very efficient technique and its application to the dynamic inelastic buckling of thin-walled structures is highly appreciated. Kwon and Hancock [4] developed the spline FSM to handle local, distortional and overall buckling modes in post-buckling range. Interaction of two types of column buckling (failure) in thin-walled structures, local and global (Euler) column buckling, may generate an unstable coupled mode, rendering the structure highly sensitive to imperfections. The geometrically nonlinear harmonic coupled FSM (HCFSM) [5,6] is also one of the many procedures that can be applied to analyse the large deflection and buckling-mode interaction of thin-walled structures. An analysis of the buckling-mode interaction is carried out in [7], taking into account the viscoelastic behaviour of material. Furthermore, to address these issues, strips with nonuniform characteristics in the longitudinal direction and transverse stiffeners have been used in geometrical non-linear static analysis of prismatic shells [8]. A very important effort for improvement of this approach is done via parallelization. Application of parallelization was necessary to overpower the increased computational complexity, caused by coupling of all series terms in the HCFSM formulation, by the usage of the high number of modes, as well as by the need to calculate a large number of different integrals. It was inevitable to combine different parallelization models -MPI and OpenMP and use different sources of parallel resources (multicomputer, multiprocessormulticore, GPGPU based, or cloud available environment). Our results show that rational usage of hybrid parallelization approach to speed up the critical parts of our code contributes in significant performance improvements [9,10]. If uniformly compressed thin-walled structures undergo inelastic deformations, these structures generally sustain two sources of non-linearity (geometrical non-linearity due to large deflection and material non-linearity due to inelastic behaviour). Due to the slender nature of the cross-sections, their behaviour is inevitably complex, with several parallel buckling phenomena influencing performance and limit states. The analysis presented in this paper is based on the RDA [11,12]. The RDA is a type of inelastic analysis, which transforms one category of very complicated material non-linear problems to simpler linear dynamic problems by using modal analysis [13]. In this paper, we present a new approach in which the damages are investigated by RDA to solve the inelastic buckling effective stresses in thin-walled structures. The uniaxial RDA dynamic modulus is based on a concept of a complex modulus of viscoelastoplastic (VEP) materials. This modulus is used to obtain one simple continuous modulus function and a stress-strain curves which were verified with experiments. The key global parameters, such as the creep coefficient, Poisson's ratio and the damage variable, are functionally related [14]. For the analysis of thin-walled structures using the FSM, an inelastic isotropic 2D constitutive matrix is derived starting from the uniaxial state of stress [7]. The damage identification method requires the use of the FSM to model the structure in its undamaged state as well as information on its dynamic properties when in damaged state. Finally, the effective stress, as defined in the damage mechanics [15], is used. For the more complex series part of displacement function Ym, determined by the end conditions (i.e. both ends clamped), it's not possible to find the exact analytical solution of eigenfunctions [16,17], due to the hyperbolic sine and cosine functions which are embedded within their characteristic equations. In [18] it's shown that with each increasing mode the characteristic equation root-finding error grows exponentially for any non-trivial edge boundary condition, and once the higher modes have been reached the accuracy is lost. It has also been demonstrated that the said large root-finding errors may lead to severe accuracy issues when computing basis functions and their integrals, especially when higher modes are involved. A hybrid method [18] has been proposed that works around these issues by coupling multiple approaches to minimize the numerical errors throughout the entire computation pipeline. The idea underlying the approaches presented in this paper is that the influence of mode interactions appears if the dynamic properties in undamaged state as well as in damaged state are used together. The mode interactions have great influence on damage variables and inelastic buckling effective stresses in thin-walled columns.

Buckling
The non-linear strain-displacement relations in FSM can be predicted by combination of plane elasticity and the Kirchhoff bending plate theory. This has been accomplished in [5], by using the second-order terms of Green-Lagrange strains. However, since longitudinal loading is assumed here (see Fig. 1), the second-order terms are only necessary for the longitudinal normal strain where u0 and v0 are, respectively, displacements in the middle surface in x and y directions, and w is displacement in z direction. In FSM, which combines elements of the classical Ritz method and the finite element method (FEM), the general form of the displacement function, can be written as a product of polynomials and trigonometric functions where Ym(y) are functions from the Ritz method, Nk(x) are interpolation functions from FEM [1] and km q is a vector representing the m-th term nodal displacements. r is an integer specifying the number of series terms chosen for approximation and c represents the number of nodal lines of a strip. The most commonly used series are the basis functions which are derived from the solution of the beam vibration differential equation where a is length of beam (strip) and µ is a parameter [1]. The general solution of Eq.
(3) is with the coefficients C1, etc., to be determined by the boundary conditions. These have been worked out in the ref. [18] for the various boundary conditions but the basis function is listed below for a simply supported strip only We define the local Degrees Of Freedom (DOFs) as displacements u0, v0 and w, and the transverse slope amplitude w x ϕ = ∂ ∂ of a nodal line (DOFs=4), Fig. 1.
The total potential energy of a strip for plane stress and bending according to the linear theory is designated Π and is expressed with respect to the local DOFs [5] ( ) where U is the strain energy in which U p is the plane stress energy and U b is the bending energy. The conventional stiffness block matrices are, respectively  We introduce matrices, which are referred to as the strain matrices for the plane stress and bending, respectively Once the approximate functions for the plane stress Au are known it is possible to define the following matrices and vectors Hence, are the corresponding approximate functions for displacements u0, v0 and w in the middle surface, respectively, while u u q , v u q and w q represent vectors of displacement parameters in the nodal lines. The potential energy due to external surface loads pu and pw can be written simply as In order to obtain the equilibrium equations, the principle of minimum total potential energy is invoked Well known elements of the property matrices D p and D b for the plane stress and bending, respectively, for the orthotropic elastic material are y y x x y x x y xy x y x y x y x y In the above expressions, Ex, Ey, µx, µy, and G are the elastic constants while t is the thickness of the strip, Fig. 1.
Consider the simply supported flat shell strip. The strip is subjected to an initial stress σy, which varies linearly from side 1 to side 2 as shown in Fig. 1, but is constant along the longitudinal axis Considering the assigned stress distribution, from the non-linear strain tensor we include only the term given by Eq. (1). The total potential energy of a strip is defined as the sum of its strain energy, potential energy due to nodal line forces, as well as the additional potential energy due to the initial stresses. As far as linear stability is concerned, the nodal forces Q are zero and it's therefore possible to derive the eigenvalue equation [5] ( ) where is the corresponding geometric stiffness matrix [7], λ is the eigenvalue (the load factor is compression positive), and b is the eigenvector (buckling mode). Based on Eq. (16) for one finite strip we can form the eigenvalue equations for a system of finite strips (mesh). The eigenvalue problem is to extract the solution pairs im λ and bim for all DOFs i, and all series terms m = 1,...r. The buckling stresses can then be obtained as follows

Vibration
The stability and instability of a static equilibrium state are conventionally defined in terms of free motions of the system following an infinitesimal and once-and-for-all disturbance from the equilibrium state. The linear equations of free vibration in the matrix form are as follows [5]

393
M is the mass matrix and q  are accelerations in the nodal lines of a finite strip. We consider natural frequencies of vibration under an assumption that each DOF executes harmonic motion in phase with all other DOFs. Therefore This is the basic statement of the vibration problem. If there is neither geometric nor material non-linearity, then neither K nor M are a function of ω , and Eq. (20) is a standard eigenvalue problem. For a mesh of finite strips we obtain where im ω are the natural frequencies. For each im ω there is a corresponding vibration We find it very important, to point out that a duality exists between the free vibration and buckling under a uniform end compression, where the eigenvectors are identical. Thus, corresponding eigenvalues for a mesh of finite strips are related by where ρ is the mass density.

RDA -a short overview
This section can be thought of as an extension to Section 2.2. The purpose of developing a mathematical model for the rheological behaviour of solids is to permit realistic results to be obtained from mathematical analyses of damaged structures under various conditions, such as sinusoidal loading The governing differential equation between strain and stress that includes their first and second time derivatives has already been derived by the first author of this paper in [13] According to the basic RDA equations, the mass and stiffness are respectively The viscous damping is (the minimum one) that determines the return to equilibrium condition of a system without oscillations. Hence where * ϕ and ϕ are the VE and VEP creep coefficients, respectively. Consequently, the key result from the RDA complementary solution is the following relation between natural frequency and delay time T D , [11] 1 D T ω = . (25) with 0 σ being a constant and A σ a variable component of the cycle. σ ω is the stress frequency, Q σ ω ω = . Such stressing conditions have come to be known in rheology as dynamic loading and are used in the analysis of the fatigue behaviour [11]. In accordance with RDA, the fatigue strength under constant stress amplitude is is the ratio of minimum stress to maximum with the signs taken into account.
An approach based on energy considerations, which includes the rate of release of VEP energy, has already been adopted in the RDA theory. The VEP energy release rate is derived in [11] as follows According to the damage mechanics [22], the damage variable D may be characterized by the variation in Young's modulus E(D). Under an assumption that variation from undamaged to the damaged state is equal to the RDA dynamic modulus, E(D)=(1-D)EH=ER, the damage variable was obtained from the creep coefficient and relative frequency in [14] as follows Hence, the dependence between the damage variable and its associated VEP energy release rate determines the evolution of the damage according to the RDA, Fig. 3.  [11].

МЕЂУНАРОДНА КОНФЕРЕНЦИЈА
In this paper, the damage variable 0 under the quasi-static loading (δ→0) is used that was experimentally proved in [14] The assumption of isotropic damage is often sufficient to give a good prediction of the carrying capacity, the number of cycles or the time to local failure in structural components [22].

Effective stress
According to interpretations in the rheology, the viscoelastic Poison ratio is defined in a time or frequency domain [19]. Since it is a dimensionless value that grows in a defined domain to a defined boundary of incompressibility, it is associated with an increase in material damage. A similar property has a creep coefficient, which is also a dimensionless value and can be interpreted as a variable of damage. The RDA approach, which considers both the initial (undamaged) and damaged (viscoelastic) state of cylindrical rod for the analysis of the influence of Poisson ratio on the creep coefficient, has already been described by the first author in [11]. The energy equation results from application of the principle of conservation of energy to steady transfer of energy, where the equation of continuity, which results from the principle of conservation of mass is used, Fig. 4 ( ) (32) In this approach, based on Bernoulli's energy theorem and assuming that εE =σE/EH =0.001, the creep coefficient is expressed by the formula [11] ( )  . .
Hence, the mathematical expression for obtaining the value of the creep coefficient from the viscoelastic Poisson ratio is given by Poisson ratio, as defined in an idealized purely elastic material, is an elastic constant. However, such a material is hypothetical because it does not exist in the strict sense of the word. In any mechanical deformation the deformation energy is not only stored elastically, but part of it is invariably dissipated by viscous forces, in accordance with the second law of thermodynamics. This dissipation is responsible for the time dependence of the mechanical properties of any real material. Hence, the viscoelastic Poisson ratio in time domain chosen can be defined as follows [19]

399
Consequently, in order to solve the viscoelastic problem of structures, the dynamic modulus must be used in the conventional stiffness matrix K of a finite strip.
Bifurcation problem given by Eq. (16) belongs to the scope of linear analysis, and determination of critical load reduces to solution of eigenvalue problem. Beside that, Eq. (20) can be used to determine the natural frequencies and mode shapes of a structure, and it forms the basis for many damage detection methods that make use of modal properties in damage detection and identification [20]. The residual force is given as follows [21] ( ) where k, and m are, respectively, the stiffness and mass of the undamaged column. d ω is the natural frequency, which must be obtained from the viscoelastic eigenvalue problem. Hence, the damage variable used in the paper in the point of bifurcation is obtained from the natural frequencies as follows Finally, the effective buckling stress is given using the damage mechanics [22] where d σ is the inelastic buckling stress of damaged structure, which can be obtained either from viscoelastic eigenvalue problem, or from duality, Eq. (22).

Finite strip software for inelastic eigenvalue problem (MKTB)
The finite strip modelling software is composed of two loosely-coupled software modules (Fig. 6). The MKTB module is in charge of inelastic eigenvalue problem computation. It produces extensive quantity of numeric results that cannot be comprehended without an appropriate visualisation module.

Finite strip modelling supported by data visualisation
As the amount of output data from the MKTB module is overwhelming, to be able to look at required aspects of it, some aggregate representation is needed, and one of best ways to do it is to use visualisation. Depending of the wanted analysis, parts of the data need to be filtered by some criteria and displayed as 2D or 3D graph, or as multiple graphs on the same picture. For example, displaying 3D model of the beam requires extracting initial cross section layout and displacement data for each strip, all that for selected length and mode. Other types of graphs require, for example, extracting specific data for all available lengths and/or modes. It must be noted that our research requires almost constant adjustments and additions to data visualisation. To be able to respond quickly, our visualisation software needs to be flexible enough for on-demand adjustments to existing graphs, quick additions and adjustments of graph modification parameters, and a way to quickly expand it with new types of graphs. Visualisation is based on Qt framework and Python classes that are designed for rapid non-visual GUI development, making adding or removing elements a matter of one or few lines of code. The main window of visualisation software is divided into tabs, each representing a different view of the data. Every tab has two main parts, one larger that holds the graphics, and one smaller that contains numerous options for graphics customization (Fig. 7). Two helper classes for graphics representation are developed, currently both based on matplotlib (the actual back-end can be replaced relatively easily), one for 2D and one for 3D graphics. Adding a new graph type to the program consists of deriving a new class from main tab class, overriding of its plot method, to describe new data representation, and adding object of that class to the application window. Main tab class defines the main tab appearance (graphics part and options part), and has numerous methods for dynamic

401
addition of attributes and methods for their handling to the derived class. For example, if there needs to be an element to choose some values forms combo-box, all one needs to do is to call addLabelComboBox method which will dynamically add several attributes and methods to the derived tab class: • a label with text "File" to the options panel, followed by combo-box • an attribute that will hold data for currently selected string from combo-box • an attribute that will hold currently selected string from combo-box • an attribute that will hold all the Qt graphical elements that are added • a helper function that can be used to set what strings should be included as options for combo-box, and what optional data should be associated with each string (by default, data is the same as string)

Figure 7: Graphics customization.
All of the mentioned attributes and methods are available in derived tab class immediately after the call of addLabelComboBox. By default, after each change of the selected item, there will be a call to tab's update method, that in turn calls plot method. That way, there is no need for explicit update of the graphics display. Of course, custom callback functions can be defined for each element. Optional sync parameter for each added element can be used to synchronize its value with elements from other tabs that have the same name. This can be used, for example, to synchronize selected file on multiple tabs, so that each of the tabs will show its data based on the same source. The same traits described for combo-box are available for the rest of elements that can be added to options panel (radio and check buttons, numbers, text fields, sliders, etc.) making defining new tabs a relatively automated process.

Numerical applications
In this paper the simply supported ideally straight thin-walled wide-flange H-section column is analyzed. The column of length a=2310 mm consists of a web and two flanges of side b=152 mm [23]. The thickness t is 6.35 mm. The column is compressed axially. The bending stiffness EI=89.6 kNm 2 about the weak axis of the analysed column was computed from the information provided by the manufacturer for each section following the methodology developed by Barbero and Tomblin [24]. This information includes the type of fibres and matrix material, the local orientation of the fibres and the fibre content in the cross-section. Tab. 2 shows the flange and web bending stiffness components for 152 mm x 152 mm x 6.35 mm pultruded WF H-section column. The material properties in Table 3 for the implementation of HCFSM are obtained as detailed in [7]. These material parameters are also used for the FSM eigenvalue computations within this paper. The buckling behaviour of ideally straight column to axial loading and small transverse concentrated load has already been investigated by the HCFSM [7]. Milašinović [5] developed a harmonic coupled FSM by including the Green-Lagrange strain terms in the formulations, thus bending and membrane are coupled in the geometric stiffness to give more accurate buckling behaviour prediction, especially in large-deflection and elastic The interaction between the modes, shown in Fig. 8, is analysed taking into account the governing dynamic RDA modulus along the longitudinal axis. The HCFSM algorithm for geometric non-linear analysis has been described in [9] and [10]. The semi-analytical FSM eigenvalue analysis of a structure is simpler as it is explained in Section 4.1. H-section column is divided into 14 finite strips with 15 nodal lines. The main role has eigenvalues and corresponding eigenvectors. Fig. 9 presents first four local eigenvectors on length of 200 mm, respectively, while Fig. 10 presents first four global eigenvectors on length of 2310 mm. Though elastic buckling information for thin-walled members is not a direct predictor of capacity or collapse behaviour on its own, both the mode and the load are important proxies for the actual behaviour. In current design codes, such as AISI S100, New Zealand/Australia, and European Union, the design formulae are calibrated through the calculation of elastic critical buckling loads to predict the ultimate strength, thus the ability to calculate the associated elastic buckling loads is of great importance. Moreover, the buckling mode shapes are commonly employed into non-linear collapse modelling as initial geometric imperfections.  Viscoelastic buckling stress lags behind the elastic buckling stress across all modes, which is a consequence of the viscoelastic behaviour of materials. The viscoelastic behaviour is characterized by the delay time T D . As the length of the column is larger, the observed lag increases. Due to lag, the same column length does not always correspond to the same mode at the elastic and viscoelastic critical stress. However, at the column length of 2310 mm, the computed elastic critical stress of the first mode corresponds to test results [23]. Also, the computed viscoelastic stress at 2310 mm of the first mode has been previously confirmed in [7]. Fig. 12 shows the natural frequencies computed, as explained in Section 2.2, for the discussed column lengths. The frequency curves have jumps at the end of examined column length intervals, in contrast to the buckling curves. Such jumps within frequency curves are noticeable in both the elastic and viscoelastic solutions. Also, there is a noticeable lag between the viscoelastic and the elastic natural frequencies, analogue to the described critical stress behaviour.      Fig. 15 compares the values of viscoelastic natural frequencies as solved by the eigenvalue problem and the frequencies as computed by applying the critical stresses in Eq. (22). Fig. 16 presents the relative errors between the solutions as obtained by the eigenvalue problem and the solutions as computed by Eq. (22). The relative error of critical stresses is higher that the relative error of natural frequencies, because the buckling curve is curve of a higher order than the natural frequencies curve. The maximum relative error of only 1.6%/0.8% (at the approx. 957 mm) confirms the physical duality approach and effectively shows that it is sufficient to solve the eigenvalue problem only for the natural frequencies, while the critical stresses may be accurately computed through the equation of physical duality. It is well recognized that when using exact theory or semi-analytical FSM, as opposed to the more usual approximate finite element method, there exists a duality between buckling and natural vibration of an elastic structure, even though the former is a static problem whereas the latter is a dynamic one. The solution procedures for both problems are analogous because both can be essentially represented by a transcendental eigenvalue problem, the eigenvalues being natural frequencies in vibration problems and load factors in buckling problems. Thus the elastic buckling of a structure is often regarded as the degenerated case of its natural vibration at zero frequency. In this paper duality is proved for the both elastic and viscoelastic (or damage) structure using the RDA.     Fig. 18 contains a delay time curve on the lengths of columns from 200 to 900 mm. It is obvious that delay time curve is a sawtooth wave function, divided into segments between modes. Delay time segment peaks shorten when column lengths transition from lower to higher modes. However, the delay time is suddenly increased once the column transitions from the last local mode to the first global mode, because the total column length is engaged in the first global mode. This also causes a sharp increase both in the damage variable and the inelastic effective stress. Fig. 19 shows the damage variables (Eq. 38) on all column lengths from 200 to 4000 mm. The damage variables fall off for lengths where the column transitions from lower to higher modes. The damage variables fall off is greatest when transition occurs from the first to the second local buckling mode. The other falls become smaller with each successive transition to the next higher local mode. At the column length of the 957 mm the column transition from the highest reached local buckling mode into the first global buckling mode occurs. Simultaneously, the damage variable jumps to the maximum value of 1  , which means that the column does not survive this transition and fails instead. This behaviour correlates with the delay, as previously described.   It should be noticed, Figs. 19 and 20, that at lengths larger than 957 mm there are no more abrupt changes in any of the described physical properties, although the delay time constantly increases, as shown in Fig. 17. This is due to the absence of mode interaction between elastic and viscoelastic columns at lengths higher than 957 mm. After this length, the columns stay into the same (first) global buckling mode. In conclusion, this proves that the leaping changes of discussed physical properties are solely a consequence of the columns transitioning from one mode to another, whether it's local or global modes. Figs. 21 (a)-(b) show different zoomed-in views of the damage variables on column lengths from 200 to 900 mm, to better display the shape of the function. Fig. 21 (a) shows the function dropping off around mode transitions, whereas Fig. 21 (b) shows the Results of buckling effective stresses for the four observed column lengths are summarized in Table 4. Interaction between two local or local and global buckling modes occurs in intermediate length columns with near coincident critical stresses. In this case the columns do not have all elastic and inelastic characteristics in the same mode. If column is in the mode interaction of two local modes (a = 278 mm), the effective stress is smaller than elastic stress. If column is in the local and global mode interaction (a = 925 mm), the effective stress is much larger than elastic stress. Hence, the interaction between this two buckling modes (local and global) induced instantaneous catastrophic failure. If there is no mode interaction, columns have all elastic and inelastic characteristics in the same mode, local (a = 242 mm) or global (a = 2310 mm), and effective stresses converge to elastic stresses.

CONCLUSIONS
Numerical examples presented in this paper show the influence of mode interactions on effective stresses if the dynamic properties in undamaged state as well as in damaged state are used together in the inelastic buckling analysis of thin-walled columns. Numerical examples were computed by application of extensive hybrid parallelization. The mode interactions have great influence on damage variables and effective stresses around mode transitions. In its usual semi-analytical form, the FSM cannot be used to solve the mode interaction problem explained in this paper, because this technique ignores the important influence of interaction of the buckling modes when applied only for undamaged state of structure.