1. Introduction
Piezoelectric materials can be found in a wide variety of products from kitchen appliances to automobile industries and health care devices. The extensive use of piezoelectric materials in many engineering applications is due to its unique physical characteristics resulting from its specific material properties. Piezoelectric materials are crystalline materials in which mechanical and electrical states interact with each other. The mechanical and electrical states interact in a material with charged particles with no center of symmetry. Three types of microscopic level charge disorientation produce the macroscopic level polarization in a dielectric material, including electric and ionic polarization and dipole reorientation [
1]. In the electronic polarization, the applied electric field deforms the electronic cloud, but in the case of ionic polarization, the applied electric field pushes the anion to the anode and the cations to the cathode. In the dipole reorientation, an applied electric field reorients the dipole directions in the separate domain of the material in one direction to polarize the material. The latter method is valid for the macroscopic non-polarized material consisting of many polarized volumes. A widely-used class of piezoelectric materials is a polycrystalline structure of lead zirconate titanate, Pb(Zr
Ti
)O
, also known as PZT which is a solid solution of two materials, lead titanate PbZrO
and lead zirconate PbTiO
. The Zr and Ti ions in a lattice of an unpolarized PZT material are randomly distributed to neutralize the electrical polarity. The application of a high electric field at a temperature close to the Curie point aligns the Zr and Ti ions in lattices to polarize the materials. In another class of piezoelectric materials, molecular dipoles of polar polymers such as Polyvinylidene Fluoride (PVDF) are reoriented in one direction.
The ability of piezoelectric materials to generate electric potential from a mechanical strain and vice versa provides opportunities to use these materials as sensors or actuators found in many types of transducers. A sensor transforms mechanical energy to electrical output for monitoring purposes, and an actuator normally converts electrical energy to mechanical motions. The widely-discussed use of piezoelectric materials as sensors is found in the sonar (sound navigation and ranging) system, a sound-based navigation system, in which the electro-acoustic transducer converts the acoustic signals to the electrical signals. Ultrasound equipment in the medical industry uses the piezoelectric actuator to convert the electrical signal to mechanical vibration. A quartz wrist watch uses a piezoelectric actuator to convert electrical energy from a battery to mechanical vibration, which moves the second, minute, and hour arms of the watch. More recently, researchers have been exploring piezoelectric materials in their innovations such as heart rate monitors [
2], piezoelectric-based toe-heaters for frostbite protection [
3], and conversion of ocean wave energy to electrical energy [
4]. Piezoelectric materials are also widely used in many droplet-on-demand applications including inkjet head printers, spraying devices, etc. For its wide ranging applications, understanding the characteristics of piezoelectric materials is essential in the process of the design and optimization of those devices.
Characterizing piezoelectric material behaviors is typically done through experiment with a specific design and prototype. Besides experimental characterization, there exist different analytical and numerical methods to model and analyze the responses and characteristics of piezoelectric material-based devices. Equivalent circuits [
5], spring models [
6], and thermal analogy [
7] are analytical methods to model the simplified piezoelectric materials. For three-dimensional analysis, the finite element method has been a popular choice for modeling and simulations of piezoelectric material responses. Since the pioneering work of [
8] on the application of the finite element method for piezoelectric materials, there have been numerous development of FEM tools to model piezoelectric structures ranging from solid three-dimensional elements to one- and two-dimensional elements such as beam, shell, and plates; including recent work on developing linear and quadratic shell elements [
9,
10]. In [
11], a comprehensive review of the development of FEM was presented, classifying the number of different approaches based on element types, as well as the number of degrees of freedom. While the majority of solid three-dimensional FEM models focus on static and modal analysis of piezoelectric structures, there are fewer developments for unsteady analysis. This is even more desirable in Fluid-Structure Interaction (FSI) simulations when piezoelectric structure models are required to couple with flow models for analysis of FSI responses.
With this motivation, this paper presents an implementation of the Finite Element Method (FEM) to solve the piezoelectric problems and its applications to the design of piezoelectric actuators. In this work, the FEM implementation is based on the standard nodal approach extended for piezoelectric materials using three-dimensional tetrahedral and hexahedral elements. While hexahedral elements are mostly preferred in FEM analysis due to their superior accuracy and robustness, tetrahedral elements allow handling of complex geometries often encountered in practical applications. To account for electrical-mechanical coupling in piezoelectric materials, an additional degree of freedom for electrical potential is added to each node in those elements together with their standard displacement unknowns. The implementation of the piezoelectric model is based on an existing open source FEM code, OOFEM [
12]. This presents the first significant contribution of the current work to the development of three-dimensional FEM as a general purpose analysis tool for piezoelectric materials.
For any numerical analysis, it is important to quantify errors associated with a model for quality assessment of its effectiveness and predictability. Apart from known errors present in the model, there are uncertainties from various sources including input parameters and material properties. There has been a great emphasis on verification and validation together with uncertainty quantification for numerical computations over the past few decades, including the pioneering work on the introduction of the Grid Convergence Index (GCI) [
13] to a more general framework presented in [
14,
15]. Carrying out the error analysis and uncertainty quantification certainly provide confidence in the predictability of computational models. However, this process is cumbersome and time consuming, thus discouraging practitioners from adhering to those guidelines. In this work, a detailed framework of verification and validation was followed to access the convergence, as well as the accuracy of the proposed approach. Errors from uncertain parameters are also considered to establish bounds of prediction and investigate sources of discrepancy between numerical and experimental data. Uncertainties from different sources are identified, characterized, and quantified by propagating them through the model using the general framework proposed in [
15].
The rest of the paper is organized as follows. In the next section, we present the mathematical foundation for the piezoelectric material modeling including a discussion on the finite element formulation for the dynamic piezoelectric problem. The numerical method is then validated with experimental data and analytical results for a number of problems ranging from a simple unimorph beam to three-layer plates. The last section of the paper presents an application of the developed method for analyzing the responses of piezoelectric actuators used in synthetic jet devices. It includes validation with experimental data, as well as quantification of uncertainties in material properties on response predictions.
4. Uncertainty Quantification of the Piezoelectric Composite Plate Actuator’s Response
Composite plate piezoelectric actuators are found in many electroacoustic devices, including microphones, micropumps, synthetic jets, etc. The actuator comprises a circular plate of piezoelectric material bonded with a metallic material plate, and actuation is achieved by applying electrical potential on an electrode attached to the piezoelectric layer.
Figure 8a shows a typical setup of a composite piezoelectric unimorph actuator (APC 850) [
18,
19]. In this current setup, a piezoelectric plate of thickness
and radius
was bonded with a shim disc of thickness
and radius
. The shim plate was fully clamped, and the material properties of APC 850 were taken from
Table 5 extracted from [
18,
20]. The computational domain was taken as one quarter of the composite plate in which symmetric conditions were applied at the cut planes.
Figure 8b shows a block-structured mesh of the computational domain where it is characterized by
number of points in the radial and azimuthal directions of blocks on the piezoelectric plate,
number of points (radial and azimuthal) on the blocks of the shim, and
number of layers extruded in the thickness direction for each plate.
4.1. Characterization and Comparison of Actuator Responses
It is important to characterize the responses of the actuator for design purposes. Prasad et al. [
18,
19] presented several models for response predictions of unimorph actuators including a lumped element model and the analytical model based on plate element theory. It is known that the maximum deflection of the actuator is dependent on the plate thickness ratio (
), radius (
), and material properties. Here, FEM was applied for static analysis of the actuator under voltage load of 5 V. A mesh sensitivity study is shown in
Figure 9 for the radial deflection profile with clear convergence of numerical prediction when the mesh is refined.
Table 6 shows the convergence index with the mesh refinement. At the finest mesh, the maximum displacement was about
away from the asymptotic value of the predicted displacement. The order of convergence estimated according to Equation (
20) is:
which is approximately the same as the theoretical convergence order of the FEM formulation. It can be seen that the asymptotic FEM result was larger compared with the experimental and analytical data in [
18]. As explained in [
18], the difference could be due to the fact that the silver electrode layer, as well as the glue between plates were not modeled in the FEM analysis. Looking further into the experimental data, it was found that the device (APC 850) used in a later work by the same author [
19] showed a slightly different response experimentally and analytically.
Figure 10 shows experimental data from two different references for the same reported device. It can be seen that the difference on deflection of the disc was very large between the two configurations where maximum deflection at the center of the disc increased
from 0.0937
m–0.1184
m per unit applied voltage. For these two configurations, the only significant difference in the setup was the thickness of the shim layers. The current FEM results showed an increase in response amplitude with reducing shim thickness, as observed in the experiment. This is further illustrated in
Figure 11, where maximum displacements are shown with changing of the shim and plate thickness ratio in the variation of their radius. The maximum displacement also increased with increasing radius ratio. This behavior could help in explaining the discrepancies between current FEM analysis with experimental data. However, it led to uncertainties in the analysis, requiring detailed quantification of actuator responses under uncertain parameters.
4.2. Uncertainty Quantification of Actuator Responses
From the comparison between numerical prediction with experimental data for maximum displacement of the actuator, it can be seen that the error was consistently about 6.3%. The purpose of this analysis was to quantify uncertainty in the current numerical prediction in order to explain the discrepancy. Typically, the process of uncertainty quantifications includes the following steps [
14]: identification of all sources of uncertainties, characterization of input uncertainties, estimation of numerical errors, propagation of input uncertainties through the model, evaluation of model from error, and finally, the total uncertainty.
4.2.1. Sources of Uncertainties
There were two main sources of uncertainty in the model inputs, and they include the material properties and geometrical parameters of the actuator. From the supplier report on material properties of APC 850 [
20], it is highlighted that there were uncertainties in the physical properties of the disc depending on the shapes and thickness of the plates. For some quantities such as piezoelectric coefficient
, the error could be
from the reported value. The thickness of the PZT and shim layer was also a source of uncertainty due to the error in the manufacturing process, as well as uncertainties with measurement. While uncertainty in thickness can be classified as aleatory (irreducible uncertainty), the uncertainty in material properties is reducible with more refined measurement techniques. Aleatory uncertainty is normally characterized by a Probability Density Function (PDF), and epistemic uncertainty can be represented as an interval-valued quantity.
Uncertainties were introduced to four input parameters of shim plate thickness (
), PZT plate thickness (
), and piezoelectric coefficient (
,
). In this study, the responses of the actuator measured by the maximum deflection were investigated with the above uncertainty. The piezoelectric properties of materials were normally found to be in a wide range of values with the same material, especially for the piezoelectric coefficient in this case.
Table 7 presents a list of uncertain variables with their mean values, as well as the probability distribution. Here, uncertainties in physical properties and geometrical parameters are modeled as a random variable with a normal distribution characterized by the mean and standard deviation.
Sensitivity analysis was conducted for maximum center deflection with uncertain variables, as shown in
Figure 12, using the Latin hypercube sampling approach [
21]. It can be seen from the scatter plots that the response maximum amplitude was more correlated with thicknesses and piezoelectric coefficient
, while correlation with other parameter
was weak. This is evidently shown in the results of the correlation coefficients for all the parameters in
Table 7, where the strongest correlation with
was observed for different approaches and correlation measures.
4.2.2. Quantification of Numerical Errors
The main numerical error in this current computation was discretization error. To evaluate the discretization error in the Uncertainty Quantification (UQ) framework, numerical simulations were performed on two grid levels for each of the samples. The results from the coarse and fine grids can be used to estimate discretization errors using Richardson’s extrapolation, as stated in Equation (
21). Here, for each sampled case, this was carried out on two grid levels of
and
(
Table 6).
Figure 13 shows the distribution of the maximum displacement and estimated numerical error for 200 samples using the Latin hypercube sampling technique. The mean and standard deviation of numerical uncertainty was
× 10
(
m/V) and
× 10
(
m/V), respectively.
4.2.3. Propagation of Uncertainties through the Model
Depending on the types of uncertainties, there are different approaches to quantify uncertainty in model predictions [
15]. Sampling methods can be used for Uncertainty Quantification (UQ) in which the response characteristic can be obtained from a probability distribution of uncertain inputs. Apart from this sampling-based approach, other UQ methods including polynomial chaos expansion, interval analysis, and reliability were explored to determine the uncertainty of the response.
Table 8 presents the moments of the response (maximum deflection) obtained from using different UQ methods. It can be seen that the mean response predicted by all the methods was consistent, ranging from 0.1031–0.1054
m/V, while the prediction of the standard deviation was more spread using different methods. The global interval analysis approach normally used for quantification of epistemic uncertainties provided a similar response prediction as the LHS-based approach, with much fewer evaluations.
Figure 14 depicts the Cumulative Distribution Function (CDF) of maximum deflection where it is less likely (
) to obtain a response of
, while the probability of obtaining maximum deflection of
m/V is very likely (80–90%). Results from this UQ confirmed the variations in the responses measured from experiments and provided confidence to the numerical simulations.
4.2.4. Model Uncertainty
The model uncertainty is measured from the error between the model prediction of the quantity of interest and the true value obtained from experiment observation at the same input conditions. With uncertainties in model input, a CDF of the quantity of interest can be obtained when propagating those uncertainties through the model, as shown in
Figure 14. Likewise, one can also reconstruct a CDF of experimental observation taking into account the uncertainties with measurement procedures. Here, a crude reconstruction of experimental CDF was performed from two available data points. Following [
22], a model validation metric was defined as the intersection area between the two CDF curves:
Among various measures of model errors, the area validation metric provides a good estimation for evidence of disagreement between the model and experiment. It can be computed from a limited number of experimental data; thus making it more suitable for many engineering applications. From the current results of CDF shown in
Figure 14, the model error is estimated to be
× 10
m/V.
From the above estimation of uncertainties, for any given condition, a corrected prediction of maximum displacement can be computed by including both model error and numerical uncertainties:
This corrected prediction takes into account various uncertainties in model input, numerical discretization errors, as well as model errors.