Deflection of Piezoelectric Actuator

Open Live Script

This example shows how to solve a coupled elasticity-electrostatics problem.

Piezoelectric materials deform under an applied voltage. Conversely, deforming a piezoelectric material produces a voltage. Therefore, analysis of a piezoelectric part requires the solution of a set of coupled partial differential equations with deflections and electrical potential as dependent variables.

In this example, the model is a two-layer cantilever beam, with both layers made of the same polyvinylidene fluoride (PVDF) material. The polarization direction points down (negativey-direction) in the top layer and points up in the bottom layer. The typical length to thickness ratio is 100. When you apply a voltage between the lower and upper surfaces of the beam, the beam deflects in they-direction because one layer shortens and the other layer lengthens.

The equilibrium equations describe the elastic behavior of the solid:

$- \nabla \cdot σ = f$

Here, $σ$ is the stress tensor, and $f$ is the body force vector. Gauss's Law describes the electrostatic behavior of the solid:

$\nabla \cdot D = ρ$

$D$ is the electric displacement, and $ρ$ 是自由电荷分布。结合these two PDE systems into this single system:

$- \nabla \cdot {\begin{array}{c} σ \\ D \end{array}} = {\begin{array}{c} f \\ - ρ \end{array}}$

For a 2-D analysis, $σ$ has the components $σ_{11}, σ_{22},$ and $σ_{12} = σ_{21}$ , and $D$ has the components $D_{1}$ and $D_{2}$ .

The constitutive equations for the material define the stress tensor and electric displacement vector in terms of the strain tensor and electric field. For a 2-D analysis of an orthotropic piezoelectric material under plane stress conditions, you can write these equations as

${\begin{array}{c} σ_{11} \\ σ_{22} \\ σ_{12} \\ D_{1} \\ D_{2} \end{array}} = [\begin{array}{ccccc} C_{11} & C_{12} & e_{11} & e_{31} \\ C_{12} & C_{22} & e_{13} & e_{33} \\ G_{12} & e_{14} & e_{34} \\ e_{11} & e_{13} & e_{14} & - E_{1} \\ e_{31} & e_{33} & e_{34} & - E_{2} \end{array}] {\begin{array}{c} ϵ_{11} \\ ϵ_{22} \\ γ_{12} \\ - E_{1} \\ - E_{2} \end{array}}$

$C_{i j}$ are the elastic coefficients, $E_{i}$ are the electrical permittivities, and $e_{i j}$ are the piezoelectric stress coefficients. The piezoelectric stress coefficients in these equations conform to conventional notation in piezoelectric materials where thez-direction (the third direction) aligns with the "poled" direction of the material. For the 2-D analysis, align the "poled" direction with they-axis. Write the strain vector in terms of thex-displacement $u$ andy-displacement $v$ :

${\begin{array}{c} ϵ_{11} \\ ϵ_{22} \\ γ_{12} \end{array}} = {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{array}}$

Write the electric field in terms of the electrical potential $ϕ$ :

