Finite element models in vibration analysis of two - Dimensional functionally graded beams

MINISTRY OF EDUCATION AND TRAINING VIETNAM ACADEMY OF SCIENCE AND TECHNOLOGY GRADUATE UNIVERSITY SCIENCE AND TECHNOLOGY ----------------------------- TRAN THI THOM FINITE ELEMENT MODELS IN VIBRATION ANALYSIS OF TWO-DIMENSIONAL FUNCTIONALLY GRADED BEAMS Major: Mechanics of Solid code: 9440107 SUMMARY OF DOCTORAL THESIS IN MATERIALS SCIENCE Hanoi – 2019 The thesis has been completed at: Graduate University Science and Technology – Vietnam Academy of Science and Technology. Supervisors: 1. Assoc. Prof. Dr. Nguyen Dinh Kien 2. Assoc. Prof. Dr. Nguyen Xuan Thanh Reviewer 1: Prof. Dr. Hoang Xuan Luong Reviewer 2: Prof. Dr. Pham Chi Vinh Reviewer 3: Assoc. Prof. Dr. Phan Bui Khoi Thesis is defended at Graduate University Science and Technology- Vietnam Academy of Science and Technology at , on . Hardcopy of the thesis be found at : - Library of Graduate University Science and Technology - Vietnam national library 1PREFACE 1. The necessity of the thesis Publications on vibration of the beams are most relevant to FGM beams with material properties varying in one spatial direction only, such as the thickness or longitudinal direction. There are practical circumstances, in which the unidirectional FGMs may not be so appropriate to resist multi-directional variations of thermal and mechanical loadings. Optimiz- ing durability and structural weight by changing the volume fraction of FGM’s component materials in many different spatial directions is a mat- ter of practical significance, being scientifically recognized by the world’s scientists, especially Japanese researchers in recent years. Thus, struc- tural analysis with effective material properties varying in many different directions in general and the vibration of FGM beams with effective mate- rial properties varying in both the thickness and longitudinal directions of beams (2D-FGM beams) in particular, has scientific significance, derived from the actual needs. It should be noted that when the material properties of the 2D-FGM beam vary in longitudinal direction, the coefficients in the differential equation of beam motion are functions of spatial coordinates along the beam axis. Therefore analytical methods are getting difficult to analyze vibration of the 2D-FGM beam. Finite element method (FEM), with many strengths in structural analysis, is the first choice to replace traditional analytical methods in studying this problem. Developing the finite element models, that means setting up the stiffness and mass ma- trices, used in the analysis of vibrations of the 2D-FGM beam is a mat- ter of scientific significance, contributing to promoting the application of FGM materials into practice. From the above analysis, author has selected the topic: Finite element models in vibration analysis of two-dimensional functionally graded beams as the research topic for this thesis. 2. Thesis objective This thesis aims to develop finite element models for studying vibra- tion of the 2D-FGM beam. These models require high reliability, good convergence speed and be able to evaluate the influence of material pa- rameters, geometric parameters as well as being able to simulate the effect of shear deformation on vibration characteristics and dynamic responses of the 2D-FGM beam. 3. Content of the thesis 2Four main research contents are presented in four chapters of the the- sis. Specifically, Chapter 1 presents an overview of domestic and for- eign studies on the 1D and 2D-FGM beam structures. Chapter 2 pro- poses mathematical model and mechanical characteristics for the 2D- FGM beam. The equations for mathematical modeling are obtained based on two kinds of shear deformation theories, namely the first shear de- formation theory and the improved third-order shear deformation theory. Chapter 3 presents the construction of FEM models based on different beam theories and interpolation functions. Chapter 4 illustrates the nu- merical results obtained from the analysis of specific problems. Chapter 1. OVERVIEW This chapter presents an overview of domestic and foreign regime of re- searches on the analysis of FGM beams. The analytical results are dis- cussed on the basis of two research methods: analytic method and nu- merical method. The analysis of the overview shows that the numerical method in which FEM method is necessary is to replace traditional ana- lytical methods in analyzing 2D-FGM structure in general and vibration of the 2D-FGM beam in particular. Based on the overall evaluation, the thesis has selected the research topic and proposed research issues in de- tails. Chapter 2. GOVERNING EQUATIONS This chapter presents mathematical model and mechanical characteris- tics for the 2D-FGM beam. The basic equations of beams are set up based on two kinds of shear deformation theories, namely the first shear defor- mation theory (FSDT) and the improved third-order shear deformation theory (ITSDT) proposed by Shi [40]. In particular, according to ITSDT, basic equations are built based on two representations, using the cross- sectional rotation θ or the transverse shear rotation γ0 as an independent function. The effect of temperature and the change of the cross-section are also considered in the equations. 2.1. The 2D-FGM beam model The beam is assumed to be formed from four distinct constituent mate- rials, two ceramics (referred to as ceramic1-C1 and ceramic2-C2) and two metals (referred to as metal1-M1 and metal2-M2) whose volume fraction 3varies in both the thickness and longitudinal directions as follows: VC1 = ( z h + 1 2 )nz [ 1− ( x L )nx] VC2 = ( z h + 1 2 )nz ( x L )nx VM1 = [ 1− ( z h + 1 2 )nz][ 1− ( x L )nx] VM2 = [ 1− ( z h + 1 2 )nz]( x L )nx (2.1) Fig. 2.1 illustrates the 2D-FGM beam in Cartesian coordinate system (Oxyz). C1 0 M1 L, b, h Z C2 M2 X h b y z Fig. 2.1. The 2D-FGM beam model In this thesis, the effective material properties P (such as Youngs modulus, shear modulus, mass density, etc.) for the beam are evaluated by the Voigt model as: P =VC1PC1 +VC2PC2 +VM1PM1 +VM2PM2 (2.2) When the beam is in thermal environment, the effective properties of beams depend not only on the properties of the component materials but also on the ambient temperature. Then, one can write the expression for the effective properties of the beam exactly as follows: P(x,z,T ) = {[ PC1(T )−PM1(T ) ]( z h + 1 2 )nz +PM1(T) }[ 1− ( x L )nx] + {[ PC2(T )−PM2(T ) ]( z h + 1 2 )nz +PM2(T ) }( x L )nx (2.4) 4For some specific cases, such as nx = 0 or nz = 0, or C1 and C2 are identical, and M1 is the same as M2, the beam model in this thesis re- duces to the 1D-FGM beam model. Thus, author can verification the FEM model of the thesis by comparing with the results of the 1D-FGM beam analysis when there is no numerical result of the 2D-FGM beam. Its important to note that the mass density is considered to be temperature- independent [41]. The properties of constituent materials depend on temperature by a nonlinear function of environment temperature [125]: P = P0(P−1T−1 +1+P1T +P2T 2 +P3T 3) (2.7) This thesis studies the 2D-FGM beam with the width and height are linear changes in beam axis, means tapered beams, with the following three tapered cases [138]: Case A : A(x) = A0 ( 1− c x L ) , I(x) = I0 ( 1− c x L ) Case B : A(x) = A0 ( 1− c x L ) , I(x) = I0 ( 1− c x L )3 Case C : A(x) = A0 ( 1− c x L )2 , I(x) = I0 ( 1− c x L )4 (2.9) 2.2. Beam theories Based on the pros and cons of the theories, this thesis will use Timo- shenko’s first-order shear deformation theory (FSDT) [127] and the im- proved third-order shear deformation theory proposed by Shi (ITSDT) [40] to construct FEM models. 2.3. Equations based on FSDT Obtaining basic equations and energy expressions based on FSDT and ITSDT theory is similar, so Section 2.4 presents in more detail the process of setting up equations based on ITSDT. 2.4. Equations based on ITSDT 2.4.1. Expression equations according to θ From the displacement field, this thesis obtains expressions for strains and stresses of the beam. Then, the conventional elastic strain energy, UB 5is in the form UB = 1 2 L∫ 0 [ A11ε2m +2A12εmεb +A22ε2b −2A34εmεhs−2A44εbεhs +A66ε2hs +25 ( 1 16B11− 1 2h2 B22 + 1 h4 B44 ) γ20 ] dx (2.27) where A11, A12, A22, A34, A44, A66 and B11, B22, B44 are rigidities of beam and defined as: (A11, A12, A22, A34, A44, A66)(x,T) = ∫ A(x) E(x,z,T )(1, z, z2, z3, z4, z6)dA (B11, B22, B44)(x,T) = ∫ A(x) G(x,z,T )(1, z2, z4)dA (2.28) The kinetic energy of the beam is as follow: T = 1 2 L∫ 0 [ I11(u˙20+ w˙ 2 0)+ 1 2 I12u˙0(w˙0,x+5 ˙θ)+ 1 16 I22(w˙0,x +5 ˙θ)2 − 10 3h2 I34u˙0(w˙0,x + ˙θ)− 56h2 I44(w˙0 + ˙θ)(w˙0+5 ˙θ)+ 25 9h4 I66(w˙0,x + ˙θ)2 ] dx (2.29) in which (I11, I12, I22, I34, I44, I66)(x)= ∫ A(x) ρ(x,z) ( 1, z, z2, z3, z4, z6 ) dA (2.30) are mass moments. The beam rigidities and mass moments of the beam are in the follow- ing forms: Ai j = AC1M1i j − ( AC1M1i j −AC2M2i j )( x L )nx Bi j = BC1M1i j − ( BC1M1i j −B C2M2 i j )( x L )nx (2.31) 6with AC1M1i j , BC1M1i j are the rigidities of 1D-FGM beam composed of C1 and M1; AC2M2i j , BC2M2i j are the rigidities of 1D-FGM beam composed of C2 and M2. Noting that rigidities of 1D-FGM beam are functions of z only, the explicit expressions for this rigidities can easily be obtained. 2.4.2. Expression equations according to γ0 Using a notation for the transverse shear rotation (also known as clas- sic shear rotation), γ0 = w0,x+θ as an independent function, the axial and transverse displacements in (2.13) can be rewritten in the following form u(x,z, t) = u0(x, t)+ 1 4 z ( 5γ0−4w0,x ) − 5 3h2 z 3γ0 w(x,z, t) = w0(x, t) (2.35) Similar to the construction of basic equations according to θ , the thesis also receives basic equations expressed in γ0. 2.5. Initial thermal stress Assuming the beam is free stress at the reference temperature T0 and it is subjected to thermal stress due to the temperature change. The initial thermal stress resulted from a temperature ∆T is given by [18, 70]: σTxx =−E(x,z,T)α(x,z,T)∆T (2.41) in which elastic modulus E(x,z,T ) and thermal expansion α(x,z,T) are obtained from Eq.(2.4). The strain energy caused by the initial thermal stress σTxx has the form [18, 65]: UT = 1 2 L∫ 0 NTw20,xdx (2.42) where NT is the axial force resultant due to the initial thermal stress. σTxx: NT = ∫ A(x) σTxxdA =− ∫ A(x) E(x,z,T )α(x,z,T)∆TdA (2.43) The total strain energy resulted from conventional elastic strain energy UB, and strain energy due to initial thermal stress UT [70]. 2.6. Potential of external load 7The external load considered in the present thesis is a single moving constant force with uniform velocity. The force is assumed to cause bend- ing only for beams. The potential of this moving force can be written in the following form V =−Pw0(x, t)δ [ x− s(t) ] (2.44) where δ (.) is delta Dirac function; x is the abscissa measured from the left end of the beam to the position of the load P, t is current time calcu- lated from the time when the load P enters the beam, and s(t) = vt is the distance which the load P can travel. 2.7. Equations of motion In this section, author presents the equations of motion based on ITSDT with γ0 being the independent function. Motion equations for beams based on FSDT and ITSDT with θ is independent function that can be obtained in the same way. Applying Hamiltons principle, one obtained the motion equations system for the 2D-FGM beam placed in the temper- ature environment under a moving force as follows: I11u¨0 + 1 4 ( 5γ¨0−4w¨0,x ) I12− 5 3h2 I34γ¨0− [ A11u0,x + 1 4 A12 ( 5γ0,x−4w0,xx ) − 5 3h2 A34γ0,x ] ,x = 0 (2.51) I11w¨0 + [ I12u¨0 + 1 4 ( 5γ¨0−4w¨0,x ) I22− 5 3h2 I44γ¨0 ] ,x − [ A12u0,x + 1 4 A22 ( 5γ0,x−4w0,xx ) − 5 3h2 A44γ0,x ] ,xx = ( NTw0,x ) ,x −Pδ [ x− s(t) ] (2.52) 81 4 I12u¨0 + 1 16 I22 ( 5γ¨0−4w¨0,x ) − 1 3h2 I34u¨0− 1 3h2 I44 (5 2 γ¨0− w¨0,x ) + 5 9h4 I66γ¨0− [ 1 4 A12u0,x + 1 16A22 ( 5γ0,x−4w0,xx ) − 1 3h2 A34u0,x − 1 3h2 A44 (5 2 γ0,x−w0,xx ) − 5 9h4 A66γ0,x ] ,x +5 ( 1 16B11− 1 2h2 B22 + 1 h4 B44 ) γ0 = 0 (2.53) Notice that the coefficients in the system of differential equations of motion are the rigidities and mass moments of the beam, which are the functions of the spatial variable according to the length of the beam and the temperature, thus solving this system using analytic method is diffi- cult. FEM was selected in this thesis to investigate the vibration charac- teristics of beams. Conclusion of Chapter 2 Chapter 2 has established basic equations for the 2D-FGM beam based on two kinds of shear deformation theories, namely FSDT and ITSDT. The effect of temperature and the change of the cross-section is consid- ered in establishing the basic equations. Energy expressions are presented in detail for both FSDT and ITSDT in Chapter 2. In particular, with ITSDT, basic equations and energy expressions are established on the cross-sectional rotation θ or the transverse shear rotation γ0 as indepen- dent functions. The expression for the strain energy due to the tempera- ture rise and the potential energy expression of the moving force are also mentioned in this Chapter. Equations of motion for the 2D-FGM beam are also presented using ITSDT with γ0 as independent function. These energy expressions are used to obtain the stiffness matrices and mass ma- trices used in the vibration analysis of the 2D-FGM beam in Chapter 3. Chapter 3. FINITE ELEMENT MODELS This chapter builds finite element (FE) models, means that establish expressions for stiffness matrices and mass matrices for a characteristic element of the 2D-FGM beam. The FE model is constructed from the energy expressions received by using the two beam theories in Chapter 2. Different shape functions are selected appropriately so that beam ele- ments get high reliability and good convergence speed. Nodal load vector 9and numerical procedure used in vibration analysis of the 2D-FGM beam are mentioned at the end of the chapter. 3.1. Model of finite element beams based on FSDT This model constructed from Kosmatka polynomials referred as FBKo in this thesis can be avoided the shear-locking problem. In addition, this model has a high convergence speed and reliability in calculating the nat- ural frequencies of the beam. However, the FBKo model with 6 d.o.f has the disadvantage that the Kosmatka polynomials must recalculate each time the element mesh changes, thus time-consuming calculations. The FE model uses hierarchical functions, referred as FBHi model in the the- sis, which is one of the options to overcome the above disadvantages. Recently, hierarchical functions are used to develop the FEM model in 1D-FGM beam analysis (such as Bui Van Tuyen’s thesis). Based on the energy expressions received in Chapter 2, the thesis has built FBKo model and FBHi model using the Kosmatka function and hierarchical interpola- tion functions, respectively. The process of building FE models is similar, Section 3.2 will presents in detail the construction of stiffness and mass matrices for a characteristic element based on ITSDT. 3.2. Model of finite element beams based on ITSDT With two representations of the displacement field, two FEM models corresponding to these two representations will be constructed below. For convenience, in the thesis, FEM model uses the cross-sectional rotation θ as the independent function is called TBSθ model, FEM model uses the transverse shear rotation as the independent function is called TBSγ . 3.2.1. TBSθ model Different from the FE model based on FSDT, the vector of nodal dis- placements for two-node beam element (i, j), using the high order shear deformation theory in general and ITSDT in particular, has eight compo- nents: dSθ = {ui wi wi,x θi u j w j w j,x θ j}T (3.28) The displacements u0, w0 and rotation θ are interpolated from the nodal displacements as u0 = NudSθ , w0 = NwdSθ , θ = Nθ dSθ (3.29) where Nu, Nw and Nθ are, respectively, the matrices of shape functions 10 for u0, w0 and θ . Herein, linear shape functions are used for the axial displacement u0(x, t) and the cross-section rotation θ (x, t), Hermite shape functions are employed for the transverse displacement w0(x, t). With the interpolation scheme, one can write the expression for the de- formation components in the form of a matrix through a nodal displace- ments vector (3.28) as follows εSθm = u0,x = BSθm dSθ εSθb = 1 4 (5θ,x +w0,xx) = BSθb dSθ εSθhs = 5 3h2 (θ,x +w0,xx) = B Sθ hs dSθ εSθs = θ +w0,x = BSθm dSθ (3.33) In (3.33), the strain-displacement matrices BSθm , BSθb , BSθhs and BSθs are as follows BSθm = { − 1 l 0 0 0 1 l 0 0 0 } BSθb = 1 4 { 0 − 6l2 + 12x l3 − 4 l + 6x l2 − 5 l 0 6 l2 − 12x l3 − 2 l + 6x l2 5 l } BSθhs = 5 3h2 { 0 − 6l2 + 12x l3 − 4 l + 6x l2 − 1 l 0 6 l2 − 12x l3 − 2 l + 6x l2 1 l } BSθs = { 0 − 6xl2 + 6x2 l3 1− 4x l + 3x2 l2 l− x l 0 6x l2 − 6x2 l3 − 2x l + 3x2 l2 x l } (3.34) The elastic strain energy of the beam UB in Eq.(2.27) can be written in the form UB = 1 2 nE ∑(dSθ )T kSθ dSθ (3.9) where the element stiffness matrix kSθ is defined as kSθ = kSθm +kSθb +k Sθ s +kSθhs +k Sθ c (3.35) 11 in which kSθm = l∫ 0 ( BSθm )T A11 BSθm dx ; kSθb = l∫ 0 ( BSθb )T A22 BSθb dx kSθs = 25 l∫ 0 ( BSθs )T( 1 16B11− 1 2h2 B22 + 1 h4 B44 ) BSθs dx kSθhs = l∫ 0 ( BSθhs )T A66 BSθhs dx kSθc = l∫ 0 [( BSθm )T A12 BSθb − ( BSθm )T A34 BSθhs − ( BSθb )T A44 BSθhs ] dx (3.36) One write the kinetic energy in the following form T = 1 2 nE ∑( ˙dK)T m ˙dK (3.13) in which the element consistent mass matrix is in the form m = m11uu +m 12 uθ +m 22 θθ +m 34 uγ +m 44 θγ +m 66 γγ +m 11 ww (3.37) with m11uu = l∫ 0 NTu I11Nudx ; m12uθ = 1 4 l∫ 0 NTu I12(Nw,x + 5Nθ )dx m22θθ = l∫ 0 1 16(N T w,x + 5NTθ )I22(Nw,x + 5Nθ )dx ; m34uγ =− 5 3h2 l∫ 0 NTu I34(Nw,x +Nθ )dx m44θγ =− 5 12h2 l∫ 0 (NTw,x + 5NTθ )I44(Nw,x +Nθ )dx m66γγ = 25 9h4 l∫ 0 (NTw,x +NTθ )I66(Nw,x +Nθ )dx ; m11ww = l∫ 0 NTwI11Nwdx (3.38) are the element mass matrices components. 12 3.2.2. TBSγ model With γ0 is the independent function, the vector of nodal displacements for a generic element, (i, j), has eight components: dSγ = {ui wi wi,x γi u j w j w j,x γ j}T (3.39) The axial displacement, transverse displacement and transverse shear rotation are interpolated from the nodal displacements according to u0 = Nu dSγ , w0 = Nw dSγ , γ0 = Nγ dSγ (3.40) with Nu,Nw and Nγ are the matrices of shape functions for u0,w0 and γ0, respectively. Herein, linear shape functions are used for the axial displacement u0(x, t) and the transverse shear rotation γ0, Hermite shape functions are employed for the transverse displacement w0(x, t). The con- struction of element stiffness and mass matrices are completely similar to TBSθ model. 3.3. Element stiffness matrix due to initial thermal stress Using the interpolation functions for transverse displacement w0(x, t), one can write expressions for the strain energy due to the temperature rise (2.42) in the matrix form as follows UT = 1 2 nE ∑dT kTd (3.44) where kT = l∫ 0 BTt NTBtdx (3.45) is the stiffness due to temperature rise. For different beam theories, the element stiffness matrix due to temperature rise has the same form (3.45). The only difference is that the differenc