Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Analysis of Thermo-Piezoelectricity Problems by Meshless Method

Analysis of Thermo-Piezoelectricity Problems by Meshless Method In this paper meshless method based on the local Petrov-Galerkin approach is pesented for the solution of boundary value problems for coupled thermo-electro-mechanical fields. Transient dynamic governing equations are considered in analysis of the problems. Material properties of piezoelectric materials are influenced by a thermal field. It is leading to an induced nonhomogeneity and the governing equations are more complicated compared to a homogeneous counterpart. Twodimensional analyzed domain is divided into small circular subdomains surrounding nodes that are randomly spread over the whole domain. A unit ep function is used as the te functions in the local weak-form. The derived local integral equations (LIEs) have boundary-domain integral form. The moving lea-uares (MLS) method is adopted for the approximation of the physical quantities in the LIEs and afterwards to obtain a syem of ordinary differential equations (ODE) for unknown nodal quantities. To solve this syem of ODE, Houbolt finite-difference scheme is applied as a time-epping method. INTRODUCTION Many engineering applications, smart ructures and devices take advantage of piezoelectric materials. Piezoelectric materials are often referred to as smart materials. They are extensively utilized as transducers, sensors and actuators in many engineering fields. Piezoelectric materials have usually anisotropic properties. Except this complication electric and mechanical fields are coupled each other and the governing equations are much more complex than those in the classical elaicity. In order to solve the boundary or the initial-boundary value problems for piezoelectric solids, efficient computational methods are required. Moly, the finite element method (FEM) as shown in [9], [11], [13] and the boundary element method (BEM) [7], [10], [14] and [25] are applied to solve general piezoelectric problems. Meshless methods have been also successfully applied to piezoelectric problems in [16] and [20]. In certain piezoelectric materials pyroelectric effect can be observed. This means they are temperature sensitive, i.e. an electric charge or voltage is generated when exposed to temperature variations. Thus a coupling of thermo-electro-mechanical fields is needed to be taken into account if a temperature load is considered in a piezoelectric solid. Mindlin [17] for the fir time introduced the theory of thermo-piezoelectricity. The physical laws for thermo-piezoelectric materials have been explored by Nowacki [19]. In the work of Dunn [8], micromechanics models for effective thermal expansion and pyroelectric coefficients of piezoelectric composites were udied. Qin [22] offered a review on fracture of thermo-piezoelectric materials. Boundary value problems for coupled fields are complex, thus analytical methods can be only applied to simple problems of thermo-piezoelectricity, e.g. [26], [27] and [32]. The analysis and design process of smart engineering ructures with integrated piezoelectric actuators or sensors requires powerful calculation tools. Up to now the finite element methods (FEM) provides an effective technique, as shown in [33] in a homogeneous medium. Rao and Sunar [24] inveigated the piezothermoelectric problem of intelligent ructures with diributed piezoelectric sensors and actuators and concluded that the inclusion of the thermal effects may help improve the performance characteriics of the syem. Material properties under a thermal load are influenced by temperature. However, mo inveigations in piezothermoelaicity were done under the assumption of the temperature independent material properties. Bert and Birman [5] showed that the piezoelectric coefficients are ress and electric-field dependent. If this phenomenon is considered the material properties are continuously varying with Cartesian coordinates. In this so called induced nonhomogeneity the governing equations are more complicated than in a homogeneous counterpart. Some relative simple problems of coupled electro-mechanical fields in continuously nonhomogeneous solids have been successfully solved in previous works [34] and [35]. Aouadi [1] made the fir attempt to solve induced nonhomogeneity problem in thermo-piezoelectricity for an infinite and half-space. A general numerical solution for induced nonhomogeneity problem in thermo-piezoelectricity is not available according to the be of the author's knowledge. The meshless local Petrov-Galerkin (MLPG) method is considered as a fundamental base for the derivation of many meshless formulations, since trial and te functions can be chosen from different functional spaces. The MLPG method introduced by Atluri et al. [2], Atluri [3] and Sladek et al. [28], using a Heaviside ep function as the te functions, has been applied to solve 2-D homogeneous piezoelectric problems in paper by Sladek et al. [29]. Recently, meshless method (MLPG) was applied to analyze continuously nonhomogeneous piezoelectric solids under a mechanical or electrical load [30]. Extension of the MLPG method to 2-D thermo-piezoelectric solids with induced nonhomogeneity is given in the present paper. The coupled thermo-electro-mechanical fields in thermo-piezoelectricity are described by partial differential equations, mechanical fields are described by the equations of motion with an inertial term. Maxwell´s equation for the electrical field has a quasi-ic character and thermal field is described by the heat conduction equation, which has a diffusive character. Nodal points are introduced and spread on the analyzed domain and each node is surrounded by a small circle for simplicity, but without loss of generality. Numerical integrations over simple shape of subdomains like circles can be easily carried out. The integral equations have a very simple nonsingular form. The spatial variations of the displacements, the electric potential and the temperature are approximated by the Moving Lea-uares (MLS) scheme [4]. After performing the spatial integrations, a syem of ordinary differential equations for the unknown nodal values is obtained. The boundary conditions on the global boundary are satisfied by the collocation of the MLS-approximation expressions for the displacements, the electric potential and the temperature at the boundary nodal points. Finally, the syem of ordinary differential equations is solved by the Houbolt finite-difference scheme [12]. Several numerical examples are introduced to verify the accuracy and the efficiency of the proposed MLPG method. GOVERNING EQUATIONS FOR THERMO-ELECTROMECHANICAL FIELDS The governing equations for thermo-piezoelectricity in continuously nonhomogeneous solids under the quasi-electroic assumption are given according to Mindlin [18] by the equation of motion for displacements, the fir Maxwell`s equation for the vector of electric displacements and heat conduction equation: ij , j (x, ) X i (x, ) (x)ui (x, ) D j , j ( x, ) R ( x, ) 0 (1) (2) kij (x), j (x, ) ,i (3) (x)c(x) (x, ) S (x, ) 0 vij, x, i, ui, Di, Xi, R and S are ress, time, temperature difference, displacement, electric displacement, density of body force vector, volume density of free charges and density of heat sources, respectively. Adding that t, kij and c are the mass density, thermal conductivity tensor and specific heat, respectively. The dots over a quantity indicate the time derivatives. A ic problem can be considered formally as a special case of the dynamic one, by omitting the acceleration üi(x,x) in the equations of motion (1) and the time derivative terms in equation (3). Both cases are analyzed here. The heat generation by mechanical and electrical fields can be neglected for mo materials because the inverse thermoelaic and pyroelectric effects are very weak. Assuming this, the conitutive relations representing the partial-coupling of the mechanical, electrical and thermal fields are: ij (x, ) cijkl (x) kl (x, ) ekij (x) Ek (x, ) ij (x) (x, ) D j (x, ) e jkl (x) kl (x, ) (4) (5) h jk (x) Ek (x, ) p j (x) (x, ) 18 VOLUME 14, No. 4, 2010 cijkl(x), ejkl(x), hjk(x) and pj(x) are the elaic, piezoelectric, dielectric and pyroelectric material tensors in a continuously nonhomogeneous piezoelectric medium, respectively. The ress-temperature modulus cij(x) can be expressed through the iffness coefficients and the coefficients of linear thermal expansion akl For the mechanical field following essential and natural boundary conditions are assumed: u ui (x, x ) = ui (x, x ) on Cu, u ti (x, x ) = vij n j = t i (x, x ) on Ct, for the electrical field: ij cijkl kl (6) u u } (x) = } (x) on C p, ni Di (x) = Q (x) on Cq, and finally for the thermal field: The rain tensor fij and the electric field vector Ej are related to the displacements ui and the electric potential } by following expressions 1 ij ui , j u j ,i 2 E j , j u i (x, x ) = i (x, x ) on Ca, u q (x, x ) = kij i,j (x, x ) ni (x) = q (x, x ) on Cb Cu is the part of the global boundary with prescribed displacements, and on Ct, Cp, Cq, Ca, and Cb the traction vector, the electric potential, the surface charge density, the temperature and the heat flux are prescribed, respectively. (7) (8) Considering plane problems, the conitutive equations are frequently written in terms of the second-order tensor of elaic conants [15]. Many piezoelectric solids are transversely isotropic. Under the plane rain condition with f33=f31=f32= E3=0, the conitutive equations (4) and (5) are reduced in this case to FORMULATION OF LOCAL INTEGRAL EQUATIONS The MLPG method is applied to conruct a weakform over the local fictitious subdomains Xs, which is a small region taken for each node inside the global domain, as shown in [3]. The local subdomains overlap each other, and cover the whole global domain X. The local subdomains could be of any geometrical shape and size. In the present paper, the local subdomains are considered to be of the circular shape. The local weak-form of the governing equation (1) can be written as 11 c11 c 22 12 12 0 c12 c22 0 0 11 0 e21 E 0 22 0 e22 1 E c66 212 e15 0 2 11 E C(x) 22 L(x) 1 E2 212 D1 0 D e 2 21 0 e22 11 e15 h 22 11 0 2 0 12 (9) ij , j (x, ) (x)ui (x, ) 0 E1 h22 E2 * X i (x, ) uik (x) d 0 (11) 11 p1 E G (x) 22 H (x) 1 P(x) p2 E2 212 u ik (x) is a te function. Applying the Gauss divergence theorem to eq. (11) and rearranging, one can obtain (10) ij * * (x, )n j (x)uik (x)d ij (x, )uik , j (x)d s * ik X (x, ) (x)u (x, ) u i i ( x) d 0 (12) c11 c12 0 c12 c22 0 c13 11 11 c23 22 22 0 33 0 uXs is the boundary of the local subdomain which consis of three parts uXs = Ls , C , Csu [3]. Here Ls is the local boundary that is totally inside the global domain, C is the part of the local boundary which coincides with the global traction boundary, i.e., C = uXs + Ct and similarly Csu is the part of the local boundary that coincides with the global displacement boundary, i.e. Csu = uXs + Cu. Heaviside ep function is chosen as the te func* tion u ik (x) in each subdomain at x s * uik (x) ik 0 at x s w* (x) is a te function. Applying the Gauss divergence theorem to the local weak-form (16) and considering the Heaviside ep function for the te function w* (x) one can obtain Ls sa (x)c(x)(x, )d (17) and using this, the local weak-form (12) is converted into the following local boundary-domain equations Ls su t i ( x, ) d (x)u (x, )d ti (x, )d (13) It is important to note that the local integral equations (13) are valid for both homogeneous and nonhomogeneous linear piezoelelectric solids. Nonhomogeneous material properties are included in eq. (13) through the elaic and piezoelectric tensor of the material coefficients in the traction components. The local weak form of the governing equation (2) is given in the similar way as There is no requirement in the MLPG method on the te and the trial functions to be necessarily from the same functional spaces. For internal nodes, the te function is chosen as the Heaviside ep function with its support on the local subdomain. The trial functions, on the other hand, are chosen to be the moving lea-uares (MLS) approximation over a number of nodes spread within the domain of influence. The approximated functions for the mechanical displacements, the electric potential and the temperature can be written according to Atluri [4] as: t t u h (x, x ) = UT (x) $ u ( x ) = / za (x) u a ( x ) ^ h (x, ) a (x) a ( ) n j, j ( x, ) R ( x, ) v * ( x ) d 0 (14) ^ h (x, ) a (x) a ( ) (18) o* (x) is a te function. Again, as in previous case, applying the Gauss divergence theorem to the local weak-form (14) and considering the Heaviside ep function for the te function o* (x), one will obtain Ls sp Q ( x, ) d Q ( x, ) d R ( x, ) d (15) Q(x,x)=Dj (x,x)nj. The local weak-form of the diffusion equation (3) can be written as ij (x), j (x, ) (x)c(x) (x, ) ,i (16) t t t the nodal values ua(x), } a(x) and i a(x) are fictitious parameters for the displacements, the electric potential and the temperature, respectively, and za(x) is the shape function associated with the node a. The number of nodes n used for the approximation is determined by the weight function wa(x). In the MLS approximation a 4th order spline-type weight function is applied [3]. It can be easily shown, that the C1 - continuity is ensured over the entire domain, therefore the continuity conditions of the tractions, the electric charge and the heat flux will be satisfied. The traction vector ti(x,x) at a boundary point xduXS is approximated in terms of the same nodal t values ua(x) as ^ t h (x, ) N(x)C(x) B a (x)u a ( ) n W (x, ) w (x) d 0 20 VOLUME 14, No. 4, 2010 ^ N(x)L(x) A a (x) a ( ) t u / z ( g ) } ( x ) = } ( g , x ) for g ! C a a n (23) (24) ^ N(x) a (x) a ( ) (19) t u / z ( g ) i ( x ) = i ( g , x ) for g !C a a the matrix N(x) is related to the normal vector n(x) on uXS by n 0 n2 N ( x) 1 0 n2 n1 Discretized forms of the unknown quantities in the local boundary-domain integral equations (13), (15) and (17) are obtained in the view of the MLSapproximation (19) - (21) as and the matrices Ba and Aa are represented by the gradients of the shape functions as ^ N(x)C(x)B a (x)d u a ( ) L s su z,a1 0 z,a1 B a (x) = > 0 z,a2 H, A a (x) = < a F z,2 z,a2 z,a1 The electrical charge Q(x,x) can be approximated in similar way by n ^ a (x)d u a ( ) s n ^ N(x)L(x)A a (x)d a ( ) L s su ^ Q h (x, ) N1 (x)G (x) B a (x)u a ( ) -/ d Ls + Csu t # N (x) c (x) z (x) dC n i ( x ) = a a ^ N1 (x)H (x) A a (x) a ( ) t ( x, ) d (25) ^ N1 (x)P(x) a (x) a ( ) (20) the matrices G and H are defined in eq. (10) and N1(x)= [n1 n2]. The heat flux q(x,x) is approximated by ^ q h (x, ) kij ni ,aj (x) a ( ) ^ N1 (x)K (x) A a (x) a ( ) (21) ^ N1 (x)G (x)B a (x)d u a ( ) Ls sp n ^ N1 (x)H (x)A a (x)d a ( ) Ls sp n ^ N1 (x)P(x) a (x)d a ( ) L s sp Q ( x, ) d R ( x, ) d (26) k k K (x) 11 12 k12 k22 Using the approximation formula (18) and fulfilling the essential boundary conditions at those nodal points on the global boundary, the displacements, the electrical potential and the temperature are prescribed, discretized form of the boundary conditions is obtained as ^ N1 (x)K (x)A a (x)d a ( ) Ls sa n ^ c a (x)d a ( ) s (27) t u / z ( g ) u ( x ) = u ( g , x ) for g ! C a a (22) Equations (25), (26) and (27) are considered on the sub-domains adjacent to interior nodes as well as to the boundary nodes on C, C and C. Complete syem of ordinary differential equations for the computation of the nodal unknowns can be obtained by collecting the discretized local boundary-domain integral equations together with the discretized boundary conditions for the displacements, the electrical potential and the temperature. t As nodal unknowns the fictitious parameters u a(x), t t } a(x) and i a(x) are considered. aa General variation of material properties with Cartesian coordinates was assumed in previous formulations. If material parameters are dependent on temperature, one can write for a general material parameter A: A(x) f ( (x)) ^ N1 (x)G ( k 1) (x)B a (x)d u a ( k ) ( ) Ls sp n ^ N1 (x)H ( k 1) (x)A a (x)d a ( k ) ( ) Ls sp n ^ N1 (x)P ( k 1) (x) a (x)d a ( k ) ( ) Ls sp ( x, ) d R ( x, ) d Q (32) (28) ^ N1 (x)K ( k 1) (x)A a (x)d a ( k ) ( ) Ls sa Using equation (28) to replace material parameters in conitutive equations (4) and (5), a nonlinear expression is obtained. This results in necessity of application of iterative approach in each time ep. The linearized conitutive equations in the k-th iteration ep are then given as ( (k ( ijk ) (x, ) cijkl1) ( ) klk ) (x, ) (k ( ekij 1) ( ) Ek( k ) (x, ) ijk 1) ( ) ( k ) (x, ) k ( D (j k ) (x, ) e(jkl1) ( ) klk ) (x, ) k h (jk 1) ( ) Ek( k ) (x, ) p (jk 1) ( ) ( k ) (x, ) n ^ (x) c(x) a (x)d a ( k ) ( ) s (33) (29) (30) Applying now equations (29) and (30) ed above, we obtain a set of ordinary differential equations for the k-th iteration ep in form ^ N(x)C( k 1) (x)B a (x)d u a ( k ) ( ) Ls su n n ^ (x) a (x)d u a ( k ) ( ) s Material parameters are taken at a reference temperature for the 1 iteration ep. The iteration process is ended if the difference of the Sobolev-norms for temperatures in two successive eps is smaller than a prescribed tolerance. Now, we should explain solution of the complete syem of ordinary differential equations. The syem of above given local integral equations can be rearranged in such a way that all known quantities are on the r.h.s. Thus, in matrix form the syem becomes Ax Bx Cx Y (34) In the present work, the Houbolt method is applied. In the Houbolt finite-difference scheme [12], the "acceleration" (ü=) is expressed as ^ N(x)L( k 1) (x)A a (x)d a ( k ) ( ) L s su 2x 5x 4x x 2 2 (35) -/ d # N (x) c Ls + Csu ^ k - 1h t (x) z (x) dC n ia(k) ( x ) = Dx is the time ep. The backward difference method is applied for the approximation of "velocities" t ( x, ) d (31) x x (36) Subituting eqs. (35) and (36) into eq. (34), we get 22 VOLUME 14, No. 4, 2010 the following syem of algebraic equations for the unknowns xx+Dx 1 1 2 2 A B C x 2 (5A B )x 1 A 2 4x x 2 Y (37) The value of the time ep has to be appropriately selected with respect to material parameters (wave velocities) and time dependence of the boundary conditions. te the present computational method. The material coefficients of the panel are considered for PZT5H according to Qin and Mai [23]: c11=12,6.1010Nm-2, c12=5,3.1010Nm-2, c22=11,7.1010Nm-2, c44=3,53.1010Nm-2, e15=17Cm-2, e21=-6,5Cm-2, e22=23,3Cm-2, h11=15,1.10-9C(Vm)-1, h22=13.10-9C(Vm)-1, t=7500kg/m3, k11=50W/Km, k22=75W/Km, a11=0,88.10-51/K, a22=0,5.10-51/K, p1=0, p2=-5,4831.10-6C/Km2, c=420Wskg-1K-1. Plane rain condition is assumed in the analysis. The mechanical displacement, the electrical potential and the thermal field on the finite uare panel with a size a x m x 1m are approximated by using 121 (11x11) equaly-spaced nodes. The circular local sub-domains are considered, each with a radius rloc=0,08. The Sobolev-norm is calculated for the purpose of error analysis. The relative error of the temperature in the considered time interval [0,T] is defined as NUMERICAL EXAMPLES Numerical results are presented in this section for selected problems of thermoelaicity in piezoelectric materials. In order to te the accuracy of the present meshless method a unit uare panel under a sudden heating on the top side is analyzed as the fir example (Fig. 1). The following analytical solution is available for uncoupled thermoelaicity in an isotropic material [6]: ( x2 , ) 1 exp 2n 1 n0 (1) (2n 1) 4a 2 num exact exact (39) (2n 1) x2 cos 2a (38) 2 d and time T is defined through normalized parameter Tk22/tca2=1,3. Fig. 2 presents numerical results for the temperature at the bottom side and the mid-line of the panel. The temperature is normalized by the intensity of the thermal shock i0=0. Obtained results are compared with the analytical results and an excellent agreement is observed. The relative error of the temperature, r, at both lines is less than 0,5%. For the total number of 441 nodes, the relative error r=0,15% has been obtained. The temperature diribution is not influenced by mechanical and electrical fields. Numerical results for the temporal variation of the direct ress v11 are presented in Fig. 3, an excellent agreement can be seen again between the presented and FEM results computed at both considered lines. The FEM results have been obtained by ANSYS code with 3600 quadratic eight-node elements and 1000 time increments. The resses are normalized 1/ 2 a is the side-length of the panel and l=k22/ tc is the diffusivity coefficient. x2 Q=0 t1=t2=0 , =H(-0), Q=0 q= 0 u1=t2=0 Q=0 q=0 22 1 u2=t1=0 ,q=0 , =0 x1 Fig. 1 Uniformly heated piezoelectric uare plate Homogeneous material properties are selected to by the ic value v11 =-1,2959.106Pa. MLPG: x2 /a = 0.5 0.2 0 -0.2 0 0.2 0.4 0.6 0.8 x2/a = 0. analytical x2/a = 0.5 x2/a = 0. k22/ca2 Fig. 2 Time variation of the temperature on two different lines parallel to x1- axis 11/11 MLPG: x2 /a = 0.5 x2/a = 0. FEM: x2/a = 0.5 x2/a = 0. The temporal variation of the electrical potential W is shown in Fig. 4. The potential is normalized by the ic value W=1,4861.104V at the top side of the panel if the pyroelectric material vector QS =(p1,p2) is vanishing. For p2=-5,4831.10-6C/Km2 we have obtained the ic value W=1,455.104V. One can observe that the pyroeletric parameter is only slightly decreasing the electrical potential in the considered problem. In order to control the mechanical quantities (displacements and resses) the inverse piezoelectric effect can be utilized. So as a next load case, we have applied a uniform electrical displacement D0 on the top side of the panel additionally to the thermal load. ionary boundary conditions are considered. In Fig. 5 one can observe the influence of the electrical displacement on the ress component v11. The increasing electrical displacement on the top side of the panel significantly reduces the normal ress. The electrical displacement is expressed through the normalized quantity m = D0e22 / h22c11a11i0. Negligible influence of the pyroelectric parameters on the ress values can be observed. 0 -0.2 p2 =0. p2 = -5.48e-6 -0.2 0 0.2 0.4 0.6 0.8 11/11c110 -0.4 -0.6 -0.8 -1 -1.2 k22/ca Fig. 3 Time variation of the v11 ress -1.4 0.00E+00 5.00E-05 1.00E-04 1.50E-04 Fig. 5 Variation of the v11 ress with the normalized electrical displacements MLPG: x2 /a = 0.5 x2/a = 1. FEM: x2/a = 0.5 x2/a = 1. -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 k22/ca2 Fig. 4 Time variation of the electric potential at two different lines parallel to x1- axis As the next example, an edge crack in a finite piezoelectric rip is analyzed. The geometry of the rip is given in Fig. 6 with the following values: a=0,5m, a/w=0,4 and h/w=1,2. Only a half of the rip is modeled due to the symmetry of the problem with respect to the x1 - axis. We have used 930 nodes equidiantly diributed for the MLS approximation of the physical quantities. On the both lateral sides of the rip thermal load is applied. 24 VOLUME 14, No. 4, 2010 x2 900 t2=t1=0 , Q=0 , q=0 =t Q=0 t1=t2=0 =0 Q=0 t1=t2=0 62 q=0 w =0 u2=t1=0 31 32 1 t1=t2=0 Q=0 a x1 Fig. 6 Edge crack in a finite rip under a thermal shock, half of the rip depicted due to a symmetry For cracks in homogeneous and linear piezoelectric solids the asymptotic behaviour of the field quantities has been given by Sosa [31] and Pak [21]. Using these expressions, one can obtain the following expression for the intensity factors [10] as KI u Re() 1 2 lim r 0 8r K IV the crack line because of the immediate electromechanical interaction. In the dynamic case the ress field is coupled not only to the immediate electric field, but also to inertia forces [9]. Numerical results for the normalized ress intensity factor fI=KI/ ra a11c11i0 are presented in Fig. 7. On the right lateral side of the finite specimen a vanishing heat flux is considered in the transient heat conduction case. We have also considered an induced nonhomogeneity caused by the dependence of the thermal expansion coefficient on the temperature, akl= akl0(1-mi), m=0,005K-1. It can be well observed that for a large time inant the SIF is vanishing in Fig. 7. The electrical displacement intensity factor (EDIF) for considered boundary value problem is vanishing since the thermal processes changes are relatively slow. FEM 1/2 MLPG: homog. induced nonhom. fI=KI/ 11c110(a) (40) the matrix B is determined by the material properties in paper by Garcia-Sanchez et al. [10]. Note that the symmetry conditions of the displacements and the potential with respect to the crack plane are utilized. Contrary to the ic case, the SIF is not vanishing for the thermal shock prescribed on the left lateral side of the rip. From the Maxwell`s equations, it is known that the velocity of electromagnetic waves is equal to the speed of light, which is much greater than the velocity of elaic and thermal waves. Hence, the use of quasi-ic approximation in eqs. (1) to (3) is juified for the interaction of electrical, mechanical and thermal fields. The response of the electric fields is immediate, while that of the elaic and thermal ones are taken as finite because of the finite velocity of both waves. On the other hand, in a ic case, the response of all the mechanical (rain, ress), thermal and electrical fields is immediate. Thus, the SIF is vanishing in such a case since the resses are zero ahead the crack-tip on k22/ca Fig. 7 Time variation of the ress intensity factor for an edge crack CONCLUSION A meshless local Petrov-Galerkin method (MLPG) is proposed for the solution of boundary value problems for coupled thermo-electro-mechanical fields. Transient dynamic governing equations are considered. The material properties of a piezoelectric material are influenced by a thermal field. It is leading to an induced nonhomogeneity and so to mathematical complications of the governing equations compared to a homogeneous counterpart. The analyzed domain is divided into small overlapping circular subdomains. A unit ep function is used as the te functions in the local weak-form of the governing partial differential equations. For the approxi- mation of the physical field quantities the Moving Lea-uares (MLS) scheme is used. The proposed method is a truly meshless method, which requires neither domain elements nor background cells in either the interpolation or the integration. The present method represents an alternative numerical tool for many exiing computational methods like the FEM or the BEM. The main advantage of the present method is its simplicity. The present formulation is promising for numerical analyses of multi-field problems like piezoelectric or thermoelaic problems. ACKNOWLEDGEMENT The authors gratefully acknowledge the support by the Slovak Science and Technology Assiance Agency regiered under number APVV-0427-07 and the Slovak Grant Agency VEGA-2/0039/09. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Acta Mechanica Slovaca de Gruyter

