3D FEM analysis of stress concentrations in a rectangular thick plate with PZT inclusions under bend

Abstract

Keywords: Piezoelectric Inclusions; Stress Concentration; 3D Finite Elements Method

Introduction

The piezoelectric property is a property seen in some materials without central symmetry in their crystalline internal structure. When an external force is applied to these materials, that is, when the material changes shape, polarization occurs in its structure. In other words, there is a situation that can create an electric current in the structure

(direct effect). Conversely, when the material is placed in an electric field, mechanical deformation occurs in its structure (indirect effect). These effects are related to the change in polarization density inside the material. The response of piezoelectric materials to mechanical action depends on the crystal orientation determined by the arrangement of atoms in the crystal structures and the angular values ​​between the

*Corresponding author. *E-mail address: ubabuscu@yildiz.edu.tr This paper was recommended for publication in revised form by Regional Editor Ahmet Selim Dalkilic Published by Yıldız Technical University Press, İstanbul, Turkey Copyright 2021, Yıldız Technical University. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

applied mechanical stress and the direction of the elongation [1]. The dual properties of piezoelectric materials are utilized in many advanced technological applications such as actuators, sensors, transducers and generators. Recently, piezoelectric materials, which have a wide research area, have generally been used by combining (integrating) them with non-piezoelectric materials. There are many studies in the literature on piezoelectric materials containing discontinuities in various geometric forms. Wang [2] develop a solution for an infinite, piezoelectric medium containing a piezoelectric, ellipsoidal inclusion by using Green’s function technique. A piezoelectric ellipsoidal inhomogeneity embedded in a non-piezoelectric elastic matrix was analyzed via equivalent inclusion method by Fan and Qin [3]. Xu and Rajapakse [4] considered coupled elastic and electric fields in piezoelectric solids with arbitrary shaped defects (cavities, inclusions, cracks, etc.) under the assumptions of plane strain or stress conditions by a boundary integral equation based on closed form Green’s functions. Effect of cracks on structures were studied by using finite element method by Gonenli, Das and Ozturk [5,6,7]. Xiao and Bai [8] investigated stress analysis for a circular piezoelectric fiber sensor embedded in a non-piezoelectric elastic material based on Eshelby’s equivalent inclusion method. Spherically symmetric deformation of an inclusionmatrix problem, which consists of an infinite isotropic matrix and a spherically uniform anisotropic piezoelectric inclusion was studied in [9]. Mishra et al. [10] analyzed an elliptical piezoelectric inclusion embedded in an infinite piezoelectric matrix in the framework of linear piezoelectricity by using the conformal mapping technique in [10] and in the framework of antiplane piezoelectricity in [11]. The piezoelectricity problems containing arbitrary elliptical inhomogeneities in an infinite medium was studied by Lee et al. [12]. The Green’s function for an infinite piezoelectric medium with an elliptic hole filled with air or without air

under the generalized line dislocation and the generalized line force was given by Zhou et al. [13]. They discussed the interaction of an arbitrary distributed small crack and an elliptic hole in plane piezoelectric medium subjected to uniform electro-mechanical loads at infinity. Yang and Gao [14] addressed the plane problems of multiple piezoelectric circular inclusions in a non-piezoelectric matrix based on the complex potentials. The problem of the stress state in an orthotropic piezoelectric body with a triaxial ellipsoidal inclusion under homogeneous force and electric loads was studied by Kirilyuk and Levchuk [15]. The problem was solved by the Eshelby method of equivalent inclusion generalized to the case of a piezoelectric orthotropic space. In [16], the electrostatic analysis of a thick plate containing a prismatic inclusion made of piezoelectric material under the effect of axial tensile force was investigated. In summary, the aforementioned works focused on problems involving piezoelectric inclusion(s) and sensors in various forms for infinite media or plane problems. In this study, a mathematical model for the three-dimensional stress analysis of rectangular plates under the effect of bending and containing rectangular prismatic inclusion(s) made of piezoelectric material, is presented with the help of the piecewise homogeneous body model and the 3D exact equations of the electro-elasticity theory, which are numerically examined using the three-dimensional finite element method. Using the mathematical model presented, the effect of the size, location of the piezoelectric inclusion/ inclusions and volume fraction as well as the coupling effect between the electrical and mechanical fields on the static analysis of the rectangular plate are studied. Theory It is assumed that the rectangular prismatic inclusion(s) in the structure are located along the posterior thickness of the plate (Ox3 axis) and are parallel to the Ox3 coordinate axis (Figure 1).

