In this blog post we discuss piezoelectric transducers and the best way to model them with finite element analysis (FEA).
What is a piezoelectric transducer?
Piezoelectric transducers are a special type of transducer which use the piezoelectric effect to convert mechanical deformation into an electric voltage (or the reverse).
This piezoelectric effect is very useful in ultrasonic sensors which are one type of piezoelectric transducer. Bats use the vibration of its vocal cords to generate ultrasound waves. Those waves bounce on objects and come back, indicating the location of the object to the bat. In an ultrasonic piezoelectric transducer, the wave is generated not by vocal cords, but by a vibrating piezoelectric crystal excited by an electric signal.
In which industries are piezoelectric transducers used?
Piezoelectric transducers are used heavily in sonar, medical imaging, NDE (NonDestructive Engineering) and signal processing systems.
Here is an example of an ultrasonic piezoelectric transducer used for medical imaging.
Those transducers are currently used for diagnostic imaging, Doppler velocity measurement and many other specialty applications (intracavity, biopsy, etc.) including disease treatment (lithotripsy, hyperthermia, tissue ablation).
In the last 20 years, engineers in the ultrasound industry have done a remarkable job developing and refining these devices using a combination of semianalytical design procedures and prototype experiments.
What is the traditional way to simulate those piezoelectric transducers?
The traditional way to design those sonar, medical imaging and NDE systems is to use 1D models and a lot of experimental prototypes.
Those methods are limited and 2D/3D Finite Element modeling can offer a lot more insights in the design of those systems.
Unfortunately, the simulation of those systems using the traditional implicit algorithms for frequency or timedomain integration has limited a lot the scope.
Why is implicit algorithm a limiting factor for the simulation of such devices?
In general, implicit algorithms are best suited to linear static problems, steady state vibrations, and low frequency dynamics. A much better choice for transient phenomena, linear or nonlinear, is to use an explicit time domain algorithm.
Said simply, the implicit algorithm requires the resolution of a costly system of equation at each time step, making the computation really costly for systems which require a very small time step to obtain an accurate solution… such as wave propagation problems in ultrasonic transducers.
Side Note:
OnScale uses a mixed explicit/implicit algorithm which performs a direct time integration of the 2D or 3D electromechanical equations. This special method increases the speed at which such 2D or 3D piezoelectric can be calculated while maintaining a very good level of results accuracy. The computational advantages over conventional implicit algorithms is typically a factor of 100 or more in speed and 1001000 in model size on a same workstation. The cloud amplify increases even the speed even more by providing an almost infinite amount of calculation resources scalable at will.
How is the finite element method applied and what are the assumptions taken?
The finite element method reduces the electromechanical partial differential equations (PDEs) over the model domain to a system of ordinary differential equations (ODEs) in time that can be solved.
For that, various theorems and principles are used such as:
• The principle of virtual work
• Weak form equations
• The Galerkin’s method
• Weighted residuals
• Conservation and balance laws
• Etc…
Behavior of a single flextensional model driven in water, showing radiating wave pattern (left). View of mode shape in the single flextensional model driven at 4.5 kHz (right).
The result is that spatial derivatives in the PDEs are reduced to a summation of “elemental” systems of linear algebraic equations on the unknown field values at nodes of the finite element discretization.
The finite elements used are 4 node quadrilaterals in 2D and 8 nodes hexahedrons in 3D and the unknown field over an element is represented by low order shape functions determined by nodal (corner) values, i.e., bilinear in 2D and trilinear in 3D.
Using a minimum of 15 elements per wavelength limits wave dispersion errors to less than 1%.
Of course, experiment based validation has been used to demonstrate that those choices offer the most robust basis for large scale wave propagation analysis in structural and isotropic or anisotropic continuum models.
What is the system of ordinary differential equations (ODEs) which is solved?
The electromechanical finite element equations are derived from the piezoelectricity constitutive relations and the equations of mechanical and electrical equilibrium
The global system of ODEs is the following:
Equation (1), (2) and (3) govern the couple behavior of the elastic, electric and acoustic fields.
Global unknowns u, Phi , and Psi are, respectively, the elastic displacement vector, the electric potential vector, and the velocity potential vector, while F is the applied force vector and Q is the charge vector. These vectors are defined by field values at all nodes in the model.
Coefficients M, C, and K denote the various uncoupled and coupled “mass,” “damping,” and “stiffness” matrices, respectively.
Should we solve those equations in the Frequency domain or in the Time domain?
Frequency domain solutions and Time domain solutions have different ways to consider those equations.
Let’s resume those differences in the following table:
Frequency Domain Solutions  Time Domain Solutions 
Assume timeharmonic behavior, effectively removing time as an independent variable. Because of that, the unknowns become: This leads to a system of equations independent of time. This leads to a system of equations independent of time. 
Timedomain solutions assume general Integration can be done using either an implicit or an explicit method

