Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 54, 2, pp. 449-462, Warsaw 2016 DOI: 10.15632/jtam-pl.54.2.449 SOME ASPECTS OF DYNAMIC COUPLED RESPONSE OF FUNCTIONALLY GRADED THIN-WALLED COLUMNS WITH SQUARE CROSS-SECTIONS UNDER IN-PLANE PULSE COMPRESSION Zbigniew Kołakowski Lodz University of Technology, Department of Strength of Materials, Łódź, Poland; e-mail: zbigniew.kolakowski@p.lodz.pl Andrzej Teter Lublin University of Technology, Department of Applied Mechanics, Lublin, Poland; e-mail: a.teter@pollub.pl The present paper deals with a dynamic coupled response of functionally graded columns with a quadratic cross-section subjected to an in-plane pulse loading. An Al-TiC metal- -ceramic material is applied. It is assumed that functionally graded materials (FGMs) are subject to Hooke’s law. The thin-walled structures are simply supported at the ends. This study is devoted to the stability problem of rectangular dynamic pulse load. The effects of temperature, wave propagation and damping are neglected. In order to obtain the equations of motion of individual plates, the classic laminate plate theory (CLPT) has been modified in such a way that it additionally accounts for all components of inertial forces. A plate model is adopted for the structures. The problem of an interaction of the global mode with the local ones is concerned (i.e., a three-modes approach). Attention has been focused on some unexpected aspects related to dynamic interactive buckling of columns having two axes of the cross-section symmetry. In the present study, a new approach to the description of this phenomenon, based onKoiter’s theory, has been applied. Keywords: FGM, dynamic response, interactive buckling, thin-walled structures, compres- sion, pulse load 1. Introduction The dynamic buckling or dynamic response takes place when a pulse load of a mean amplitude and a pulse duration comparable to the fundamental natural flexural vibration period occurs in compression of the thin-walled column. In this case, effects of damping can be neglected in practice.When the amplitude of load is high, then the structure can vibrate very strongly or can move divergently, which is caused by dynamic buckling. One can determine the critical ampli- tude of load using various criteria. In the literature on this problem, a lot of criteria concerning dynamic stability have been adopted. Themost widely used is the Budiansky-Hutchinson crite- rion (Kubiak, 2007, 2013), in which it is assumed that the dynamic stability loss occurs when the maximum structure deflection grows rapidly at a small variation in the load amplitude. Other criteria have been discussed in many papers: Ari-Gur and Simonetta (1979), Petry and Fahlbusch (2000), Kubiak (2007, 2013), Teter (2011). Dynamic global and local buckling instabilities of component functionally graded plates (the so-called FG plates) of structures subjected to conservative loads have been taken into account. The problem of an interaction of the global mode with the local ones is very interesting. The concept of interactive buckling involves the general asymptotic theory of stability. Among all versions of the general nonlinear theory, Koiter’s theory (van der Heijden, 2009) of conserva- tive systems is the most popular one (Kołakowski et al., 1999; Teter and Kołakowski, 2004; Kołakowski andKubiak, 2005; Kołakowski and Królak, 2006). 450 Z. Kołakowski, A. Teter In the present study, the classical laminate plate theory (CLPT) (Jones, 1999; Reddy, 2004) is employed to obtain the governing equations of the thin FG plate equilibrium. In order to obtain the equations of individual plates for the asymptotic analytical-numerical method, the nonlinear theoryof compositeplateshasbeenmodified in suchaway that it additionallyaccounts for forces of inertia. The differential equations of motion have been obtained from Hamilton’s Principle, taking into account Lagrange’s description, full Green’s strain tensor, the second Piola-Kirchhoff’s stress tensor and all components of inertia forces. The study is based on the numerical method of the transition matrix using Godunov’s orthogonalization (Kołakowski et al., 1999; Teter and Kołakowski, 2004; Kołakowski and Kubiak, 2005; Kołakowski and Królak, 2006). A plate model of the column has been adopted in the study to describe global buckling, which leads to lowering the theoretical value of the limit load. The solution method assumed in this study enables analyzia of interactions of all buckling modes. The nonlinear equations of dynamic instability are solved with the modified Runge-Kutta method. The nonlinear analysis of Functionally Gradient plates and shells devoted to basic types of loads is covered in the monograph by Hui-Shen (2009). The shear deformation effect is em- ployed in the framework of Reddy’s higher order shear deformation theory (HSDT) (Reddy, 2000; Reddy, 2004). Reddy (2000) presents a comparison of applications of the first order shear deformation theory (FSDT) and the classical lamination plate theory (CLPT) to functionally graded plates. The discrepancy between both theories is of 2% in the calculated deflections of the plates under analysis. The buckling and postbuckling problem of FG plates is discussed, for example, in papers by Reddy (2000), Samsam Shariata et al. (2005), Kołakowski et al. (2015). Due to the complexity of buckling problems of FG structures under compoundmechanical and thermal loads, the finite element method (FEM) is the only solution possible in many cases. Therefore, in the literature, one can find many papers which present results of a solution to different problems of FG structure buckling obtained with an application of the FEM, see for example: Birman and Byrd (2007), Panda and Ray (2008), Na and Kim (2009), Kołakowski et al. (2015). 2. Formulation of the problem Long thin-walled prismatic columns of length l, composed of plane rectangular plate segments interconnected along longitudinal edges and simply supported at both ends, are considered. All materials the FG plates are made of are subject to Hooke’s law. The material properties are assumed to be temperature independent. A plate model is adopted for the structures. Wave propagation and damping effects have been neglected in the present study, as it is done in the majority of works devoted to dynamic stability. In the present study, the classical laminate plate theory (CLPT) (Jones, 1999; Reddy, 2004) is employed to obtain the governing equations of the thin FG structure equilibrium (Kołakowski and Królak, 2006; Panda and Ray, 2008). For the plate component, precise geometrical relationships are assumed in order to enable consideration of both out-of-plane and in-plane bending of the plate (Kubiak, 2007, 2013; Teter, 2007, 2011; Kołakowski et al., 1999; Teter andKołakowski, 2001, 2004, 2005; Kołakowski andKubiak, 2005; Kołakowski andKrólak, 2006) ε=    εx εy 2εxy    =    u,x+ 1 2 (w2,x+v 2 ,x+u 2 ,x) v,y + 1 2 (w2,y+u 2 ,y +v 2 ,y) u,y+v,x+w,xw,y+u,xu,y +v,xv,y    κ=    κx κy κxy    =    −w,xx −w,yy −2w,xy    (2.1) whereu, v,w are components of the displacement vector of the plate in thex, y, z axis direction, respectively, and the plane xy overlaps the mid-plane before its buckling. Some aspects of dynamic coupled response of functionally graded... 451 According to the rule of mixture, the properties of the functionally graded material (E – Young’s modulus, ν – Poisson’s ratio, ρ – density) can be expressed as follows E(z) =Em+(Ec−Em) (z h + 1 2 )q ν(z)= νm+(νc−νm) (z h + 1 2 )q ρ(z) = ρm+(ρc−ρm) (z h + 1 2 )q (2.2) where the indices m and c refer to the metal and ceramic material, respectively, and q is the fraction exponent. Using the classical laminate plate theory (CLPT) (Jones, 1999; Reddy, 2004), the stress and moment resultants (N,M) are defined as (Jones, 1999; Teter andKołakowski, 2005, Kołakowski and Królak, 2006) { N M } = [ A B B D ]{ ε κ } (2.3) whereA,B,D are extensional, coupling and bending stiffnessmatrices, respectively, for the FG structure. Their components are listed below A11 =A22 = h/2∫ −h/2 E(z) 1−ν2(z) dz A12 =A21 = h/2∫ −h/2 E(z)ν(z) 1−ν2(z) dz A66 = h/2∫ −h/2 E(z) 2[1+ν(z)] dz A16 =A61 =0 B11 =B22 = h/2∫ −h/2 E(z) 1−ν2(z) z dz B12 =B21 = h/2∫ −h/2 E(z)ν(z) 1−ν2(z) z dz B66 = h/2∫ −h/2 E(z) 2[1+ν(z)] z dz B16 =B61 =0 D11 =D22 = h/2∫ −h/2 E(z) 1−ν2(z) z2 dz D12 =D21 = h/2∫ −h/2 E(z)ν(z) 1−ν2(z) z2 dz D66 = h/2∫ −h/2 E(z) 2[1+ν(z)] z2 dz D16 =D61 =0 (2.4) Due to the presence of the nontrivial submatrix B, the coupling between extensional and bending deformations exists as it is in the case of unsymmetrical laminated plates (Jones, 1999; Kołakowski et al., 2015). An extensional force results not only in extensional deformations, but also bendingof theFGplate.Moreover, such aplate cannot be subjected to themomentwithout suffering simultaneously from extension of the middle surface. The nonlinear problem of dynamic stability has been solved with the asymptotic perturba- tion method. Let λ be a load factor. The displacement fieldsU and the sectional force fieldsN (Koiter’s type expansion for the static buckling problem (van der Heijden, 2009)) have been expanded into power series with respect to the dimensionless amplitude of the r-thmode deflec- tion ζr (normalized in the given case by the condition of equality of the maximum deflection to the thickness of the first component plate h1) (Kołakowski et al., 1999; Teter and Kołakowski, 2001, 2004, 2005; Teter, 2007, 2011; Kołakowski and Kubiak, 2005; Kołakowski and Królak, 2006; Kubiak, 2007, 2013; van der Heijden, 2009) U=λ(t)U0+ ζr(t)Ur +ζr(t)ζq(t)Uqr + . . . N=λ(t)N0+ζr(t)Nr+ ζr(t)ζq(t)Nqr+ . . . (2.5) where the pre-buckling static fields are U0, N0, the first nonlinear order fields are Ur, Nr (the eigenvalues and eigenvectors problems) and the second nonlinear order fields – Uqr, Nqr, 452 Z. Kołakowski, A. Teter respectively. The range of indexes is [1,J], where J is the number of interacting modes. The summation is supposed on the repeated indexes. If the structure contains the geometric imperfections U ∗ (only linear initial imperfections determined by the shape of r-th buckling modes i.e., U ∗ = ζ∗rUr), then the total energy of the structures can be written in the form (Schokker et al., 1996; Kołakowski et al., 1999; Teter and Kołakowski, 2004; Teter, 2007, 2011; Kołakowski and Kubiak, 2005; Kołakowski and Królak, 2006; Kubiak, 2007, 2013) Π =− 1 2 σ2(t)a0+ 1 2 J∑ r=1 arζ 2 r(t) ( 1− σ(t) σr ) + 1 3 J∑ p J∑ q J∑ r apqrζp(t)ζq(t)ζr(t) + 1 4 J∑ r brrrrζ 4 r(t)− J∑ r σ(t) σr arζ ∗ rζr(t)+ 1 2 J∑ r mrζ 2 r,t(t) (2.6) Then, Hamilton’s principle leads to the following Lagrange’s equation (i.e. equations of mo- tion) ζr,tt+Ω 2 rζr+ω 2 r ( apqrζpζq+ brrrrζ 3 r − σ σr ζ∗r + . . . ) =0 for r=1, . . . ,J (2.7) where ζr is the dimensionless amplitude of the r-th buckling mode, σr, ωr, ζ ∗ r – critical stress instead of the load parameter λr of the r-th bucklingmode, circular frequency of free vibrations and dimensionless amplitude of the initial deflection corresponding to the r-th buckling mode, respectively. In equations of motion (2.7), the inertia forces of the pre-buckling state and the second order state have been neglected (Sridharan et al., 1984; Schokker, 1996; Warmiński and Teter, 2012). The coefficients in equilibrium equations (2.7) are given in papers: Sridharan et al. (1984), Kołakowski (1996), Schokker (1996), Kołakowski et al.) (1999), Teter and Kołakowski (2004), Teter (2007, 2011), Kołakowski and Kubiak (2005), Kołakowski and Królak (2006), Kubiak (2007, 2013),Warmiński andTeter (2012), Kołakowski andMania (2013). In the former parts of this paper, in relationships (2.7), λr has been replaced with σr, whereas λ with σ, correspondingly. The initial conditions have been assumed in the form ζr(t=0)=0 ζr,t(t=0)=0 (2.8) In equation (2.7), the quantity Ω2r depending on the values of σ/σr can take the following values: • if σ/σr < 1, then Ω 2 r =ω 2 r(1−σ/σr)> 0, whereΩr can be called the equivalent angular velocity. The linear general solution to equation (2.9) in the case of an ideal structure (that is to say, for ζ∗r =0) are trigonometric functions; • if σ/σr =1, thenΩ 2 r =0; • if σ/σr < 1, then Ω 2 r < 0, where Ωr can be called the equivalent growing function and then the linear general solution to equation (2.7) for ζ∗r =0 are hyperbolic functions. The nonlinear static stability (i.e., for ζr,tt in (2.7)) of thin-walled multilayer structures in the first order approximation of Koiter’s theory is solved with themodified analytical-numerical method (ANM) presented by Kołakowski and Królak (2006). The analytical-numerical method (ANM) should consider also the second order approximation of the theory in the analysis of postbuckling of elastic composite structures. The second order postbuckling coefficients were estimated with the semi-analytical method (SAM) (Kołakowski, 1996) modified by the solution method in Kołakowski and Mania (2013). In the semi-analytical method (SAM) for static pro- blems, one postulates to determine approximated values of the coefficients brrrr in (2.7) on the Some aspects of dynamic coupled response of functionally graded... 453 basis of the linear buckling problem. This approach allows the values of the coefficients apqr in (2.7) – according to the applied nonlinear Byskov and Hutchinson theory (Byskov and Hut- chinson, 1977) – to be determined precisely. The natural frequencies have been determined analogously as in Teter and Kołakowski (2003), whereas the problem of interactive dynamic buckling (2.7) has been solved by means of the Runge-Kutta numerical method modified by Hairer andWanner. In thepresentpaper,weassumethat thebucklingmodesare the sameas thevibrationmodes, so the solution to the eigenvalue problem is sought for various values of the m-th harmonic. Values of the natural frequencies are determined taking into account all components of inertia forces. For static problems, Koiter and van der Neut (Kubiak, 2013) have proposed a technique in which an interaction of the overall modewith two local modes having the samewavelength (i.e., a three-mode approach J =3) has been considered. The fundamental local mode is henceforth called “primary” and the nontrivial highermode (having the samewavelength as the “primary” one), corresponding to themode triggered by the overall long-wavemode, is called “secondary”. In total energy, the coefficients of the cubic terms ζ1ζ 2 2, ζ1ζ 2 3 and ζ1ζ2ζ3 (where ζj is the amplitude of the r-th bucklingmodeand the index is: 1 for the globalmode, 2 for theprimary local buckling mode, and 3 for the secondary local mode) are the key terms governing the interaction. In the analysis of the column with doubly symmetric cross-sections, the coefficients of ζ1ζ 2 2 and ζ1ζ 2 3 terms – the coefficients apqr of non-linear system of equations (2.7) – vanish. In the paper byKołakowski et al. (2015), nonlinearKoiter’s theory has been used to explain the effect of the imperfection sign (sense) on local postbuckling equilibrium paths of plates made of functionally gradedmaterials (FGMs). In the case of the FG plate, nonzero first-order sectional inner forces that cause an occurrence of nonzero postbuckling coefficients responsible for sensitivity of the system to imperfections appear. It results in the fact that postbuckling equilibrium paths of plate structures made of FGMs are unsymmetrically stable. This explains differences in the plate responsedependenceon the imperfection sign (sense).On the other hand, in Kołakowski and Mania (2015), an analysis of the influence of the imperfection sign on the dynamic postbuckling equilibrium paths of the FG square plates has been continued. In the numerical calculations of dynamic interactive buckling, a rectangular shape of the in-plane pulse loading (i.e., σ(t) = σD for 0 ¬ t ¬ T1 and σ(t) = 0 for T < t1) equal to the fundamental period of the natural frequency T1 =2π/ω1 is considered. The prebuckling solution to the FG plate consisting of homogenous fields is assumed as in Kołakowski andKrólak (2006) u(0) = ( l 2 −x ) ∆ v(0) = y∆ A12 A22 w(0) =0 (2.9) where∆ is the actual loading. This loading of the zero state is specified as a product of the unit loading and the scalar load factor ∆. Taking into account relationship (2.3), the inner sectional forces of the prebuckling (i.e., unbending) state for the assumed homogeneous field of displacements (2.9) are expressed by the following relationships before the redistribution of forces in the plate due to plate deformations (Kołakowski et al., 2015; Kołakowski andMania, 2015; Kołakowski, 2016) N (0) x =− ( A11− A212 A22 ) ∆ M (0) x =− ( B11−B12 A12 A22 ) ∆ N (0) y =0 M (0) y =− ( B12−B22 A12 A22 ) ∆ N (0) xy =0 M (0) xy =0 (2.10) 454 Z. Kołakowski, A. Teter The assumed field of displacements and the field of inner forces, corresponding to it for the prebuckling state, fulfil the equilibrium equations for the zero state as an identity. In Kołakowski and Królak (2006), an unbending, prebuckling state, i.e., a distribution field of the zero state according to (2.9), has been assumed. Dependence (2.3) for the zero state (i.e., prebuckling) takes the form N0 = { N(0) M(0) } = [ A B B D ]{ ε (0) 0 } (2.11) It results in an occurrence of the nonzero inner sectional forces (2.10) N (0) x , M (0) x , M (0) y in the FGplate for the zero state. Special attention shouldbe paid to the fact that nonzeromagnitudes of the sectional momentsM (0) x andM (0) y appear due to the effects of deformations of themiddle surface (i.e., membrane deformations) resulting from the nontrivial coupling submatrix B and deformations of the middle surface, and not due to an appearance of curvatures of the middle surface.Thesemomentsaffect obviously thevalues of critical loads and thevalues ofpostbuckling coefficients (Kołakowski et al., 2015; Kołakowski andMania, 2015). 3. Analysis of the results Prismatic thin-walled beam-columns with quadratic cross-sections subjected to an axial pulse compression have been considered. The columns aremade of the same FGM subject to Hooke’s law. A schematic view of the column and its dimensions are shown in Fig. 1. Fig. 1. Thin-walled cross-section of the beam-column Identically as in Kołakowski et al. (2015), Kołakowski and Mania (2015) two different ma- terials the columns are made of have been considered, namely: • isotropic material – steel – with the following material properties: Young’s modulus E =210GPa, Poisson’s ratio ν =0.3, and density ρ=7850kg/m3; • Al-TiC functionally gradedmaterial (FGM).The componentmaterial properties are given: Al Young’s modulus: Em = 69GPa, TiC Young’s modulus: Ec = 480GPa, Al Poisson’s ratio: νm =0.33,TiCPoisson’s ratio: νc =0.20,Al density:ρm =2700kg/m 3, TiCdensity: ρc =4920kg/m 3. The fraction exponent q (2.2) is equal to 1.0. Three variants of the beam-column structure have been considered, namely: • variant I – isotropic column (the so-called reference variant); • variant II – QCMFG column (ceramics inside the cross-section, metal outside); • variant III – QMCFG column (metal inside the cross-section, ceramics – outside). The isotropic beam-column (variant I) has symmetrical stable postbuckling local equilibrium paths,whereas theFGcolumn (variants II and III) hasnonsymmetrical stable postbuckling local Some aspects of dynamic coupled response of functionally graded... 455 equilibrium paths (Kołakowski et al., 2015; Kołakowski and Mania, 2015). Due to the above aspects, two variants of the FG beam-column –QCM andQMC – have been considered. For each of these three structural variants, values of critical stresses σr: the global (Euler) buckling mode σ1, the lowest primary local buckling mode σ2 and the secondary local buckling mode σ3 have been determined, respectively. Moreover, numbers of halfwaves m corresponding to the critical loads have been given. The global mode occurs at m = 1 and the local modes at m > 1. Also, values of frequencies of free vibrations ωr and a period of free vibrations Tr corresponding to the buckling modes under analysis have been defined as well. An interactive dynamic response of beam-columns to the load whose duration corresponds to the fundamental period of flexural free vibrations of unloaded columns T1 =2π/ω1 has been investigated. The assumed time of the rectangular load pulse duration corresponds to the beam-column quasi-static load for local modes (i.e., T2 > 2T1, T3 > 2T1). The tracing time of the struc- ture dynamic response assumed as t∗ = 1.5T1 has been analysed. The dynamic load factor DLF is defined as a ratio of the dynamic load to the minimal critical value of the static load, DLF =σD/min(σ1,σ2,σ3)=σD/σ2. In Budiansky and Hutchinson (1966), it has been shown for the two-mode approach that when the frequencies of free vibrations differ at least twice, then in this case the dynamic term corresponding to a higher frequency of free vibrations can be neglected in the equations of motion.On the other hand, inKołakowski (2016), attention has beenpaid to unexpected aspects of interactive dynamic buckling in the case of beam-columns having one axis of symmetry of the cross-section. Here, in order to show qualitatively different dynamic responses, a two- or three-mode approach could be applied. It has allowed the uncoupling of the equations ofmotion for the case of the inner combined resonance, when (ω1+ω2)/ω3 ≈ 1 (Nayfeh andMook, 1979). It is not possible touse this approach in thepresent studybecause the interaction occurs only via the coefficient of the cubic term ζ1ζ2ζ3 in total energy (2.6) governing themode interaction. Therefore, an approach described in Budiansky and Hutchinson (1966) has been applied. The secondary local bucklingmode in the case of static issues is a supplementarymode that enables one to account for the effect of the primary local mode on the global flexural buckling mode. In the further part of the study, it has been assumed that when “complete” three dynamic equations of motion (2.7) are considered, then the case is referred to as case I. When the dynamic term is neglected for the primary local buckling mode (i.e., for r = 2), that is to say, when ζ2,tt = 0 is assumed in equations (2.7) – it is case II, whereas when the term ζ3,tt = 0 is neglected for the secondary local mode (i.e., for r=3) – case III. When the dynamic term is neglected in (2.7), which corresponds to case II, one of the equ- ations is a static equation of the third order with respect to the amplitude of deflection ζ2, which is solved on the basis of analytical formulae. Due to this, wrong conditions in the equ- ations could be avoided and the numerical solution could be stabilised. It means that in the formula for total energy of the system (2.6), the expression for kinetic energy 0.5m2ζ 2 2,t corre- sponding to the primary local mode (case II) could be neglected. Analogously, in case III – the expression for the energy 0.5m3ζ 2 3,t corresponding to the secondary local mode (i.e., ζ3) could be omitted. Further on in the study, the following index notations of dimensionless amplitudes of deflec- tion have been introduced for the cases under consideration, namely: case I – ζr for r=1,2,3; case II – ζ̂r for r = 1,2,3; case III – ζ̃r for r = 1,2,3. In the calculations of dynamic stability, the following values of imperfections have been assumed: ζ∗1 =1.0, ζ ∗ 2 =0.2, ζ ∗ 3 =0.1. Contrary to other works by the authors (e.g., Kołakowski, 1996, 2016; Kołakowski et al., 1999; Teter and Kołakowski, 2001, 2004, 2005; Teter, 2001, 2007, 2011; Kołakowski and Kubiak, 2005; Koła- kowski andKrólak, 2006; Kubiak, 2007; Kołakowski andMania, 2013; Kubiak, 2013), the most unfavourable combination of initial imperfection signs ζ∗r (for r=1,2,3) has not been assumed 456 Z. Kołakowski, A. Teter tomake the effects of nonsymmetrical stable postbuckling equilibrium path and the influence of the connection between the adjacent plates of the cross-section more visible. Values of the critical dynamic load factorsDLFcr have beendetermined fromtheBudiansky- -Hutchinson criterion (Budiansky amd Hutchinson, 1977; Kubiak 2007, 2013) in which it is assumed that the loss of dynamic stability occurswhen thevelocitywithwhich thedisplacements grow is the highest for a certain force amplitude. The values of DLFcr presented correspond with some accuracy to the maximum values of deflections max(ζ1) within the applicability of the assumed theory (i.e., the totalmaximumdeflection of the column is at least a hundred times as high as the cross-section wall thickness), and not to asymptotic values (Kubiak, 2013). 3.1. Variant I – isotropic columns Variant I has been assumed as the reference one because in the case of an isotropic column, the plates constituting the cross-section are characterized by the symmetry of cross-sectional uniformity (i.e., coupling stiffness matrixB=0 in (2.3)). Table 1. Solutions to the eigenproblem for an isotropic column – variant I (reference) r σr [MPa] m [–] ωr [rad/s] Tr [ms] 1 1185.6 1 761 8.256 2 72.28 16 3013 2.085 3 103.66 16 3609 1.741 First, the eigenvalues of theproblemgiven inTable 1havebeendetermined, namely: values of critical stresses corresponding to the global Euler buckling σ1, the lowest primary local buckling mode σ2 and the secondary local buckling mode σ3, the frequencies of free vibrations and the periods corresponding to them, as well as a number of halfwavesm of the eigenmodes along the longitudinal direction. The ratios of frequencies of free vibrations are equal to: ω2/ω1 = 3.96, ω3/ω1 =4.74, ω3/ω2 =1.197, (ω1+ω2)/ω3 =1.045, respectively. In Fig. 2a, absolute values of the maximum amplitudes of global defections (for r=1) (i.e., ζ1, ζ̂1, ζ̃1) as a function of DLF for cases I-III under consideration, whereas in Fig. 2b – for local deflections of the modes ζ2, ζ̂3 and ζ̃2 for cases I, II and III are shown, correspondingly. Fig. 2. Maximum dimensionless global deflections as a function ofDLF for an isotropic column (variant I): (a) global modes, (b) local modes Quantative differences in the component amplitudes of global deflections (i.e., ζ1, ζ̂1 and ζ̃1) and the local ones (i.e., ζ2, ζ̂3 and ζ̃2) can be observed.As can be easily seen in these two figures, the maximum values of the amplitudes of the global and local deflection for each case under Some aspects of dynamic coupled response of functionally graded... 457 investigation occur for the same values of DLF (i.e., DLF ≈ 1.69 for case I, DLF ≈ 1.44 for case II, andDLF ≈ 1.56 for case III, respectively). When DLF = 1.434, then the amplitude of dynamic load is σD = σ3 and, moreover, σ3/σ2 = 1.434 (Table 1). When t ¬ T1, then for 1.0 ¬ DLF < 1.434 and Ω 2 1 > 0, Ω 2 2 < 0, Ω23 > 0, whereas forDLF > 1.434, Ω 2 1 > 0,Ω 2 2 < 0,Ω 2 3 < 0, respectively. Themaximumvaluemax(ζ1) ismore than twice ashighasmax(ζ̂1) andmax(ζ̃1),whereas the differences betweenmax(ζ̂1) andmax(ζ̃1) are equal to approx. 20%.The differences between the maximumvalues of local amplitudes (i.e., max(ζ2), max(ζ̃2) andmax(ζ̂3)) are not so significant. The value max(ζ2) is more than 30% higher than the max(ζ̃2) and max(ζ̂3). The maximum values of global deflections are by one order of magnitude higher than for local deflections. As can be seen in Fig. 2 for case II, the deflections ζ̂1 and ζ̂3 forDLF > 1.7 aremany times lower than the deflections in cases I and III. For case II, the lowest value of DLFcr and the lowest value max(ζ̂1)≈ 200 for DLF ≈ 1.44 have been obtained (Fig. 2a), which correspond to: σD ≈σ3, σ3/σ2 =1.434 andΩ 2 3 ≈ 0. When t¬T1 and forDLF =1.56 in case III, we have max(ζ̃1)≈ 250,Ω1 =723,Ω2 =2259 and Ω3 = 1081, and for DLF = 1.69 in case I, we have max(ζ1)< 500, Ω1 = 720, Ω2 = 2498 andΩ3 =1518. When T1 0, Ω 2 2 < 0 and Ω 2 3 > 0, whereas for DLF > 1.4307, we have, Ω 2 1 > 0, Ω 2 2 < 0 and Ω23 < 0, correspondingly. The maximum values max(ζ1) and max(ζ̃1) are twice as high as max(ζ̂1), and max(ζ1)/max(ζ̃1)≈ 1.5. The differences between the maximum values of local deflections, i.e., ζ2, ζ̂3 and ζ̃2, are not so high and they do not exceed 1.5 times. The values max(ζ1), max(ζ̂1), max(ζ̃1) for globalmodes are higher by one order ofmagnitude thanmax(ζ2),max(ζ3),max(ζ̂3), max(ζ̃2) for local modes. In Fig. 3, for the duration of load (i.e., for t ¬ T1), we have Ω1 = 1224, Ω2 = 3235 and Ω3 =1251 for DLF =1.51, whereas for DLF =1.81 – Ω1 =1213, Ω2 =4125 and Ω3 =2849, respectively. For T1 500 andmax(ζ2)> 25. ForDLF =1.4383, we have σ3/σ2 = 1.4383. When t ¬ T1, then for 1.0 < DLF < 1.4383, we have Ω 2 1 > 0, Ω22 < 0 and Ω 2 3 > 0, whereas at DLF > 1.4383, Ω 2 1 > 0, Ω 2 2 < 0 and Ω 2 3 < 0, respectively. The lowest values of maximum deflections are for case II, that is to say, for max(ζ̃1) and max(ζ̂3), analogously as for variant II. Fig. 4. Maximum dimensionless global deflections as a function ofDLF for variant III: (a) global modes, (b) local modes In Fig. 4, under the dynamic loading (for t ¬ T1) for DLF = 1.45, we have Ω1 = 1231, Ω2 = 3015 and Ω3 ≈ 0, whereas for DLF = 1.85 – Ω1 = 1218, Ω2 = 4157and Ω3 = 28488, respectively. At T1 < t ¬ t ∗, we have (ω1 +ω2)/ω3 = 1.069. The critical values of DLFcr are equal to: case I – DLFcr = 1.56, case II – DLFcr = 1.69, case III – DLFcr = 1.40. When t ¬ T1, then for the assumed values of DLF, we have Ω 2 1 > 0, Ω 2 2 < 0 and Ω 2 3 > 0. For the three values ofDLF under consideration and, simultaneously, for the three cases, themaximum global deflection during the response tracing time (i.e., 0¬ t¬ t∗) occurs for case II –max(ζ̂1). When the time functions for variants II and III are compared, some insignificant differences can be observed. The maximum values of global deflections (i.e., max(ζ1), max(ζ̂1), max(ζ̃1)) take place for T1 1.45, we have Ω21 > 0, Ω 2 2 < 0 and Ω23 < 0. It means that the amplitudes of local deflections grow exponentially and become dominating in a short time. A dynamic interaction of the globalmode and the local ones results in a very dramatic increase in the amplitudes of deflections. 460 Z. Kołakowski, A. Teter For the time functions under consideration for all variants and cases and the assumed valu- es of load coefficients DLF, the amplitudes of global deflections attain their maximum values after the pulse load finishes (e.g., for T1 < t ¬ t ∗). For this time range, free unsteady vibra- tions occur in the transient process. To the authors’ knowledge, the only phenomenon that can explain this effect is the inner combined resonance which takes place accordingly to the theory of vibration (Nayfeh and Mook, 1979) for steady processes, correspondingly, in the cases when (ω1+ω2)/ω3 ≈ 1.0. The three cases under analysis complywith the above-mentioned condition. A change in the tracing time of the dynamic response t∗ can exert an influence on themaximum values of amplitudes of deflections. In the study, a long enough tracing time has been assumed, as it is equal to t∗ = 1.5T1. However, shortening of this time (e.g., up to t ∗ = 1.3T1) will not exert any influence on the general conclusions drawn here. 4. Conclusions Dynamic interactive buckling of thin-walled FGMcolumnswith a square cross-section is discus- sed. Three cases are considered. Case I corresponds to analysis of the dynamic response of the FG structure for “complete” equations of motion (2.7), case II refers to the situation when the dynamic effect is neglected for the primary local bucklingmode (i.e., r=2), whereas case III – for the secondary local mode (i.e., r=3). Such an approach has been pointed out byBudiansky and Hutchinson (1966). When the structure has two axes of symmetry of the cross-section, the interactive buckling occurs only via the coefficient of the cubic term ζ1ζ2ζ3 in total potential energy (2.6). The solutionmethod applied in this studyhas not been so efficient as in the case of the trapezoidal cross-section with one axis of symmetry. InKołakowski (2016), however, a quite different approach to the solution of the problemhas been assumed.Thus, the above conclusions should be subject to further, thorough investigations. Acknowledgements This publication is a result of the research works carried out within the project subsidized over the years 2012-2015 by the state funds designated for the National Science Centre (NCN UMO- -2011/01/B/ST8/07441). References 1. Ari-Gur J., Simonetta S.R., 1979, Dynamic pulse buckling of rectangular composite plates, Composites Part B, 28B, 301-308 2. Birman V., Byrd L.W., 2007, Modeling and analysis of functionally materials and structures, Applied Mechanics Review, 60, 195-216 3. Budiansky B., Hutchinson J.W., 1966, Dynamic buckling of imperfection-sensitive structures, Proceedings of the 11th International Congress of Applied Mechanics, Goetler H. (Ed.), Munich, 636-651 4. Byskov E., Hutchinson J.W., 1977, Mode interaction in axially stiffened cylindrical shells, AIAA Journal, 15, 7, 941-948 5. Hui-Shen S., 2009,Functionally GradedMaterials –Nonlinear Analysis of Plates and Shells, CRC Press, Taylor & Francis, London 6. Jones R.M., 1999,Mechanics of Composite Materials, 2nd ed. Taylor & Francis, London 7. Kołakowski Z., 1996, A semi-analytical method for the analysis of the interactive buckling of thin-walled elastic structures in the second order approximation, International Journal of Solids Structures, 33, 25, 3779-3790 Some aspects of dynamic coupled response of functionally graded... 461 8. Kołakowski Z., 2016, Some aspects of interactive dynamic stability of thin-walled trapezoidal FGM beam-columns under axial load,Thin-Walled Structures, 98, 431-442 9. Kołakowski Z., Królak M., 2006, Modal coupled instabilities of thin-walled composite plate and shell structures,Composite Structures, 76, 303-313 10. Kołakowski Z., KrólakM.,Kowal-MichalskaK., 1999,Modal interactive buckling of thin- -walled composite beam-columns regardingdistortional deformations, International Journal of En- gineering Science, 37, 1577-1596 11. Kołakowski Z., Kubiak T., 2005, Load-carrying capacity of thin-walled composite structures, Composite Structures, 67, 417-426 12. Kołakowski Z., Mania J.R., 2013, Semi-analytical method versus the FEM for analysis of the local post-buckling,Composite Structures, 97, 99-106 13. Kołakowski Z., Mania J.R., 2015, Dynamic response of thin FG plates with a static unsym- metrical stable postbuckling path,Thin-Walled Structures, 86, 10-17 14. Kołakowski Z., Mania R.J., Grudziecki J., 2015, Local nonsymmetric postbuckling equili- brium path in thin FGM plate, Eksploatacja i Niezawodność (Maintenance and Reliability), 17, 135-142 15. Kubiak T., 2007, Criteria of dynamic buckling estimation of thin-walled structures,Thin-Walled Structures, 45, 10/11, 888-892 16. Kubiak T., 2013, Static and Dynamic Buckling of Thin-Walled Plate Structures, Springer 17. Na K.S., Kim J.H., 2009, Volume fraction optimization of functionally graded composite panels for stress reduction and critical temperature,Finite Elements in Analysis and Design, 45, 845-851 18. Nayfeh A.H., Mook D.T., 1979,Nonlinear Oscilations, Wiley, NewYork 19. Panda S., Ray M.C., 2008, Nonlinear finite element analysis of functionally graded plates inte- grated with patches of piezoelectric fiber reinforced composite, Finite Elements in Analysis and Design, 44, 493-504 20. Petry D., Fahlbusch G., 2000, Dynamic buckling of thin isotropic plates subjected to in-plane impact,Thin-Walled Structures, 38, 267-283 21. Reddy J.N., 2000, Analysis of functionally graded plates, International Journal for Numerical Methods in Engineering, 47, 663-684 22. Reddy J.N., 2004,Mechanics of Laminated Composite Plates and Shells, CRCPress, London 23. Samsam Shariata B.A., Javaherib R.B.A., Eslami B.A., 2005, Buckling of imperfect func- tionally graded plates under in-plane compressive loading,Thin-Walled Structures, 43, 1020-1036 24. Schokker A., Sridharan S., Kasagi A., 1996, Dynamic buckling of composite shells,Compu- ters and Structures, 59, 1, 43-55 25. Sridharan S., Benito R., 1984, Columns: static and dynamic interactive buckling, Journal of Engineering Mechanics ASCE, 110, 1, 49-65 26. Teter A., 2007, Static and dynamic interactive buckling of isotropic thin-walled closed columns with variable thickness,Thin-Walled Structures, 45, 10/11, 936-940 27. Teter A., 2011, Dynamic critical load based on different stability criteria for coupled buckling of columns with stiffened open cross-sections,Thin-Walled Structures, 49, 589-595 28. TeterA.,KołakowskiZ., 2001,Lowerboundestimationof load-carryingcapacityof thin-walled structures with intermediate stiffeners,Thin-Walled Structures, 39, 8, 649-669 29. Teter A., Kołakowski Z., 2003, Natural frequencies of thin-walled structures with central intermediate stiffeners or/and variable thickness,Thin-Walled Structures, 41, 291-316 30. Teter A., Kołakowski Z., 2004, Interactive buckling and load carrying capacity of thin-walled beam-columns with intermediate stiffeners,Thin-Walled Structures, 42, 211-254 462 Z. Kołakowski, A. Teter 31. TeterA.,Kołakowski Z., 2005,Buckling of thin-walled composite structureswith intermediate stiffeners,Composite Structures, 69, 4, 421-428 32. Van der Heijden A.M.A. (ed.), 2009, WT Koiter’s Elastic Stability of Solids and Structures, Cambridge University Press 33. Warmiński J., Teter A., 2012, Non-linear parametric vibrations of a composite column under uniform compression, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 226, 8, 1921-1938 Manuscript received Jule 29, 2015; accepted for print September 2, 2015