(a) (b) Figure 1. The geometry and loading of (a) plate containing a single PZT inclusion (Problem 1) b) plate containing double PZT inclusions (Problem 2).

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

The solution region of the problems discussed are: For Problem 1: (1)

component, is the displacement vector’s component, is the electrical displacement vector’s component, is the electric is the electrical field vector’s component, is the piezoelectric constant, is the dielecpotential, is the component of tric (permittivity) constant, and elastic stiffness matrix. Compact matrix notation is used to write constitutive equations in matrix form. This notation reduces the pairs of indices ij or kl to a single index such as p or q as follows:

Here, for Problem 1 (Problem 2), ( and ) shows the region where the inclusion(s) are located, and or shows the region outside the inclusion(s) (matrix). The plate is simply supported on all its lateral surfaces, that the electrical potential is zero, and that it is subject to a uniformly distributed pressure force only from its upper surface. In addition, it is assumed that no mechanical and/ or electrical loads act on the lateral surfaces of the plate or the inclusion and matrix interfaces. It is accepted that ideal contact conditions are provided at the interface between the inclusion contained in the plate and the matrix, and that the PZT inclusion is grounded from all surfaces (ie, there is no electrical displacement change in the normal direction). It should be noted that the mathematical modeling will be given by assuming that the material of both the inclusion and the surrounding region (matrix) is a piezoelectric material in order to be more general. However, in numerical calculations, the case where the matrix material is an elastic material will be taken into account. In mathematical modeling, the superscript (1) will be used to represent the dimensions of the matrix material and the superscript (2) and (3) will be used for the inclusion material. Accordingly, the equilibrium equations provided in the solution region are:

If the constitutive equations of the material are arranged according to the new index notation: (5) or in matrix form:

the following boundary conditions are satisfied on all the lateral surfaces:

In Eqs. (3) and (4), when there is a single (double) inclusion in the structure, the parameter m takes the values m=1,2 (m=1,2,3) and all quantities referring to the matrix (inclusion) material will be marked by the superis script 1 (2 and 3) and the following notation is used: is the Green strain tensor the stress tensor component,

and the following contact conditions at the interfaces are satisfied:

In Eq. (8), k=2 (k=2,3) for Problem 1 (Problem 2). (k=1,2,3) is the component of the unit outward normal

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

vector of surface S, and S represents the surface of the inclusions and is equal to Eq. (9) for Problem 1, and Eq. (10) for Problem 2.

The boundary value problem is solved numerically with the help of three-dimensional (3D) finite element modeling. The solution of the problems discussed was made for the half (1/2) region by making use of the symmetry of the structural element and the loading with respect to . The solution region is discretized into 8-node rectangular prismatic finite elements separately for the matrix and the inclusion material in accordance with the geometry (Figure 2). The electrical potential (φ) and the displacements in the direction of the three axes are taken as unknowns in the nodes. Shape functions, defined by the product of linear Lagrangian interpolations, are used at the nodes of the rectangular prismatic sample element (Figure 2) [17]. Displacement and electrical potential functions in each finite element within the framework of the finite element method

