Abstract
The prediction of the overall effective properties of fibre-reinforced piezocomposites has drawn much interest from investigators recently. In this work, an algorithm used in two-dimensional (2D) analysis for calculating transversely isotropic material properties is developed. Since the finite element (FE) meshing patterns on the opposite areas are the same, constraint equations can be applied directly to generate appropriate load. The numerical results derived using this model have found a good agreement with those in the literature. The 2D algorithm is then modified and improved in such a way that it is valid for three-dimensional (3D) analysis in the case of random distributed shorts and inclusions. Linear interpolation of displacement field is employed to establish constraint equations of nodal displacements between two adjacent elements.
1. Introduction
Piezoelectric materials have unique features of converting mechanical field into electric field or vice versa. Being the most widely used “smart” material, it has drawn an increasing attention in the past several decades and much effort has been devoted to the research of the fabrications and applications of piezoelectric materials [1, 2]. These materials have been widely applied in the industry for fabricating devices including sensors, transducers, and actuators in some of the most promising areas such as aerospace industry and biotech. To enhance understanding of electromechanical coupling mechanism in piezoelectric materials and to explore their potential applications in practical engineering, numerous investigations, either analytically or numerically, have been reported over the past decades [3–11]. Due to the limitation of analytical methods which are available for problems with simple geometries and loading conditions numerical methods are usually employed in solving most practical problems [5, 12].
Since piezoelectric composites are heterogeneous on a microscopic level, the definition and study of the material properties are hard to impose on a specific area. Historically, homogenization approach [13, 14] has been applied to investigate the effective properties of such heterogeneous materials. And micromechanics approach [14] seems preferable for revealing the relations between effective properties and microstructure. The idea of homogenization method herein is to model an averaging statistical homogeneity through representative volume element (RVE) approach. It has been employed to investigate effects of different types of inclusion, such as microdefects or fibre reinforcement on the properties of the materials [15–18]. Among the investigators of effective properties of piezoelectric composites, Pettermann and Suresh [19], Berger et al. [16], and Kar-Gupta and Venkatesh [20] have investigated the overall properties of 1–3 piezoelectric materials, which contain piezoelectric rods embedded in a polymer matrix and aligned through the thickness of the device. In addition, the effect of volume fraction over the effective properties has been covered in their research. Kar-Gupta and Venkatesh [20] systematically assessed the effects of fibre distribution on the piezoelectric coefficients by comparing 14 different fibre arrangements. Recently, FE and boundary element (BE) methods [21–26] have been extensively used in homogenization process, for one of the reasons that both FE and BE methods are applicable regardless of the type of inclusion in the unit cell. However, the adoption of numerical methods raises two major issues: one is that results generated from this method are highly meshing dependent The correctness and accuracy can be affected by the modelling and meshing methods selected; another is that the RVE presents continuity in the homogeneous media, which means the deformation of such an element should keep continuous within an element and on the interfacial boundary of any two adjacent elements as well. In such context, the boundary conditions to be applied are critical to the success of modelling. Finally, it is worth mentioning that the enhancement of commercial FE software packages such as ANSYS provides great opportunities to conduct such analysis. With the aid of such packages, verifications of previous analytical results can be achieved; more importantly, given the ease of manipulating the inputs, such as the shape and configuration of the virtual models, the results of numerous analyses will guide the realistic experiments to be more result oriented and efficient.
Besides, in using RVE to evaluate materials effective properties, the condition of the periodicity of a unit cell should be considered. To this end, the periodic boundary condition is to be applied to impose a uniform deformation field and subsequently the same stress field. In this case, no gaps or overlaps should exist between adjacent cells. The general periodic boundary condition given by Havner [27] provided a strict boundary condition and has been regarded as one of the fundamentals for introducing the concept of RVE. It is noted that, based on the study of interacting periodic cells, Aboudi [13] in his work has suggested a different homogeneous boundary condition, where the displacement field is linearly distributed at the boundaries. Xia et al. [28] then indicated that this boundary condition is only valid for normal tractions, whereas the unit cell is not only overconstrained, but also against stress/strain periodicity if shear loading is applied. For the simplicity of application, a modified periodic boundary condition has been derived and used in the work afterwards [16, 24, 28].
It should be mentioned that the cost of the fabrication of laminated fibre composite is considerably high because of precise requirements and restrictions to the reaction conditions. On the other hand, random distributed short fibre or particle inclusions provide approximate advantages as fibre laminates such as high specific strength and light weight while the cost is much lower for mass production [25, 26]. For such reasons, the prediction of properties of piezoelectric composites with random distributed short fibres becomes important. Practically, the prevalent literature is available presenting the approximation of effective properties of piezoelectric composites with aligned fibre arrangement using FE method based homogenization theory on RVE. The knowledge and the modeling for the piezoelectric composites with random distributed short fibre are, however, limited.
In this study, we aim to develop a practical algorithm for the analysis of 3D piezocomposites with random inclusions. A 2D FE model is first evaluated and verified as an introduction of the modelling and processing actions adopted by ANSYS Parametric Design Language (APDL). For 3D modeling, periodic boundary conditions for a free-meshed unit cell is defined so that it is suitable for any composite unit cell with random configuration or inclusion such as particle, fibre reinforced, or microdefects which can also be generated using different algorithms.
2. Constitutive Equation, Periodic Condition, and Meshing
2.1. Effective Constitutive Relations
For the transversely isotropic composite studied here (see Figure 1), the effective constitutive relation of linear piezoelectricity, which is extensively used in the characterization of piezocomposites in this study, is defined as where is the stress tensor, the electric displacement, the electric field, the fourth rank elasticity tensor, the second order dielectric tensor, the piezoelectric constant, and a bar over a variable here stands for volume average.
(a)
(b)
The prediction of effective coefficients appeared in (1) requires the adoption of periodic boundary conditions to generate appropriate loading.
2.2. Periodic Boundary Condition
As the homogeneous medium consists of periodic unit cells, periodic boundary conditions are required to apply on the boundaries of the RVE. The general periodic conditions expressed by Havner [27] can be applied to ensure periodic displacement and subsequent stress field: where and denote the displacement and stress in three dimensions, respectively, represents any point in the periodic domain, and represents the periodicity. Applying this displacement condition to the boundary of the unit cell in Figure 1 yields which means the three-dimensional displacement vector for any pair of corresponding locations on areas , , and should be the same. A more explicit periodic boundary condition is then given as [14] where the average strain is included as an arbitrarily imposed constant strain and denotes the periodic part of displacement component, which depends on the global loadings.
Based on the boundary condition (4), Xia et al. [28] and Berger et al. [16] developed a unified periodic boundary condition:
In (5) above, for the constant terms on the right side of the equation, , , and represent the normal loads which are either traction or compression, while, , , and represent the in-plane shear load. Equation (5) is extensively used in the calculation in the later sections.
2.3. 2D FE Modelling
Figure 1 illustrates the configuration of the homogeneity of continuous fibre reinforced 1–3 composite by Berger et al. [16]. In this case, the cylindrical fibres are in square arrangement and poled along the direction. This RVE configuration will be the focus in this section.
2.3.1. Element Type and Material Property
SOLID226 in ANSYS element library is used, which is a 20-node hexagonal shaped element type with 3D displacement degree of freedom (DoF) and additional voltage degree of freedom. This element type is easy for the implementation of periodic boundary conditions. And in the later development of 3D model, this element type will be suitable as the meshing method used for irregular model is only valid with tetrahedral element.
The material properties inputs are based on Berger et al. [16], which is listed in Table 1 for the readers’ convenience.
2.3.2. Element Mesh
For meshing, the area geometry is generated first and then sweep mesh is used to further generate the volume. In this way, the meshing result on is the same. In addition, with the setting of the RVE edge line divisions, meshing results on and are also the same. Ultimately, this provides explicit convenience to imposing periodic boundary conditions. As illustrated, when dealing with the situation when volume fraction is specified, say, 0.666, the outline of the fibre circle is much closer to the RVE edge; in this case, a lower density of element as shown in Figure 2(a) is not sufficient for accurate analysis since the elements between the boundary of the RVE and the fibre have been lessened and shown distortion. When the edge division is 40 as indicated in Figure 2(b), the meshing quality is significantly improved.
As for periodic boundary condition, specific boundary conditions will be assigned to the exact opposite positions, namely, , , and −. For example, in -direction where denotes the displacement component in three dimensions and the subscripts and are the node number of any pair of nodes on opposite locations, and area, respectively.
The boundary conditions are shown in Figure 3.
Since the meshing scheme has ensured that there exists a pair of corresponding nodes at the opposite positions, the problem lies in developing a method to apply constraint equations on each pair of node for overcoming the problem of time-consuming over the node pair selection by graphical users’ interface. An internal programme has thus been designed for accomplishing the task. The procedures of the implementation are described as follows.(a)Define the areas , , and as master areas, while , , and as slave areas. Establish two arrays containing the node number and coordinates (, ) of each node (the coordinate is not necessary, because the nodes are located on and areas where the coordinate is a constant).(b)Start from the first node in master array; get the node number .(c)Use the coordinates (, ) of the node to find the node at the exact opposite location, ; ; select the node from the slave array.(d)Given the node number of the nodes on opposite location, constraint equations could be established.
The same procedures are adopted on and areas, whilst the coordinates obtained and stored will be and , respectively.
When integrating the constraint equations in three directions, special care has been taken to avoid overconstraint over the edges that connect areas , , and . Overconstraint may occur when the degree-of-freedom of one node is specified more than once. For example, when applying - in-plane shear load via constraint equations as shown in Figure 4, based on the periodic boundary conditions for , the DoF relations between nodes 1 and 2 are and while as to areas , there will be relations between nodes 2 and 3 that and ; the same problem will also occur in node 3. In this case, when applying constraint equations over , the corner nodes of the RVE will be excluded to avoid overconstraint.
3. Numerical Results of Effective Coefficients
Proper boundary conditions with strain load are specified pertaining to the calculation of different coefficients. For example, for the calculation of and , the boundary conditions are applied in such a way that, except for the average normal strain in direction , all the other mechanical and electric strain components are set to be zero. By this means, in (1), stiffness tensor and can be derived by Practically, this is achieved by setting the -displacement on , -displacement on areas, and -displacement on areas to be zero; electric field on all areas to be zero. If we adopt periodic boundary conditions discussed before to , where , the periodic boundary condition in (5) is simplified to which indicates that all nodes on area will present a displacement in direction. The calculation of other coefficients follows the same routine. Table 2 summarizes the applied loads and boundary conditions for each coefficient to be derived.
The numerical results obtained are shown in Figure 5 and compared with the FE results from Berger et al. [16]. The blue curve shows the data results from this work and the pink curve as the results in Berger et al. [16]. The curves shown in Figure 5 indicated a good match pertaining to elastic tensors, except for when the fibre volume fraction exceeds 0.444. Since the calculation of is highly dependent on the out-of-plane shear strain , and the implementation of such strain requires the original form of periodic boundary conditions rather than a normal displacement applied on one area when calculating , , , and , the accuracy of is significantly dependent on the meshing density than others.
The other observation is that, as for piezoelectric and permittivity tensors, the trend matches with that presented in the literature well.
4. RVE Modification for 3D FE Analysis
The 2D model developed in Section 3 cannot be adopted in 3D problems with anisotropic material properties because the way of applying constraint equations on a 2D model is not allowed on a random meshing scheme. In this case, a 3D model needs to be developed (Figure 6).
4.1. Modifications to the Previous Unit Cell
The most noticeable variation to the 2D model is the meshing method adopted. As in the model to be used for 3D analysis, sweep meshing method is not allowed for the irregular geometry. As a result, free mesh will be used which meshes the volume depending on the geometry of the inclusion inside the volume and the geometry patterns on the surfaces. In the free mesh method, SOLID227 is selected as the subsequent element type, which is a 10-node tetrahedral shaped element with 3D displacement DoFs and additional electric DoF as used in SOLID226 for the 2D model. The material properties of PZT and polymer are the same as that in 2D model. The geometry of the RVE is designed to be the same as the one used for applying out-of-plane loading. As discussed in Section 3, due to the restriction to the computation capacity, an edge division of 12–15 is adopted in this model to generate the most accurate results without offending the computation limitation. The geometry and meshing results are shown in Figure 7.
4.2. Implementation of Periodic Boundary Condition
4.2.1. Methodology
The most principal part of 3D model generation is the application of periodic boundary conditions in order to maintain periodicity on exact opposite location of an RVE. The nature of the 2D model offers the privilege that, at each pair of opposite locations, there exist a pair of nodes generated from sweep meshing, whereas, for 3D model, the mesh patterns on opposite areas are different. There are three scenarios: when projecting the node on one area to the other, there is one existing node on the opposite location as node 1 (1′); the projected node falls onto the edge of an element as node 2 (2′); the projected node falls inside an element as node 3 (3′), as shown in Figure 8. In order to apply periodic boundary conditions on opposite locations, for scenario 1, constraint equations are established directly as in 2D model, while, for scenarios 2 and 3, constraint equations are required to be established between the existing node and a hypothetical projected node.
The method of linking the existing node on one area to the hypothetical projected node on the opposite area is developed and illustrated in Figure 9. Validation is made to show the feasibility of the method to link the existing node on one area to the element on the opposite area that the projected node falls in. Using the derived shape functions of the selected element, a linear element displacement field is constructed. Consequently, the displacement of the projected hypothetic node is obtained in terms of the displacement of the corner nodes of the element.
Furthermore, since the projected node is always on the surface of the RVE, it means that the displacement result of the projected node is regardless of the fourth corner node inside the RVE as shown in Figure 10. The displacement component of the red node in -direction is interpolated in the same manner as in , -direction. Finally, the problem is simplified to the plane stress triangular element with 3D displacement DoFs.
4.2.2. Element Displacement Function Domain
After free meshing of the RVE, elements and corresponding nodes are generated subsequently. The same presteps described in Section 3 are taken to obtain the number and location of each node on all 6 areas of the RVE. The node numbers and locations are restored in separate arrays as original data. The major issue is that: given one certain node on the slave area , on the opposite area three nodes need to be selected to establish element displacement field.
As illustrated in Figure 11, theoretically as long as the projected node (black) falls inside the triangular element, the selected nodes (black) form and shape functions can be derived directly based on the coordinates of the three corner nodes. However, this will inevitably reduce the accuracy of the calculation because the triangle element with a greater area is used for linear interpolation. Therefore, the white element is the most suitable one for linear interpolation as it is the smallest triangle element that the path falls in, and the three corner nodes can be efficiently selected due to the fact that they are already the corner nodes of one element from meshing.
A path is applied to select the element that the projected node falls into. The path is traced by selecting nodes or specifying locations along it. Using the “Path” command, ANSYS will select the elements that the path runs through as illustrated in Figure 12. The path is defined by selecting the element with one of its nodes at the location of the projected node.
4.2.3. Tolerance of Element Selection
Some of the testing results have manifested that, in ANSYS when dealing with entity selection, such as element selection in this case, tolerance is one of the most significant factors affecting the results and needs to be considered. To be specific, if the path node falls slightly away from the edge of elements as illustrated in Figure 13, without properly defined tolerance for element selection, there is much possibility that the wrong element is selected. (It is noted that more than one elements share the edge line, such as Elements I and II and certain elements inside the RVE.)
Due to the fact that the description for manipulating the tolerance when using path to select element is vague, at this stage there has not been a way developed to ensure the exact element that the path node falls in is selected; as a result, problematic cases such as those illustrated in Figure 14 may occur.
In the first graph of Figure 14, nodes on area (top) are projected to the area (bottom), when dealing with node number 2882 on , the coordinates of the projected node or the path node are , , and ; when selecting the element using the path node, the blue element is selected while the corner node close to the path follows the coordinates , , and . The tolerance requirement reaches . In such cases, linear interpolation for displacements is conducted in a problematic element.
4.2.4. Linear Interpolation for Displacement
When an element is selected to be the domain for linear interpolation, corner nodes or corner nodes on the RVE face are further selected due to the reason that the displacement components of the node on the boundary of the RVE could be determined without considering the fourth node inside the RVE. The three-node triangle element which is called Turner triangle is shown in Figure 15(a). The area of the triangle can be calculated using the determinate,
It should be noted that the area given by (9) is a signed quantity. The corner nodes should be numbered in a counter-clockwise manner so that the area obtained is positive. Now it is assumed that the displacement varies linearly within the triangle domain; conventionally, the linear interpolation functions are selected as where and are the coordinates of an arbitrary point inside the triangle, and , , and are coefficients to be determined from the corner nodes. Since in this case displacement in -direction is treated with the same manner, the function is expanded as By transformation into matrix form, the relations can be expressed as where , , and are shape functions in terms of the corner nodes, which are given by
The coefficients , , and can be determined using (14) below; with automatically changing the corresponding subscript by adding 1 for each time, all the coefficients could then be determined:
By substituting (13) and (14) into (12), the displacement components of any point within the triangle element can be determined, given the coordinates of the point, the coordinates, and the displacement components at the three corner nodes.
The simplified linear interpolation relations are given when the unknown point lands on the edge of the triangle element. In this case, the displacement is only depending on the two nodes at the ends of the edge and can be determined through where the subsequent shape function is defined as
The relations are used when the path node locates on the edge of an element, or in the selected element, there are two nodes on the RVE face instead of three.
4.3. Verification of the FE Model
It is desirable that the 3D RVE model could be verified with the random inclusions generated using algorithms such as RSA and MC. However, since these algorithms are currently unavailable, one way available to verify the model is to recalculate the transversely isotropic matrix derived using the 2D model in the previous sections. Using regular geometry of inclusion, but free meshing to generate irregular meshing patterns similar to that of random inclusions, the stiffness tensors are recalculated afterwards. If similar results could be obtained, it can be proved that the algorithm for generating 3D model is suitably developed.
In the conducted 2D analysis, when calculating the elastic tensors , , , and , simplified periodic boundary conditions are applied. The two scenarios when applying the original form of periodic boundary conditions are the calculation of elastic tensors and . In these two cases, out-of-plane () and in-plane () shear loads are applied, respectively. These two coefficients will be calculated again using the 3D model. The boundary conditions for calculating and are the same as those listed in Table 1.
4.4. Results and Discussion
Figures 16(a) and 16(b) illustrate strain distribution and stress distribution under in-plane shear loading, respectively, for the calculation of ; Figures 16(c) and 16(d) illustrate strain distribution and stress distribution under out-of-plane shear loading, respectively, for the calculation of . Further, distortion is observed for a great proportion of the RVE. The reason for the distortion is considered to be related to the problematic element selection.
The results are listed in Table 3 together with previous results from 2D model for a comparison. Although the nodal solution distribution has seen distortion, Figure 17 has shown well agreement between the two sets of results with slight differences when volume fraction reaches 0.666 and the error rate generally less than 6%. Admittedly, the contributors to this variation may vary greatly as discussed in Section 3; however, at this stage, the algorithm used herein on 3D model is capable of predicting stiffness properties of composite with random inclusions.
There are several points that should be noted. First of all, the employment of the 10-node tetrahedral element type (Solid227) is a compromise over the irregular inclusion configuration. Theoretically, 20-node hexagonal shaped element type (Solid226) presents a more complex order of shape function along its element edges and more accurate DoFs interpolation than Solid227, so Solid226 is more capable of generating accurate stress and strain field results.
In addition, the 3D analysis requires a much larger number of elements, and the use of Solid227 has again increased the total elements involved; therefore an efficient algorithm is crucial pertaining to the computation capability. To the best knowledge of the author, the final version of both algorithms has been improved, if not optimized. However, as previously mentioned, the number of the division of the RVE edge still has to be controlled under 15 divisions per RVE edge.
5. Conclusions
An algorithm using ANSYS Parametric Design Language to predict effective properties of piezoelectric composite materials has been developed in this work. In the first part, the fibre inclusion within the RVE is of regular geometry and the analysis is simplified as two-dimensional models with effective coefficients. The conclusions can be drawn as follows.(i)With regular geometry, sweep meshing could be adopted to generate the same meshing patterns on opposite areas. The developed algorithm for applying constraint equations on exact opposite node pairs is easy to use.(ii)The calculated transversely isotropic effective coefficients are in good agreement with those in the literatures. The FE results for longitudinal elastic tensor are not very close to those in the literature when the fibre volume fraction is equal to or greater than 0.444.In the second part, the RVE model is developed in such a way that 3D analysis can be conducted for predicting effective properties of composite with random inclusions. Element shape functions are used to interpolate the degree of freedom, so that constraint equations could be applied indirectly. The algorithms for generating random inclusions are not currently available; as a result, the geometry of the RVE in this case is still regular; however, with the use of free meshing, dissimilar meshing patterns are produced for verifying the improved algorithm. It is found that the nodal solution patterns have found distortion. The distortion has been investigated that it is due to the problematic element solution before implementing linear interpolation. Since numerical results have shown validation of the model, the problematic element solution is not supposed to affect the calculation at this stage.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.