This example shows how to solve a coupled elasticity-electrostatics problem using the Partial Differential Equation Toolbox™. Piezoelectric materials deform when a voltage is applied. Conversely, a voltage is produced when a piezoelectric material is deformed.
Analysis of a piezoelectric part requires the solution of a set of coupled partial differential equations with deflections and electrical potential as dependent variables. One of the main objectives of this example is to show how such a system of coupled partial differential equations can be solved using PDE Toolbox.
PDE For a Piezoelectric Solid
The elastic behavior of the solid is described by the equilibrium equations
where is the stress tensor and is the body force vector. The electrostatic behavior of the solid is described by Gauss' Law
where is the electric displacement and is the distributed, free charge. These two PDE systems can be combined into the following single system
In 2D, has the components and and has the components and .
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 2D, orthotropic, piezoelectric material under plane stress conditions these are commonly written as
where are the elastic coefficients, are the electrical permittivities, and are the piezoelectric stress coefficients. The piezoelectric stress coefficients are written to conform to conventional notation in piezoelectric materials where the z-direction (3-direction) is aligned with the "poled" direction of the material. For the 2D analysis, we want the poled direction to be aligned with the y-axis.
Finally, the strain vector can be written in terms of the x-displacement, , and y-displacement, as
and the electric field written in terms of the electrical potential, , as
See reference 2, for example, for a more complete description of the piezoelectric equations.
The strain-displacement equations and electric field equations above can be substituted into the constitutive equations to yield a system of equations for the stresses and electrical displacements in terms of displacement and electrical potential derivatives. If the resulting equations are substituted into the PDE system equations, we have a system of equations that involve the divergence of the displacement and electrical potential derivatives. Arranging these equations to match the form required by PDE Toolbox will be the topic for the next section.
Converting the Equations To PDE Toolbox Form
The PDE Toolbox requires a system of elliptic equations to be expressed in the form
or in tensor form
where summation is implied by repeated indices. For the 2D piezoelectric system described above, the PDE Toolbox system vector is
This is an system. The gradient of is given by
The documentation for the function
assempde shows that it is convenient to view the tensor
submatrices. The most convenient form for the
input argument for this symmetric,
system has 21 rows in
and is described in detail in the PDE Toolbox documentation. It is repeated here for convenience.
For the purposes of mapping terms from constitutive equations to the form required by PDE Toolbox it is useful to write the tensor and solution gradient in the following form
From this equation the traditional constitutive coefficients can be mapped to the form required for the PDE Toolbox matrix. Note the minus sign in the equations for electric field. This minus must be incorporated into the matrix to match the PDE Toolbox convention. This is shown explicitly below.
Piezoelectric Bimorph Actuator Model
Now that we have defined the equations for a 2D piezoelectric material, we are ready to apply these to a specific model. The model is a two-layer cantilever beam that has been extensively studied (e.g. refs 1 and 2). It is defined as a "bimorph" because although both layers are made of the same Polyvinylidene Fluoride (PVDF) material, in the top layer the polarization direction points down (minus y direction) and in the bottom layer, it points up. A schematic of the cantilever beam is shown in the figure below.
This figure is not to scale; the actual thickness/length ratio is 100 so the beam is very slender. When a voltage is applied between the lower and upper surfaces of the beam, it deflects in the y-direction; one layer shortens and the other layer lengthens. Devices of this type can be designed to provide the required motion or force for different applications.
Create a PDE Model with three dependent variables
numberOfPDE = 3; pdem = createpde(numberOfPDE);
Geometry and Mesh
The simple two-layer geometry of the beam can be created by defining the sum of two rectangles.
L = 100e-3; % beam length in meters H = 1e-3; % overall height of the beam H2 = H/2; % height of each layer in meters % The two lines below contain the columns of the % geometry description matrix (GDM) for the two rectangular layers. % The GDM is the first input argument to decsg and describes the % basic geometric entities in the model. topLayer = [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']'); % Create a geometry entity and append to the PDE Model geometryFromEdges(pdem,g); figure; pdegplot(pdem, 'edgeLabels', 'on'); title 'Two-layer Piezoelectric Cantilever Beam (with edge labels)' xlabel 'X-coordinate, meters'; ylabel 'Y-coordinate, meters'; axis([-.1*L, 1.1*L, -4*H2, 4*H2]); % We need a relatively fine mesh with maximum element size roughly equal H/16 % to accurately model the bending of the beam. hmax = H/16; msh = generateMesh(pdem, 'Hmax', hmax, 'MesherVersion', 'R2013a');
Warning: Approximately 51200 triangles will be generated.
The material in both layers of the beam is Polyvinylidene Fluoride (PVDF), a thermoplastic polymer with piezoelectric behavior.
E = 2.0e9; % Elastic modulus, N/m^2 NU = 0.29; % Poisson's ratio G = 0.775e9; % Shear modulus, N/m^2 d31 = 2.2e-11; % Piezoelectric strain coefficients, C/N d33 = -3.0e-11; % relative electrical permittivity of the material relPermittivity = 12; % at constant stress % electrical permittivity of vacuum permittivityFreeSpace = 8.854187817620e-12; % F/m C11 = E/(1-NU^2); C12 = NU*C11; c2d = [C11 C12 0; C12 C11 0; 0 0 G]; pzeD = [0 d31; 0 d33; 0 0]; % The piezoelectric strain coefficients for PVDF are % given above but the constitutive relations in the % finite element formulation require the % piezoelectric stress coefficients. These are calculated on the next % line (for details see, for example, reference 2). pzeE = c2d*pzeD; D_const_stress = [relPermittivity 0; 0 relPermittivity]*permittivityFreeSpace; % Convert dielectric matrix from constant stress to constant strain D_const_strain = D_const_stress - pzeD'*pzeE; % As discussed above, it is convenient to view the 21 coefficients % required by assempde as a 3 x 3 array of 2 x 2 submatrices. % The cij matrices defined below are the 2 x 2 submatrices in the upper % triangle of this array. 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)];
Function To Return C Coefficients
The c-matrix for this N = 3 system is symmetric. From the documentation for
assempde, we see that the most convenient form for defining the c-matrix has 21 rows defining the upper triangle of the matrix.
c = @(p, t, u, t0) calcCMatPiezoActuator(p, t, c11, c12, c22, c13, c23, c33); % The function shown below is called by the PDE Toolbox routines to % return the required 21 entries in the c-matrix. type calcCMatPiezoActuator
function c = calcCMatPiezoActuator( p, t, c11, c12, c22, c13, c23, c33 ) %CALCCMATPIEZOACTUATOR C-matrix for piezoelectric actuator example % c = CALCCMATPIEZOACTUATOR( p, t, c11, c12, c22, c13, c23, c33 ) % returns the 'c' coefficient matrix for the piezoelectric actuator % example given the point and element matrices along with the % constitutive submatrices (cij) for the PVDF material. % Copyright 2012 The MathWorks, Inc. numElems = size(t,2); c=zeros(21,numElems); % % Although the material in both layers is PVDF, in the top layer % the polarization direction points down (minus y direction) and in the % bottom layer, it points up. That is, the top layer has d-coefficients % that are the negative of those in the bottom layer. % % The code below examines the y-location of the centroid of each % triangular element and assigns the correct material properties to % element depending on whether it is in the top or bottom layer. % ctop = [c11(:); c12(:); c22(:); -c13(:); -c23(:); -c33(:)]; cbot = [c11(:); c12(:); c22(:); c13(:); c23(:); -c33(:)]; % calculate y-coordinate of triangle centers yCenter=(p(2,t(1,:)) + p(2,t(2,:)) + p(2,t(3,:)))/3; for i=1:numElems if(yCenter(i) < 0) c(:,i) = cbot; else c(:,i) = ctop; end end end
Boundary Condition Definition
For this example, the top geometry edge (edge 1) has the voltage prescribed as 100 volts. The bottom geometry edge (edge 2) has the voltage prescribed as 0 volts (i.e. grounded). The left geometry edge (edges 6 and 7) have the u and v displacements equal zero (i.e. clamped). The stress and charge are zero on the right geometry edge (i.e. q = 0).
V = 100; % Set the voltage (solution component 3) on the top edge to V. voltTop = applyBoundaryCondition(pdem,'Edge',1, 'u', V, 'EquationIndex', 3); % Set the voltage (solution component 3) on the bottom edge to zero. voltBot = applyBoundaryCondition(pdem,'Edge',2, 'u', 0, 'EquationIndex', 3); % Set the x and y displacements (solution components 1 and 2) % on the left end (geometry edges 6 and 7) to zero. clampLeft = applyBoundaryCondition(pdem,'Edge',6:7, 'u', [0 0], 'EquationIndex', 1:2);
Finite Element Solution
a = 0; f = [0 0 0]'; u = assempde(pdem, c, a, f); % % For display and plotting purposes, it is convenient to reshape the % solution vector as three columns containing the x-displacement, % y-displacement, and electrical potential, respectively. % uu = reshape(u, , numberOfPDE); feTipDeflection = uu(1,2); fprintf('Finite element tip deflection is: %12.4e\n', feTipDeflection); varsToPlot = char('X-Deflection, meters', 'Y-Deflection, meters', ... 'Electrical Potential, Volts'); for i = 1:size(varsToPlot,1) figure; pdeplot(pdem, 'xydata', uu(:,i), 'contour', 'on'); title(varsToPlot(i,:)); % scale the axes to make it easier to view the contours axis([0, L, -4*H2, 4*H2]); xlabel 'X-Coordinate, meters' ylabel 'Y-Coordinate, meters' end
Finite element tip deflection is: -3.2772e-05
A simple, approximate, analytical solution was obtained for this problem in reference 1.
tipDeflection = -3*d31*V*L^2/(8*H2^2); fprintf('Analytical tip deflection is: %12.4e\n', tipDeflection);
Analytical tip deflection is: -3.3000e-05
The color contour plots of x-deflection and y-deflection show the standard behavior of the classical cantilever beam solution. The linear distribution of voltage through the thickness of the beam is as expected. There is good agreement between the PDE Toolbox finite element solution and the analytical solution from reference 1.
Although this example shows a very specific coupled elasticity-electrostatic model, the general approach here can be used for many other systems of coupled PDEs. The key to applying PDE Toolbox to these types of coupled systems is the systematic, multi-step coefficient mapping procedure described above.
Hwang, W. S.; Park, H. C; Finite Element Modeling of Piezoelectric Sensors and Actuators. AIAA Journal, 31(5), pp 930-937, 1993.
Pieford, V; Finite Element Modeling of Piezoelectric Active Structures. PhD Thesis, Universite Libre De Bruxelles, 2001.