Figure 2. Rectangular prismatic sample element, position of nodes and coordinates of nodes linear algebraic equation system. In Eq. (13), K is the coefficient (Stiffness) matrix, u is the vector containing the unknowns, and F is the right-hand side (force) vector. The reduced system is found by applying the displacements to the system of Eqs. (13) and the boundary conditions according to the electric potential given in Eq. (7). With the solution of the reduced system, the sought quantities in the nodes are obtained and the stress distributions in the structure are determined with the help of the constitutive equations given in Eq. (6). It will be attempted to determine the mutual effects of various geometric and material parameters and electrical and mechanical fields on these stress distributions.

Results And Discussion

(11) are represented in serial form. For FEM modeling of the boundary value problem, the functional that expresses the total electro-mechanical energy accumulated in the structure

(12) is used [18]. The solution of the boundary value problem handled with the help of the functional, Eq. (11) expressions, and the known Ritz technique is reduced to the solution of (13)

In solving the boundary value problems with the help of the finite element method, it is assumed that the inclusion material consists of piezoelectric (PZT) material, and the matrix material surrounding the inclusion consists of elastic material. Since the considered problem is symmetric plane, only half of the region is used. about the For the FEM modeling, 6400 rectangular prismatic finite elements, 7497 nodes and 27591 NDOFs are used in total. First, for testing the mesh sensitivity, Figure 3 is given. The parameter N, M and NZ shows the number of rectangular FEs along the Ox1 axis, Ox2 axis and Ox3 axis, respectively. Figs. 3(a-c) show the influence of (a) N, (b) M and for the case and (c) NZ on the values of in the cross section of . It approach follows from the graphs that the values of a certain limit value by increasing the parameters. These

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

Figure 3. Influence of the number of rectangular FEs along the (a) Ox3 axis i.e. NZ, (b) Ox1 axis i.e. N and (c) Ox2 axis i.e. for the case and in the cross section of . M on the values of

Table 2. The values of modulus of elasticity, density and Poisson ratio of the selected metals Metal

results confirm the reliability of the mesh sensitivity used for determining the numerical solution. Thus, for the numerical results, 6400 rectangular FEs are used in total by taking the parameters NZ=20, N=20 and M=16. A FORTRAN code is implemented for all the programs and algorithms required for the solutions of the boundary value problems using the finite element method. The definite integrals can be evaluated numerically by Gaussian quadrature with 10 sample points. The material constants used in the numerical calculations for the piezoelectric material constituting the inclusion material are given in Table 1. The piezoelectric inclusion(s) is/are embedded in non-piezoelectric metal materials. The mechanical constants of some of the metal materials used are given in Table 2. It should be noted that in all numerical calculations, the polarization direction for the PZT material is taken as Ox2. In this case, the material matrices given in Eq. (6) for the PZT material for the Ox2 poled direction [18] are as follows:

For determination of the effect of the electro-mechanical coupling on stress distributions, numerical results are determined in two cases. Case 1: The dielectric and piezoelectric constants of the . PZT plate’s material are equated to zero i.e., Case 2: The piezoelectric constants of the PZT plate . materials differ from zero i.e., In Figure 4, the stress distribution in the rectangular thick plate made of homogeneous isotropic material at different parameter values is . As the examined in the cross section of parameter increases, the length of value of the the plate in the direction of the specified axis increases. Accordingly, as a result of this parameter change, it will be possible to approach the limit values that can be obtained in the plane strain condition of this problem for values. From the graphs it follows that the increase with and converge to values of the value ([19]) that should be obtained in the case of plane strain. This gives us confidence of the programs and algorithms we created. In Figure 5, the influence of the volume ratios of the single PZT inclusion in the plate (Vd / V) on the values of

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

Figure 4. The effect of on the values of for a homogeneous-isotropic plate under bending in the cross section of .

Figure 5. The effect of Vd / Von the values of the cross section of .

tigated in the cross section of . Note that V represents the total volume of the plate and Vd represents the volume of the PZT inclusion. As can be seen from the graph, as the value of the parameter Vd / V decreases, that is, as the volume fraction of the PZT inclusion in the plate decreases, the absolute values of the electrical potential increase in the considered cross section.

