FEATool Multiphysics
v1.14
Finite Element Analysis Toolbox

This example describes how to model piezoelectric bending of a beam on the command line as an mscript model in FEATool. The test case comes from a research paper by Hwang and Park [3]. For this problem, the following coupled system of equations should be solved
\[ \nabla\cdot \left[\begin{array}{c} \sigma \\ D \end{array}\right] = \left[\begin{array}{c} f \\ \rho \end{array}\right] \]
where \(\sigma\) is a stress tensor, \(f\) represents body forces, \(D\) electric displacement, and \(\rho\) distributed charges.
By using the constitutive relations for a piezoelectric material [1], the stresses can in two dimensions be converted to the following form
\[ \nabla\cdot\left[ C \left[\begin{array}{c} \nabla u \\ \nabla v\ \\ \nabla V \end{array}\right]\right] = \left[\begin{array}{c} 0 \\ 0 \\ \rho \end{array}\right] \]
where \(C\) is an array of constitutive relations and material parameters. This leads to a system of three equations for the unknown displacements \(u\) and \(v\), and the potential \(V\).
This model is available as an automated tutorial by selecting Model Examples and Tutorials... > Multiphysics > Piezo Electric Bending of a Beam from the File menu. Or alternatively, follow the stepbystep instructions below.
Select the Plane Stress physics mode from the Select Physics dropdown menu.
The beam consists of two 120 x 1 mm rectangular subdomains. First create the upper half and then copy it with a 1 mm offset.
0.12
into the x_{max} edit field.0
into the y_{min} edit field.0.001
into the y_{max} edit field.0 0.001
into the Space separated string of displacement lengths edit field.Press OK to finish and close the dialog box.
As the beam is very thin, manually change the axis ranges to see it better.
Enter 0.005 0.13 0.002 0.002
into the Axis limits edit field.
Press OK to finish and close the dialog box.
Press the Refine Toolbar button.
nu
into the Poisson's ratio edit field.E
into the Modulus of elasticity edit field.Second order accurate finite element shape functions are used for the displacement which are more accurate and allows for a coarser grid.
Select (P2/Q2) second order conforming from the FEM Discretization dropdown menu.
Open the edit equations dialog box and modify the planestress equations to include the coupling terms for the piezoelectric stress contributions, that is e31*Vyx and e31*Vyx.
E/(1nu^2)*( (ux + nu*vy)_x + ((1nu)/2*(uy+vx))_y )  e_31*Vy_x = 0
into the Equation for u edit field.Enter E/(1nu^2)*( (nu*ux + vy)_y + ((1nu)/2*(uy+vx))_x )  e_33*Vy_y = 0
into the Equation for v edit field.
Add the Conductive Media DC physics mode to account for the electrostatics effects.
Select the Conductive Media DC physics mode from the Select Physics dropdown menu.
Again, change and modify the electrostatics equation to account for the piezoelectric effects.
Enter (e0er*Vx_x + (e0er+c_33)*Vy_y + e_31*ux_y + e_33*vy_y) = 0
into the Equation for V edit field.
Name  Expression 