Let’s consider now transient signal and time integration.
For such time domain solutions, the time cannot be considered as an independent variable and removed from the equations (like for the frequency domain solutions)
It requires a stepbystep integration in time.
It means that we will need to find a way to calculate the current time step from the previous time steps.
For that, Implicit methods and Explicit Methods can both be used.
What is the difference between using an Implicit or an Explicit integration Scheme?
Implicit Method  Implicit Method 
Implicit methods couple the current solution vector, hence, the global system of equations must be solved at each time step. Their advantage is unconditional stability with respect to time step. The disadvantage for such Implicit integration of wave phenomena is that solution accuracy requires a time step smaller than onetenth the period of the highest frequency to be resolved. 
Explicit methods decouple the current solution vector and eliminate the global system solve

*This means that there is a time step limit (CFL condition) above which the method is unstable.
If Implicit is unconditionally stable with respect to time step, why can we still use explicit?
The answer is very specific to wave propagation problems:
For such problems, the Implicit integration scheme requires a time step smaller than onetenth the period of the highest frequency to be resolved and provide accurate results.
This is close to the CFL stability limit for explicit methods and effectively removes the principal advantage of implicit integration.
How does the Explicit integration scheme work?
Explicit integration of (1) and (3) involves:
1 Diagonalizing the uncoupled mass and damping matrices, Muu, Myy , Cuu , Cyy
2 Using nodal lumping, replacing the time derivatives with finite differences
3 Integrating using a central difference scheme (2nd order accurate)
Then, we must respect the CFL condition which states that:
For stability the time step must be smaller than the shortest wave transit time across any element.
This follows from the hyperbolic (wave) nature of the original PDEs which provides the following benefit:
During time step Dt the field at a point is only influenced by the field at neighboring points within a sphere of radius Dx=cpDt where cp is the fastest local
wave speed. Therefore, for Dt<Dx/cp nodal fields are decoupled during a single time step and can be integrated independently.
Why is this faster than solving with a “more traditional” Implicit scheme ?
The disadvantage for such Implicit integration of wave phenomena is that solution accuracy requires a time step smaller than onetenth the period of the highest frequency to be resolved and the global system of equations must be solved at each timestep.
This is computationally HUGE.
With the explicit scheme, it is possible to eliminate the manipulation and solution of large systems of electromechanical equations by integrating (1) and (3) explicitly in time for u and Psi and using Phi from the previous time step.
The new updated Phi is obtained from solving (2) implicitly for the new Phi from u using a preconditioned CG iteration (diagonal scaling).
Thus, the algorithm operates on an elementbyelement basis, where elemental contributions are accumulated in intermediate global vectors.
The nodal force or charge vectors, and the algorithm processes these vectors only. This is equivalent to matrixvector arithmetic without actually forming the matrix, which facilitates vectorization and parallelization.
How about damping? Can it be modeled correctly in timedomain?
Actually, damping is inherently a timedomain phenomenon, because real world is time domain.
When damping is defined in frequency domain, it is basically nothing more than a simple multiplicator coefficient CHOSEN to provide the wanted damping result.
Because of that, it is not possible to model damping which depends on frequency in frequency domain (that would make the system of equations nonlinear and unsolvable) true?
What are the most used Damping models?
An Important issue in transducer modeling is exactly this frequency dependent material damping.
The two most convenient damping models are:
 Massproportional (viscous, µ1 /w )
 Stiffnessproportional ( µw )
Rayleigh damping is a linear combination of both.
In the time domain, constant coefficients yield damping that is inversely or directly proportional to frequency, or a linear combination.
Side Note:
OnScale also has a materialdependent, threeparameter viscoelastic damping model.
Is the choice of the Viscosity constants and relaxation time important?
Proper choice of viscosity constants and a relaxation time yields a damping maximum at the selected frequency and smooth falloff.
Viscoelastic models may be superposed to yield a discrete spectrum of relaxation times that represent specified damping behavior
over a limited frequency range, but at significant cost in memory.