In Figure 6, the influence of the volume ratios (Vd / V) of the single PZT/Al inclusion for Al matrix material on the values of the differences of the stress distribution and (b) of (a) along the Ox1 axis is investigated in the cross section of for both cases. As seen in Figure 6 (a and b), as the volume of the inclusion increases, the and (b) absolute values of (a) increase. In addition, as the volume of the inclusion increases, the difference between the stress

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

(a) Figure 7. The effect of Vd / V on the values of (a) . values obtained in Case 2 and the values obtained in Case 1 increase for . For the values of the differences of stress distribution of , inequality is obtained for Case 2, while inequality is obtained for Case 1. In Figure 7, the influence of the volume ratios (Vd / V) of the single PZT/Al inclusion for Al matrix material on the values of the differences of stress distribution of (a) and (b) along the Ox1 axis is investigated in the cross section of for both cases. In both graphs, the absolute values of stress differences obtained for Case 1 remain larger than the stress differences obtained for Case 2 for the cross section without edge effects. In other words, the fact that the inclusion material is PZT reduces the stress difference values. As the inclusion volume decreases, the difference values obtained for Case 1 and Case 2 increase. The presence of PZT inclusion in the plate causes an average of 40% change in stress values (Figure 6-Figure 7) at the interface between the matrix and PZT inclusion(s). Figure 8 shows the influence of the upward replacement of the position of the PZT inclusion (i.e., reducing of the parameter ) on the values of along the Ox2 axis for PZT-5H inclusions with Ox2 polled axis/ AL matrix for both cases in the cross section of , . From the graphs it follows that the absolute values of the electrical potential decrease with decreasing along the Ox2 axis. Figure 9 shows the effect of the upward replacement of the position of the PZT inclusion (i.e., reducing of the parameter ) on the values of (a) and (b) along the Ox1 axis for PZT-5H inclusions with Ox2

polled axis/ AL matrix for both cases in the cross section of . From these graphs it follows that the values of the stresses of decrease with decreasing for both cases, but the difference between the cases increases (Figure 9a). On the other hand, as the value of

Figure 8. The effect of on the values of for PZT-5H inclusion with Ox2 polled axis/ AL matrix for Case 2 in the cross section of , .

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

(a) on the values of (a) Figure 9. The effect of both cases in the cross section of

(a) on the values of (a) Figure 10. The effect of for both cases in the cross section of

the parameter decreases, the absolute values of the stresses of within the boundaries of inclusion increase for Case 2, while they decrease for Case 1, and therefore the difference between these cases decreases (Figure 9b). Figure 10 shows the effect of the upward replacement of the position of the PZT inclusion (i.e., reducing of the parameter ) on the values of (a) and (b) along the Ox1 axis for PZT-5H inclusions with Ox2 polled axis/ AL matrix for both cases in the cross section

of . From the graphs it follows that the absolute values of the shear stresses of increase with decreasing within the boundaries of inclusion, and the change in Case 2 is more than the change in Case 1 (Figure 10a). On the other hand, as the value of the parameter decreases, the values of the stresses of within the boundaries of inclusion decrease for both cases, but the difference between these cases increases (Figure 10b).

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

eter on the values of the electrical potential for PZT-5H inclusions with Ox2 polled axis/ AL matrix for both cases in the cross section of . It

should be noted that as the value of the ratio decreases, the PZT inclusions get closer to each other, and =0 indicates the case where the plate contains a single inclusion (Problem 1). The total volume of the two inclusions contained in the plate in Problem 2 is equal to the volume of the single inclusion considered in Problem 1. As can be seen from the graph, as the value of the parameter decreases, that is, as the PZT inclusions get closer to each other, the electric potential values in the considered section are affected very little. Figure 13 shows the effect of the distance between the inclusions along the Ox1 axis, i.e. the influence of the parameter on the values of the stresses of (a) , (b) for PZT-5H inclusions with Ox2 polled axis/ AL matrix for