nu  0.29 
E  2e9 
e0er  1.062e10 
d_31  2.2e11 
d_33  4e11 
c_33  (E*d_31*(d_31+nu*d_33))/(nu^21)+(E*d_33*(d_33+nu*d_31))/(nu^21) 
s  1 1 
e_31  s*((E*d_31)/(nu^21)(E*nu*d_33)/(nu^21)) 
e_33  s*((E*d_33)/(nu^21)(E*nu*d_31)/(nu^21)) 
The beam is considered fixed at the left end, which means constrained displacements in all directions.
Select the Fixed displacement, u radio button.
A voltage difference of 250 V is applied between the top and bottom boundaries.
Enter 250
into the Electric potential edit field.
Select Electric potential from the Conductive Media DC dropdown menu.
Plot and visualize the displacement in the ydirection. The maximum displacement it around 29.6 µm which agrees well with the theoretically expected solution.
Select the Deformation Plot check box.
Press OK to plot and visualize the selected postprocessing options.
The piezo electric bending of a beam multiphysics model has now been completed and can be saved as a binary (.fea) model file, or exported as a programmable MATLAB mscript text file, or GUI script (.fes) file.
In this section the model is defined and solved as a MATLAB mscript file. Here we will be looking at composition of two 1.2 cm long and 2 mm thick beams with opposite polarization.
First we assign names to the space coordinates (x and y), and create a grid. Note, that the shape is so simple here that we don't need to explicitly define a geometry, we can use the 'rectgrid' command directly
fea.sdim = { 'x' 'y' }; fea.grid = rectgrid( 20, 20, [ 0 12e3; 1e3 1e3] );
Next is to define the equation coefficients and constitutive relations. See the attached file for their definitions, one thing to note that to define the polarization the coefficient is multiplied with a switch expression 2*(y<0)1 which evaluates to 1 in the top half and 1 in the bottom (here we have also used the name y for the space coordinate in the 2nd direction as we defined earlier).
Now we define the dependent variables and assign 2nd order Lagrange finite element shape functions for them all
fea.dvar = { 'u' 'v' 'V' }; fea.sfun = { 'sflag2' 'sflag2' 'sflag2' };
The material parameters and constitutive relations are defined as
% Equation coefficients. Emod = 2e9; % Modulus of elasticity nu = 0.29; % Poissons ratio Gmod = 0.775e9; % Shear modulus d31 = 0.22e10; % Piezoelectric sn coefficient d33 = 0.3e10; % Piezoelectric sn coefficient prel = 12; % Relative electrical permittivity pvac = 0.88541878176211; % Elec. perm. of vacuum % Constitutive relations. constrel = [ Emod/(1nu^2) nu*Emod/(1nu^2) 0 ; nu*Emod/(1nu^2) Emod/(1nu^2) 0 ; 0 0 Gmod ]; piezoel_st = [ 0 d31 ; 0 d33 ; 0 0 ]; piezoel_sn = constrel*piezoel_st; dielmat_st = [ prel 0 ; 0 prel ]*pvac ; dielmat_sn = [ dielmat_st  piezoel_st'*piezoel_sn ]; % Populate coefficient matrices (negative sign % due to fem partial integration). c{1,1} = { constrel(1,1) constrel(1,3) ; constrel(1,3) constrel(3,3) }; c{1,2} = { constrel(1,3) constrel(1,2) ; constrel(3,3) constrel(2,3) }; c{2,1} = c{1,2}'; c{2,2} = { constrel(3,3) constrel(1,3) ; constrel(1,3) constrel(2,2) }; c{1,3} = { piezoel_sn(1,1) piezoel_sn(1,2) ; piezoel_sn(3,1) piezoel_sn(3,2) }; for i=1:4 c{1,3}{i} = [num2str(c{1,3}{i}),'*(2*(y<0)1)']; end c{3,1} = c{1,3}'; c{2,3} = { piezoel_sn(3,1) piezoel_sn(3,2); piezoel_sn(2,1) piezoel_sn(2,2) }; for i=1:4 c{2,3}{i} = [num2str(c{2,3}{i}),'*(2*(y<0)1)']; end c{3,2} = c{2,3}'; c{3,3} = { dielmat_sn(1,1) dielmat_sn(2,1) ; dielmat_sn(2,1) dielmat_sn(2,2) };
To define the finite element bilinear forms we use the following code
bilinear_form = [ 2 2 3 3 ; 2 3 2 3 ]; for i=1:length(fea.dvar) for j=1:length(fea.dvar) fea.eqn.a.form{i,j} = bilinear_form; fea.eqn.a.coef{i,j} = c{i,j}(:)'; end end
where each bilinear form has four terms, the top line in the definition defines the finite element shape for the dependent variable, and second line for the test function space (a 2 indicates xderivative, and 3 yderivative).
In this case all three source terms are zero, thus
fea.eqn.f.form = { 1 1 1 }; fea.eqn.f.coef = { 0 0 0 };
The boundary conditions are defined as follows
n_bdr = max(fea.grid.b(3,:)); fea.bdr.d = cell(3,n_bdr); fea.bdr.n = cell(3,n_bdr); i_top = 3; i_bottom = 1; i_left = 4; [fea.bdr.d{3,i_top}] = deal( 200 ); [fea.bdr.d{3,i_bottom}] = deal( 0 ); [fea.bdr.d{1:2,i_left}] = deal( 0 );
In this way we have set a potential difference of 200, and fixed the left hand side boundary. Now we parse the problem struct and solve the system with
fea = parseprob( fea ); fea.sol.u = solvestat( fea );
After the problem has been solved we can for example visualize the displacement in the ydirection with the following command
postplot( fea, 'surfexpr', 'v', ... 'isoexpr', 'u', ... 'deformexpr', { 'u', 'v' } )
An example mscript file for this model can be found in the ex_piezoelectric1 mscript file.
[1] Tseng CI. Electromechanical Dynamics of a Coupled Piezoelectric/Mechanical System Applied to Vibration Control and Distributed Sensing. Univ. of Kentucky, Lexington, 1989.
[2] Pierfort V. Finite Element Modeling of Piezoelectric Active Structures. ULB, Faculty of Applied Sciences, 2000.
[3] Hwang WS, Park HC. Finite Element Modeling of Piezoelectric Sensors and Actuators. AIAA Journal, Vol. 31, No. 5, 1993.