${\begin{array}{c} E_{1} \\ E_{2} \end{array}} = - {\begin{array}{c} \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

You can substitute the strain-displacement equations and electric field equations into the constitutive equations and get a system of equations for the stresses and electrical displacements in terms of displacement and electrical potential derivatives. Substituting the resulting equations into the PDE system equations yields a system of equations that involve the divergence of the displacement and electrical potential derivatives. As the next step, arrange these equations to match the form required by the toolbox.

Partial Differential Equation Toolbox™ requires a system of elliptic equations to be expressed in a vector form:

$- \nabla \cdot (c \otimes \nabla u) + a u = f$

or in a tensor form:

$- \frac{\partial}{\partial x_{k}} (c_{i j k l} \frac{\partial u_{j}}{\partial x_{l}}) + a_{i j} u_{j} = f_{i}$

where repeated indices imply summation. For the 2-D piezoelectric system in this example, the system vector $u$ is

$u = {\begin{array}{c} u \\ v \\ ϕ \end{array}}$

This is an $N = 3$ system. The gradient of $u$ is

$\nabla u = {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

For details on specifying the coefficients in the format required by the toolbox, see:

The c coefficient in this example is a tensor. You can represent it as a 3-by-3 matrix of 2-by-2 blocks:

$[\begin{array}{cccccc} c (1) & c (2) & c (4) & c (6) & c (11) & c (13) \\ \cdot & c (3) & c (5) & c (7) & c (12) & c (14) \\ \cdot & \cdot & c (8) & c (9) & c (15) & c (17) \\ \cdot & \cdot & \cdot & c (10) & c (16) & c (18) \\ \cdot & \cdot & \cdot & \cdot & c (19) & c (20) \\ \cdot & \cdot & \cdot & \cdot & \cdot & c (21) \end{array}]$

To map terms of constitutive equations to the form required by the toolbox, write the c tensor and the solution gradient in this form:

$[\begin{array}{cccccc} c_{1111} & c_{1112} & c_{1211} & c_{1212} & c_{1311} & c_{1312} \\ \cdot & c_{1122} & c_{1221} & c_{1222} & c_{1321} & c_{1322} \\ \cdot & \cdot & c_{2211} & c_{2212} & c_{2311} & c_{2312} \\ \cdot & \cdot & \cdot & c_{2222} & c_{2321} & c_{2322} \\ \cdot & \cdot & \cdot & \cdot & c_{3311} & c_{3312} \\ \cdot & \cdot & \cdot & \cdot & \cdot & c_{3322} \end{array}] {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

From this equation, you can map the traditional constitutive coefficients to the form required for the c matrix. The minus sign in the equations for the electric field is incorporated into the c matrix to match the toolbox's convention.

$[\begin{array}{cccccc} C_{11} & \cdot & \cdot & C_{12} & e_{11} & e_{31} \\ \cdot & G_{12} & G_{12} & \cdot & e_{14} & e_{34} \\ \cdot & \cdot & G_{12} & \cdot & e_{14} & e_{34} \\ \cdot & \cdot & \cdot & C_{22} & e_{13} & e_{33} \\ \cdot & \cdot & \cdot & \cdot & - E_{1} & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & - E_{2} \end{array}] {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

Beam Geometry

Create a PDE model. The equations of linear elasticity have three components, so the model must have three equations.

model = createpde(3);

Create the geometry and include it in the model.

L = 100e-3;% Beam length in metersH = 1e-3;% Overall height of the beamH2 = H/2;% Height of each layer in meterstopLayer = [3 4 0 L L 0 0 0 H2 H2]; bottomLayer = [3 4 0 L L 0 -H2 -H2 0 0]; gdm = [topLayer;bottomLayer]'; g = decsg(gdm,'R1+R2',['R1';'R2']'); geometryFromEdges(model,g);

Plot the geometry with the face and edge labels.

figure pdegplot(model,'EdgeLabels','on',...'FaceLabels','on') xlabel('X-coordinate, meters') ylabel('Y-coordinate, meters') axis([-.1*L,1.1*L,-4*H2,4*H2]) axissquare

Figure contains an axes object. The axes object contains 10 objects of type line, text.

Material Properties

年代pecify the material properties of the beam layers. The material in both layers is polyvinylidene fluoride (PVDF), a thermoplastic polymer with piezoelectric behavior.

E = 2.0e9;% Elastic modulus, N/m^2NU = 0.29;% Poisson's ratioG = 0.775e9;% Shear modulus, N/m^2d31 = 2.2e-11;% Piezoelectric strain coefficients, C/Nd33 = -3.0 e-11;

年代pecify relative electrical permittivity of the material at constant stress.

relPermittivity = 12;

年代pecify the electrical permittivity of vacuum.

permittivityFreeSpace = 8.854187817620e-12;% F/mC11 = E/(1 - NU^2); C12 = NU*C11; c2d = [C11 C12 0; C12 C11 0; 0 0 G]; pzeD = [0 d31; 0 d33; 0 0];

年代pecify the piezoelectric stress coefficients.

pzeE = c2d*pzeD; D_const_stress = [relPermittivity 0; 0 relPermittivity]*permittivityFreeSpace;

Convert the dielectric matrix from constant stress to constant strain.

D_const_strain = D_const_stress - pzeD'*pzeE;

You can view the 21 coefficients as a 3-by-3 matrix of 2-by-2 blocks. Thecijmatrices are the 2-by-2 blocks in the upper triangle of this matrix.

c11 = [c2d(1,1) c2d(1,3) c2d(3,3)]; c12 = [c2d(1,3) c2d(1,2); c2d(3,3) c2d(2,3)]; c22 = [c2d(3,3) c2d(2,3) c2d(2,2)]; c13 = [pzeE(1,1) pzeE(1,2); pzeE(3,1) pzeE(3,2)]; c23 = [pzeE(3,1) pzeE(3,2); pzeE(2,1) pzeE(2,2)]; c33 = [D_const_strain(1,1) D_const_strain(2,1) D_const_strain(2,2)]; ctop = [c11(:); c12(:); c22(:); -c13(:); -c23(:); -c33(:)]; cbot = [c11(:); c12(:); c22(:); c13(:); c23(:); -c33(:)]; f = [0 0 0]'; specifyCoefficients(model,'m',0,'d',0,'c',ctop,'a',0,“f”,f,'Face',2); specifyCoefficients(model,'m',0,'d',0,'c',cbot,'a',0,“f”,f,'Face',1);

Boundary Conditions

年代et the voltage (solution component 3) on the top of the beam (edge 1) to 100 volts.

voltTop = applyBoundaryCondition(model,'mixed',...'Edge',1,...'u', 100,...'EquationIndex'3);

年代pecify that the bottom of the beam (edge 2) is grounded by setting the voltage to 0.

voltBot = applyBoundaryCondition(model,'mixed',...'Edge',2,...'u',0,...'EquationIndex'3);

年代pecify that the left side (edges 6 and 7) is clamped by setting thex- andy-displacements (solution components 1 and 2) to 0.

clampLeft = applyBoundaryCondition(model,'mixed',...'Edge',6:7,...'u'[0 0],...'EquationIndex',1:2);

The stress and charge on the right side of the beam are zero. Accordingly, use the default boundary condition for edges 3 and 4.

有限元分析解决方案金宝搏官方网站

Generate a mesh and solve the model.

msh = generateMesh(model,'Hmax',5e-4); result = solvepde(model)

result = StationaryResults with properties: NodalSolution: [3605x3 double] XGradients: [3605x3 double] YGradients: [3605x3 double] ZGradients: [0x3 double] Mesh: [1x1 FEMesh]

Access the solution at the nodal locations. The first column contains thex-deflection. The second column contains they-deflection. The third column contains the electrical potential.

rs = result.NodalSolution;

Find the minimumy-deflection.

feTipDeflection = min(rs(:,2)); fprintf('Finite element tip deflection is: %12.4e\n',feTipDeflection);

Finite element tip deflection is: -3.2900e-05

Compare this result with the known analytical solution.

tipDeflection = -3*d31*100*L^2/(8*H2^2); fprintf('Analytical tip deflection is: %12.4e\n',tipDeflection);

Analytical tip deflection is: -3.3000e-05

Plot the deflection components and the electrical potential.

varsToPlot = char('X-Deflection, meters',...'Y-Deflection, meters',...'Electrical Potential, Volts');fori = 1:size(varsToPlot,1) figure; pdeplot(model,'XYData',rs(:,i),'Contour','on') title(varsToPlot(i,:))% scale the axes to make it easier to view the contoursaxis([0, L, -4*H2, 4*H2]) xlabel('X-Coordinate, meters') ylabel('Y-Coordinate, meters') axissquareend

Figure contains an axes object. The axes object with title X-Deflection, meters contains 12 objects of type patch, line.

Figure contains an axes object. The axes object with title Y-Deflection, meters contains 12 objects of type patch, line.

Figure contains an axes object. The axes object with title Electrical Potential, Volts contains 12 objects of type patch, line.

References

Hwang, Woo-Seok, and Hyun Chul Park. "Finite Element Modeling of Piezoelectric Sensors and Actuators."AIAA Journal31, no.5 (May 1993): 930-937.
Pieford, V. "Finite Element Modeling of Piezoelectric Active Structures." PhD diss., Universite Libre De Bruxelles, 2001.