Figure 11 indicates the values of the stresses of (a) , (b) , (c) and (d) along the Ox1 axis for the various materials of the matrix with PZT-5H inclusion for the Ox2 polled axis. It can be concluded that the difference between Case 1 and Case 2 is greatest for the ST matrix and smallest for the MG matrix for all stress distributions. So, the effect of piezoelectricity on the stress distributions is maximum for the ST matrix material. Figure 12 shows the effect of the distance between the inclusions along the Ox1 axis, i.e. the influence of the param-

Figure 11. The effect of various matrix materials on the values of (a) inclusion with Ox2 polled axis for both cases in the cross section of

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

Conclusion

Figure 12. The effect of on the values of for PZT-5H inclusions with Ox2 polled axis/ AL matrix for . both cases in the cross section of

(a) on the values of (a) Figure 13. The effect of for both cases in the cross section of

both cases in the cross section of . When there are two parallel inclusions in the structure, the greatest change in the stress distribution is in the starting and ending sections of the inclusions. As can be seen from the graphs given in Figure 13, as the value of the parameter decreases, that is, as the PZT inclusions get closer to each other, the places where the stress distribution changes are displaced according to the inclusions’ positions, and, as the inclusions get closer to each other, the numerical values are affected very little for both stresses. The mutual effect of inclusions does not make a significant difference for elastic/ PZT inclusions.

In this study, the electro-static analysis of the rectangular prismatic form of the composite thick plate containing PZT inclusion(s) under the effect of bending is solved numerically with the help of the Finite Element Method, within the framework of the three-dimensional exact equations of the theory of electro-elasticity. The effects of various matrix materials, the size, volume fraction and location of the piezoelectric inclusions as well as the coupling effect between the mechanical and electrical fields on the stress distribution in the structure are investigated. Based upon the discussion in this paper, the following conclusions can be established: • The absolute values of electrical potential decrease with the volume fraction of the PZT inclusion in the plate; • PZT inclusion(s) within the rectangular plate under bending force causes a decrease in the values of normal stresses and causes an increase in the values of shear

stresses at the interface between the matrix and PZT inclusion(s). The fact that the inclusion material is PZT increases the and ; and reduces stress difference values of and ; the stress difference values of As the volume of the inclusion in the structural element increases, the difference between Case 2 (considering the piezoelectric characteristics of the PZT material of inclusion) and Case 1 (considering the material and , but inclusion as elastic) increases for and ; decreases for

Sigma J Eng Nat Sci, Vol. 41, No. 6, pp. 1197−1208, December, 2023

As the PZT inclusion approaches the upper face plane of the plate for both cases, the values of the stresses of and decrease for both cases, the absolute values of the stress of decrease for Case 1 but increase for Case 2, and the values of the stresses increase; • The absolute values of the electrical potential decrease as the PZT inclusion approaches the upper face plane of the plate along the Ox2 axis; • The effect of piezoelectricity on the stress distributions is maximum for the ST matrix material; • When there are two parallel inclusions in the structure, the greatest change in stress distribution is found in the starting and ending cross sections of the inclusions. As the PZT inclusions get closer to each other, the places where there is a change in the stress distribution are displaced according to the position of the inclusions, and, as the inclusions get closer to each other, the numerical values ​​are slightly affected for the stresses; • The mutual effect of inclusions does not make a significant difference for elastic/PZT inclusions. We recall that the aim of present investigation is to provide new scientific information about the effect of a PZT inclusion(s) on the static characteristics of elastic plates. Since the exact theory is used, it is valid for both thin and thick plates and forms a reference to studies with approximate plate theory. Different geometric form of PZT inclusion(s) may be used in some future applications to decrease the stress concentrations of construction elements. Moreover, in the present study, the mathematical modeling of the problem under consideration is given for the geometrically linear case. These investigations can be extended to the study of nonlinear problems. However, under the solution of these nonlinear problems, the finite element mesh and some parts of the FEM modeling used in the present study can also be employed.

Share and Cite