Analysis of Thermo-Piezoelectricity Problems by Meshless Method

Loading next page...
 
/lp/de-gruyter/analysis-of-thermo-piezoelectricity-problems-by-meshless-method-DpvBEkFzXi

References (44)

Publisher
de Gruyter
Copyright
Copyright © 2010 by the
ISSN
1335-2393
DOI
10.2478/v10147-011-0030-z
Publisher site
See Article on Publisher Site

Abstract

In this paper meshless method based on the local Petrov-Galerkin approach is pesented for the solution of boundary value problems for coupled thermo-electro-mechanical fields. Transient dynamic governing equations are considered in analysis of the problems. Material properties of piezoelectric materials are influenced by a thermal field. It is leading to an induced nonhomogeneity and the governing equations are more complicated compared to a homogeneous counterpart. Twodimensional analyzed domain is divided into small circular subdomains surrounding nodes that are randomly spread over the whole domain. A unit ep function is used as the te functions in the local weak-form. The derived local integral equations (LIEs) have boundary-domain integral form. The moving lea-uares (MLS) method is adopted for the approximation of the physical quantities in the LIEs and afterwards to obtain a syem of ordinary differential equations (ODE) for unknown nodal quantities. To solve this syem of ODE, Houbolt finite-difference scheme is applied as a time-epping method. INTRODUCTION Many engineering applications, smart ructures and devices take advantage of piezoelectric materials. Piezoelectric materials are often referred to as smart materials. They are extensively utilized as transducers, sensors and actuators in many engineering fields. Piezoelectric materials have usually anisotropic properties. Except this complication electric and mechanical fields are coupled each other and the governing equations are much more complex than those in the classical elaicity. In order to solve the boundary or the initial-boundary value problems for piezoelectric solids, efficient computational methods are required. Moly, the finite element method (FEM) as shown in [9], [11], [13] and the boundary element method (BEM) [7], [10], [14] and [25] are applied to solve general piezoelectric problems. Meshless methods have been also successfully applied to piezoelectric problems in [16] and [20]. In certain piezoelectric materials pyroelectric effect can be observed. This means they are temperature sensitive, i.e. an electric charge or voltage is generated when exposed to temperature variations. Thus a coupling of thermo-electro-mechanical fields is needed to be taken into account if a temperature load is considered in a piezoelectric solid. Mindlin [17] for the fir time introduced the theory of thermo-piezoelectricity. The physical laws for thermo-piezoelectric materials have been explored by Nowacki [19]. In the work of Dunn [8], micromechanics models for effective thermal expansion and pyroelectric coefficients of piezoelectric composites were udied. Qin [22] offered a review on fracture of thermo-piezoelectric materials. Boundary value problems for coupled fields are complex, thus analytical methods can be only applied to simple problems of thermo-piezoelectricity, e.g. [26], [27] and [32]. The analysis and design process of smart engineering ructures with integrated piezoelectric actuators or sensors requires powerful calculation tools. Up to now the finite element methods (FEM) provides an effective technique, as shown in [33] in a homogeneous medium. Rao and Sunar [24] inveigated the piezothermoelectric problem of intelligent ructures with diributed piezoelectric sensors and actuators and concluded that the inclusion of the thermal effects may help improve the performance characteriics of the syem. Material properties under a thermal load are influenced by temperature. However, mo inveigations in piezothermoelaicity were done under the assumption of the temperature independent material properties. Bert and Birman [5] showed that the piezoelectric coefficients are ress and electric-field dependent. If this phenomenon is considered the material properties are continuously varying with Cartesian coordinates. In this so called induced nonhomogeneity the governing equations are more complicated than in a homogeneous counterpart. Some relative simple problems of coupled electro-mechanical fields in continuously nonhomogeneous solids have been successfully solved in previous works [34] and [35]. Aouadi [1] made the fir attempt to solve induced nonhomogeneity problem in thermo-piezoelectricity for an infinite and half-space. A general numerical solution for induced nonhomogeneity problem in thermo-piezoelectricity is not available according to the be of the author's knowledge. The meshless local Petrov-Galerkin (MLPG) method is considered as a fundamental base for the derivation of many meshless formulations, since trial and te functions can be chosen from different functional spaces. The MLPG method introduced by Atluri et al. [2], Atluri [3] and Sladek et al. [28], using a Heaviside ep function as the te functions, has been applied to solve 2-D homogeneous piezoelectric problems in paper by Sladek et al. [29]. Recently, meshless method (MLPG) was applied to analyze continuously nonhomogeneous piezoelectric solids under a mechanical or electrical load [30]. Extension of the MLPG method to 2-D thermo-piezoelectric solids with induced nonhomogeneity is given in the present paper. The coupled thermo-electro-mechanical fields in thermo-piezoelectricity are described by partial differential equations, mechanical fields are described by the equations of motion with an inertial term. Maxwell´s equation for the electrical field has a quasi-ic character and thermal field is described by the heat conduction equation, which has a diffusive character. Nodal points are introduced and spread on the analyzed domain and each node is surrounded by a small circle for simplicity, but without loss of generality. Numerical integrations over simple shape of subdomains like circles can be easily carried out. The integral equations have a very simple nonsingular form. The spatial variations of the displacements, the electric potential and the temperature are approximated by the Moving Lea-uares (MLS) scheme [4]. After performing the spatial integrations, a syem of ordinary differential equations for the unknown nodal values is obtained. The boundary conditions on the global boundary are satisfied by the collocation of the MLS-approximation expressions for the displacements, the electric potential and the temperature at the boundary nodal points. Finally, the syem of ordinary differential equations is solved by the Houbolt finite-difference scheme [12]. Several numerical examples are introduced to verify the accuracy and the efficiency of the proposed MLPG method. GOVERNING EQUATIONS FOR THERMO-ELECTROMECHANICAL FIELDS The governing equations for thermo-piezoelectricity in continuously nonhomogeneous solids under the quasi-electroic assumption are given according to Mindlin [18] by the equation of motion for displacements, the fir Maxwell`s equation for the vector of electric displacements and heat conduction equation: ij , j (x, ) X i (x, ) (x)ui (x, ) D j , j ( x, ) R ( x, ) 0 (1) (2) kij (x), j (x, ) ,i (3) (x)c(x) (x, ) S (x, ) 0 vij, x, i, ui, Di, Xi, R and S are ress, time, temperature difference, displacement, electric displacement, density of body force vector, volume density of free charges and density of heat sources, respectively. Adding that t, kij and c are the mass density, thermal conductivity tensor and specific heat, respectively. The dots over a quantity indicate the time derivatives. A ic problem can be considered formally as a special case of the dynamic one, by omitting the acceleration üi(x,x) in the equations of motion (1) and the time derivative terms in equation (3). Both cases are analyzed here. The heat generation by mechanical and electrical fields can be neglected for mo materials because the inverse thermoelaic and pyroelectric effects are very weak. Assuming this, the conitutive relations representing the partial-coupling of the mechanical, electrical and thermal fields are: ij (x, ) cijkl (x) kl (x, ) ekij (x) Ek (x, ) ij (x) (x, ) D j (x, ) e jkl (x) kl (x, ) (4) (5) h jk (x) Ek (x, ) p j (x) (x, ) 18 VOLUME 14, No. 4, 2010 cijkl(x), ejkl(x), hjk(x) and pj(x) are the elaic, piezoelectric, dielectric and pyroelectric material tensors in a continuously nonhomogeneous piezoelectric medium, respectively. The ress-temperature modulus cij(x) can be expressed through the iffness coefficients and the coefficients of linear thermal expansion akl For the mechanical field following essential and natural boundary conditions are assumed: u ui (x, x ) = ui (x, x ) on Cu, u ti (x, x ) = vij n j = t i (x, x ) on Ct, for the electrical field: ij cijkl kl (6) u u } (x) = } (x) on C p, ni Di (x) = Q (x) on Cq, and finally for the thermal field: The rain tensor fij and the electric field vector Ej are related to the displacements ui and the electric potential } by following expressions 1 ij ui , j u j ,i 2 E j , j u i (x, x ) = i (x, x ) on Ca, u q (x, x ) = kij i,j (x, x ) ni (x) = q (x, x ) on Cb Cu is the part of the global boundary with prescribed displacements, and on Ct, Cp, Cq, Ca, and Cb the traction vector, the electric potential, the surface charge density, the temperature and the heat flux are prescribed, respectively. (7) (8) Considering plane problems, the conitutive equations are frequently written in terms of the second-order tensor of elaic conants [15]. Many piezoelectric solids are transversely isotropic. Under the plane rain condition with f33=f31=f32= E3=0, the conitutive equations (4) and (5) are reduced in this case to FORMULATION OF LOCAL INTEGRAL EQUATIONS The MLPG method is applied to conruct a weakform over the local fictitious subdomains Xs, which is a small region taken for each node inside the global domain, as shown in [3]. The local subdomains overlap each other, and cover the whole global domain X. The local subdomains could be of any geometrical shape and size. In the present paper, the local subdomains are considered to be of the circular shape. The local weak-form of the governing equation (1) can be written as 11 c11 c 22 12 12 0 c12 c22 0 0 11 0 e21 E 0 22 0 e22 1 E c66 212 e15 0 2 11 E C(x) 22 L(x) 1 E2 212 D1 0 D e 2 21 0 e22 11 e15 h 22 11 0 2 0 12 (9) ij , j (x, ) (x)ui (x, ) 0 E1 h22 E2 * X i (x, ) uik (x) d 0 (11) 11 p1 E G (x) 22 H (x) 1 P(x) p2 E2 212 u ik (x) is a te function. Applying the Gauss divergence theorem to eq. (11) and rearranging, one can obtain (10) ij * * (x, )n j (x)uik (x)d ij (x, )uik , j (x)d s * ik X (x, ) (x)u (x, ) u i i ( x) d 0 (12) c11 c12 0 c12 c22 0 c13 11 11 c23 22 22 0 33 0 uXs is the boundary of the local subdomain which consis of three parts uXs = Ls , C , Csu [3]. Here Ls is the local boundary that is totally inside the global domain, C is the part of the local boundary which coincides with the global traction boundary, i.e., C = uXs + Ct and similarly Csu is the part of the local boundary that coincides with the global displacement boundary, i.e. Csu = uXs + Cu. Heaviside ep function is chosen as the te func* tion u ik (x) in each subdomain at x s * uik (x) ik 0 at x s w* (x) is a te function. Applying the Gauss divergence theorem to the local weak-form (16) and considering the Heaviside ep function for the te function w* (x) one can obtain Ls sa (x)c(x)(x, )d (17) and using this, the local weak-form (12) is converted into the following local boundary-domain equations Ls su t i ( x, ) d (x)u (x, )d ti (x, )d (13) It is important to note that the local integral equations (13) are valid for both homogeneous and nonhomogeneous linear piezoelelectric solids. Nonhomogeneous material properties are included in eq. (13) through the elaic and piezoelectric tensor of the material coefficients in the traction components. The local weak form of the governing equation (2) is given in the similar way as There is no requirement in the MLPG method on the te and the trial functions to be necessarily from the same functional spaces. For internal nodes, the te function is chosen as the Heaviside ep function with its support on the local subdomain. The trial functions, on the other hand, are chosen to be the moving lea-uares (MLS) approximation over a number of nodes spread within the domain of influence. The approximated functions for the mechanical displacements, the electric potential and the temperature can be written according to Atluri [4] as: t t u h (x, x ) = UT (x) $ u ( x ) = / za (x) u a ( x ) ^ h (x, ) a (x) a ( ) n j, j ( x, ) R ( x, ) v * ( x ) d 0 (14) ^ h (x, ) a (x) a ( ) (18) o* (x) is a te function. Again, as in previous case, applying the Gauss divergence theorem to the local weak-form (14) and considering the Heaviside ep function for the te function o* (x), one will obtain Ls sp Q ( x, ) d Q ( x, ) d R ( x, ) d (15) Q(x,x)=Dj (x,x)nj. The local weak-form of the diffusion equation (3) can be written as ij (x), j (x, ) (x)c(x) (x, ) ,i (16) t t t the nodal values ua(x), } a(x) and i a(x) are fictitious parameters for the displacements, the electric potential and the temperature, respectively, and za(x) is the shape function associated with the node a. The number of nodes n used for the approximation is determined by the weight function wa(x). In the MLS approximation a 4th order spline-type weight function is applied [3]. It can be easily shown, that the C1 - continuity is ensured over the entire domain, therefore the continuity conditions of the tractions, the electric charge and the heat flux will be satisfied. The traction vector ti(x,x) at a boundary point xduXS is approximated in terms of the same nodal t values ua(x) as ^ t h (x, ) N(x)C(x) B a (x)u a ( ) n W (x, ) w (x) d 0 20 VOLUME 14, No. 4, 2010 ^ N(x)L(x) A a (x) a ( ) t u / z ( g ) } ( x ) = } ( g , x ) for g ! C a a n (23) (24) ^ N(x) a (x) a ( ) (19) t u / z ( g ) i ( x ) = i ( g , x ) for g !C a a the matrix N(x) is related to the normal vector n(x) on uXS by n 0 n2 N ( x) 1 0 n2 n1 Discretized forms of the unknown quantities in the local boundary-domain integral equations (13), (15) and (17) are obtained in the view of the MLSapproximation (19) - (21) as and the matrices Ba and Aa are represented by the gradients of the shape functions as ^ N(x)C(x)B a (x)d u a ( ) L s su z,a1 0 z,a1 B a (x) = > 0 z,a2 H, A a (x) = < a F z,2 z,a2 z,a1 The electrical charge Q(x,x) can be approximated in similar way by n ^ a (x)d u a ( ) s n ^ N(x)L(x)A a (x)d a ( ) L s su ^ Q h (x, ) N1 (x)G (x) B a (x)u a ( ) -/ d Ls + Csu t # N (x) c (x) z (x) dC n i ( x ) = a a ^ N1 (x)H (x) A a (x) a ( ) t ( x, ) d (25) ^ N1 (x)P(x) a (x) a ( ) (20) the matrices G and H are defined in eq. (10) and N1(x)= [n1 n2]. The heat flux q(x,x) is approximated by ^ q h (x, ) kij ni ,aj (x) a ( ) ^ N1 (x)K (x) A a (x) a ( ) (21) ^ N1 (x)G (x)B a (x)d u a ( ) Ls sp n ^ N1 (x)H (x)A a (x)d a ( ) Ls sp n ^ N1 (x)P(x) a (x)d a ( ) L s sp Q ( x, ) d R ( x, ) d (26) k k K (x) 11 12 k12 k22 Using the approximation formula (18) and fulfilling the essential boundary conditions at those nodal points on the global boundary, the displacements, the electrical potential and the temperature are prescribed, discretized form of the boundary conditions is obtained as ^ N1 (x)K (x)A a (x)d a ( ) Ls sa n ^ c a (x)d a ( ) s (27) t u / z ( g ) u ( x ) = u ( g , x ) for g ! C a a (22) Equations (25), (26) and (27) are considered on the sub-domains adjacent to interior nodes as well as to the boundary nodes on C, C and C. Complete syem of ordinary differential equations for the computation of the nodal unknowns can be obtained by collecting the discretized local boundary-domain integral equations together with the discretized boundary conditions for the displacements, the electrical potential and the temperature. t As nodal unknowns the fictitious parameters u a(x), t t } a(x) and i a(x) are considered. aa General variation of material properties with Cartesian coordinates was assumed in previous formulations. If material parameters are dependent on temperature, one can write for a general material parameter A: A(x) f ( (x)) ^ N1 (x)G ( k 1) (x)B a (x)d u a ( k ) ( ) Ls sp n ^ N1 (x)H ( k 1) (x)A a (x)d a ( k ) ( ) Ls sp n ^ N1 (x)P ( k 1) (x) a (x)d a ( k ) ( ) Ls sp ( x, ) d R ( x, ) d Q (32) (28) ^ N1 (x)K ( k 1) (x)A a (x)d a ( k ) ( ) Ls sa Using equation (28) to replace material parameters in conitutive equations (4) and (5), a nonlinear expression is obtained. This results in necessity of application of iterative approach in each time ep. The linearized conitutive equations in the k-th iteration ep are then given as ( (k ( ijk ) (x, ) cijkl1) ( ) klk ) (x, ) (k ( ekij 1) ( ) Ek( k ) (x, ) ijk 1) ( ) ( k ) (x, ) k ( D (j k ) (x, ) e(jkl1) ( ) klk ) (x, ) k h (jk 1) ( ) Ek( k ) (x, ) p (jk 1) ( ) ( k ) (x, ) n ^ (x) c(x) a (x)d a ( k ) ( ) s (33) (29) (30) Applying now equations (29) and (30) ed above, we obtain a set of ordinary differential equations for the k-th iteration ep in form ^ N(x)C( k 1) (x)B a (x)d u a ( k ) ( ) Ls su n n ^ (x) a (x)d u a ( k ) ( ) s Material parameters are taken at a reference temperature for the 1 iteration ep. The iteration process is ended if the difference of the Sobolev-norms for temperatures in two successive eps is smaller than a prescribed tolerance. Now, we should explain solution of the complete syem of ordinary differential equations. The syem of above given local integral equations can be rearranged in such a way that all known quantities are on the r.h.s. Thus, in matrix form the syem becomes Ax Bx Cx Y (34) In the present work, the Houbolt method is applied. In the Houbolt finite-difference scheme [12], the "acceleration" (ü=) is expressed as ^ N(x)L( k 1) (x)A a (x)d a ( k ) ( ) L s su 2x 5x 4x x 2 2 (35) -/ d # N (x) c Ls + Csu ^ k - 1h t (x) z (x) dC n ia(k) ( x ) = Dx is the time ep. The backward difference method is applied for the approximation of "velocities" t ( x, ) d (31) x x (36) Subituting eqs. (35) and (36) into eq. (34), we get 22 VOLUME 14, No. 4, 2010 the following syem of algebraic equations for the unknowns xx+Dx 1 1 2 2 A B C x 2 (5A B )x 1 A 2 4x x 2 Y (37) The value of the time ep has to be appropriately selected with respect to material parameters (wave velocities) and time dependence of the boundary conditions. te the present computational method. The material coefficients of the panel are considered for PZT5H according to Qin and Mai [23]: c11=12,6.1010Nm-2, c12=5,3.1010Nm-2, c22=11,7.1010Nm-2, c44=3,53.1010Nm-2, e15=17Cm-2, e21=-6,5Cm-2, e22=23,3Cm-2, h11=15,1.10-9C(Vm)-1, h22=13.10-9C(Vm)-1, t=7500kg/m3, k11=50W/Km, k22=75W/Km, a11=0,88.10-51/K, a22=0,5.10-51/K, p1=0, p2=-5,4831.10-6C/Km2, c=420Wskg-1K-1. Plane rain condition is assumed in the analysis. The mechanical displacement, the electrical potential and the thermal field on the finite uare panel with a size a x m x 1m are approximated by using 121 (11x11) equaly-spaced nodes. The circular local sub-domains are considered, each with a radius rloc=0,08. The Sobolev-norm is calculated for the purpose of error analysis. The relative error of the temperature in the considered time interval [0,T] is defined as NUMERICAL EXAMPLES Numerical results are presented in this section for selected problems of thermoelaicity in piezoelectric materials. In order to te the accuracy of the present meshless method a unit uare panel under a sudden heating on the top side is analyzed as the fir example (Fig. 1). The following analytical solution is available for uncoupled thermoelaicity in an isotropic material [6]: ( x2 , ) 1 exp 2n 1 n0 (1) (2n 1) 4a 2 num exact exact (39) (2n 1) x2 cos 2a (38) 2 d and time T is defined through normalized parameter Tk22/tca2=1,3. Fig. 2 presents numerical results for the temperature at the bottom side and the mid-line of the panel. The temperature is normalized by the intensity of the thermal shock i0=0. Obtained results are compared with the analytical results and an excellent agreement is observed. The relative error of the temperature, r, at both lines is less than 0,5%. For the total number of 441 nodes, the relative error r=0,15% has been obtained. The temperature diribution is not influenced by mechanical and electrical fields. Numerical results for the temporal variation of the direct ress v11 are presented in Fig. 3, an excellent agreement can be seen again between the presented and FEM results computed at both considered lines. The FEM results have been obtained by ANSYS code with 3600 quadratic eight-node elements and 1000 time increments. The resses are normalized 1/ 2 a is the side-length of the panel and l=k22/ tc is the diffusivity coefficient. x2 Q=0 t1=t2=0 , =H(-0), Q=0 q= 0 u1=t2=0 Q=0 q=0 22 1 u2=t1=0 ,q=0 , =0 x1 Fig. 1 Uniformly heated piezoelectric uare plate Homogeneous material properties are selected to by the ic value v11 =-1,2959.106Pa. MLPG: x2 /a = 0.5 0.2 0 -0.2 0 0.2 0.4 0.6 0.8 x2/a = 0. analytical x2/a = 0.5 x2/a = 0. k22/ca2 Fig. 2 Time variation of the temperature on two different lines parallel to x1- axis 11/11 MLPG: x2 /a = 0.5 x2/a = 0. FEM: x2/a = 0.5 x2/a = 0. The temporal variation of the electrical potential W is shown in Fig. 4. The potential is normalized by the ic value W=1,4861.104V at the top side of the panel if the pyroelectric material vector QS =(p1,p2) is vanishing. For p2=-5,4831.10-6C/Km2 we have obtained the ic value W=1,455.104V. One can observe that the pyroeletric parameter is only slightly decreasing the electrical potential in the considered problem. In order to control the mechanical quantities (displacements and resses) the inverse piezoelectric effect can be utilized. So as a next load case, we have applied a uniform electrical displacement D0 on the top side of the panel additionally to the thermal load. ionary boundary conditions are considered. In Fig. 5 one can observe the influence of the electrical displacement on the ress component v11. The increasing electrical displacement on the top side of the panel significantly reduces the normal ress. The electrical displacement is expressed through the normalized quantity m = D0e22 / h22c11a11i0. Negligible influence of the pyroelectric parameters on the ress values can be observed. 0 -0.2 p2 =0. p2 = -5.48e-6 -0.2 0 0.2 0.4 0.6 0.8 11/11c110 -0.4 -0.6 -0.8 -1 -1.2 k22/ca Fig. 3 Time variation of the v11 ress -1.4 0.00E+00 5.00E-05 1.00E-04 1.50E-04 Fig. 5 Variation of the v11 ress with the normalized electrical displacements MLPG: x2 /a = 0.5 x2/a = 1. FEM: x2/a = 0.5 x2/a = 1. -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 k22/ca2 Fig. 4 Time variation of the electric potential at two different lines parallel to x1- axis As the next example, an edge crack in a finite piezoelectric rip is analyzed. The geometry of the rip is given in Fig. 6 with the following values: a=0,5m, a/w=0,4 and h/w=1,2. Only a half of the rip is modeled due to the symmetry of the problem with respect to the x1 - axis. We have used 930 nodes equidiantly diributed for the MLS approximation of the physical quantities. On the both lateral sides of the rip thermal load is applied. 24 VOLUME 14, No. 4, 2010 x2 900 t2=t1=0 , Q=0 , q=0 =t Q=0 t1=t2=0 =0 Q=0 t1=t2=0 62 q=0 w =0 u2=t1=0 31 32 1 t1=t2=0 Q=0 a x1 Fig. 6 Edge crack in a finite rip under a thermal shock, half of the rip depicted due to a symmetry For cracks in homogeneous and linear piezoelectric solids the asymptotic behaviour of the field quantities has been given by Sosa [31] and Pak [21]. Using these expressions, one can obtain the following expression for the intensity factors [10] as KI u Re() 1 2 lim r 0 8r K IV the crack line because of the immediate electromechanical interaction. In the dynamic case the ress field is coupled not only to the immediate electric field, but also to inertia forces [9]. Numerical results for the normalized ress intensity factor fI=KI/ ra a11c11i0 are presented in Fig. 7. On the right lateral side of the finite specimen a vanishing heat flux is considered in the transient heat conduction case. We have also considered an induced nonhomogeneity caused by the dependence of the thermal expansion coefficient on the temperature, akl= akl0(1-mi), m=0,005K-1. It can be well observed that for a large time inant the SIF is vanishing in Fig. 7. The electrical displacement intensity factor (EDIF) for considered boundary value problem is vanishing since the thermal processes changes are relatively slow. FEM 1/2 MLPG: homog. induced nonhom. fI=KI/ 11c110(a) (40) the matrix B is determined by the material properties in paper by Garcia-Sanchez et al. [10]. Note that the symmetry conditions of the displacements and the potential with respect to the crack plane are utilized. Contrary to the ic case, the SIF is not vanishing for the thermal shock prescribed on the left lateral side of the rip. From the Maxwell`s equations, it is known that the velocity of electromagnetic waves is equal to the speed of light, which is much greater than the velocity of elaic and thermal waves. Hence, the use of quasi-ic approximation in eqs. (1) to (3) is juified for the interaction of electrical, mechanical and thermal fields. The response of the electric fields is immediate, while that of the elaic and thermal ones are taken as finite because of the finite velocity of both waves. On the other hand, in a ic case, the response of all the mechanical (rain, ress), thermal and electrical fields is immediate. Thus, the SIF is vanishing in such a case since the resses are zero ahead the crack-tip on k22/ca Fig. 7 Time variation of the ress intensity factor for an edge crack CONCLUSION A meshless local Petrov-Galerkin method (MLPG) is proposed for the solution of boundary value problems for coupled thermo-electro-mechanical fields. Transient dynamic governing equations are considered. The material properties of a piezoelectric material are influenced by a thermal field. It is leading to an induced nonhomogeneity and so to mathematical complications of the governing equations compared to a homogeneous counterpart. The analyzed domain is divided into small overlapping circular subdomains. A unit ep function is used as the te functions in the local weak-form of the governing partial differential equations. For the approxi- mation of the physical field quantities the Moving Lea-uares (MLS) scheme is used. The proposed method is a truly meshless method, which requires neither domain elements nor background cells in either the interpolation or the integration. The present method represents an alternative numerical tool for many exiing computational methods like the FEM or the BEM. The main advantage of the present method is its simplicity. The present formulation is promising for numerical analyses of multi-field problems like piezoelectric or thermoelaic problems. ACKNOWLEDGEMENT The authors gratefully acknowledge the support by the Slovak Science and Technology Assiance Agency regiered under number APVV-0427-07 and the Slovak Grant Agency VEGA-2/0039/09.

Journal

Acta Mechanica Slovacade Gruyter

Published: Dec 1, 2010

There are no references for this article.