FEATool Multiphysics
v1.10 Finite Element Analysis Toolbox |

Piezoelectric Bending of a Beam

This example describes how to model piezoelectric bending of a beam on the command line as an m-script 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 step-by-step instructions below.

- To start a new model click the
**New Model**toolbar button, or select*New Model...*from the*File*menu. - Select the
**2D**radio button. Select the

**Plane Stress**physics mode from the*Select Physics*drop-down menu.- Press
**OK**to finish the physics mode selection.

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.

- Select
**Rectangle**from the*Geometry*menu. - Enter
`0.12`

into the*x*edit field._{max} - Enter
`0`

into the*y*edit field._{min} - Enter
`0.001`

into the*y*edit field._{max} - Press
**OK**to finish and close the dialog box. - Select
**R1**in the geometry object*Selection*list box. - Press the
**Copy and/or transform selected geometry object***Toolbar*button. - Select the
**Make copy of geometry object**check box. - Enter
`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.

- Select
**Axis/Grid Settings...**from the*Options*menu. - Select the
**Axis limits**radio button. Enter

`-0.005 0.13 -0.002 0.002`

into the*Axis limits*edit field.Press

**OK**to finish and close the dialog box.- Switch to
**Grid**mode by clicking on the corresponding*Mode Toolbar*button. Press the

**Refine***Toolbar*button.- Switch to
**Equation**mode by clicking on the corresponding*Mode Toolbar*button. - Enter
`nu`

into the*Poisson's ratio*edit field. - Enter
`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*drop-down menu.- Press the
**edit**button.

Open the edit equations dialog box and modify the plane-stress equations to include the coupling terms for the piezo-electric stress contributions, that is *e31*Vyx* and *e31*Vyx*.

- Enter
`-E/(1-nu^2)*( (ux + nu*vy)_x + ((1-nu)/2*(uy+vx))_y ) - e_31*Vy_x = 0`

into the*Equation for u*edit field. Enter

`-E/(1-nu^2)*( (nu*ux + vy)_y + ((1-nu)/2*(uy+vx))_x ) - e_33*Vy_y = 0`

into the*Equation for v*edit field.- Press
**OK**to finish and close the dialog box.

Add the *Conductive Media DC* physics mode to account for the electrostatics effects.

- Switch to the
**+**tab. Select the

**Conductive Media DC**physics mode from the*Select Physics*drop-down menu.- Press the
**Add Physics >>>**button. - Press the
**edit**button.

Again, change and modify the electrostatics equation to account for the piezo-electric 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.- Press
**OK**to finish equation editing and close the subdomain settings dialog box. - In order to define the piezo-electric modeling constants and expressions, press the
**Constants***Toolbar*button, or select the corresponding entry from the*Equation*menu, and enter the following variables in the*Model Constants and Expressions*dialog box. Note how the scaling factor*s =-1 1*is used to invert the coefficients between the upper and lower halves. See the reference for more details.

Name | Expression |
---|---|

nu | 0.29 |

E | 2e9 |

e0er | 1.062e-10 |

d_31 | 2.2e-11 |

d_33 | -4e-11 |

c_33 | (E*d_31*(d_31+nu*d_33))/(nu^2-1)+(E*d_33*(d_33+nu*d_31))/(nu^2-1) |

s | -1 1 |

e_31 | s*(-(E*d_31)/(nu^2-1)-(E*nu*d_33)/(nu^2-1)) |

e_33 | s*(-(E*d_33)/(nu^2-1)-(E*nu*d_31)/(nu^2-1)) |

- Switch to
**Boundary**mode by clicking on the corresponding*Mode Toolbar*button.

The beam is considered fixed at the left end, which means constrained displacements in all directions.

- Select
**3**and**6**in the*Boundaries*list box. - Select the
**Fixed displacement, v**radio button. Select the

**Fixed displacement, u**radio button.- Switch to the
**dc**tab. - Select
**2**in the*Boundaries*list box.

A voltage difference of *250 V* is applied between the top and bottom boundaries.

- Select
**Electric potential**from the*Conductive Media DC*drop-down menu. Enter

`250`

into the*Electric potential*edit field.- Select
**4**in the*Boundaries*list box. Select

**Electric potential**from the*Conductive Media DC*drop-down menu.- Press
**OK**to finish the boundary condition specification. - Switch to
**Solve**mode by clicking on the corresponding*Mode Toolbar*button. - Press the
**=***Toolbar*button to call the solver. After the problem has been solved FEATool will automatically switch to postprocessing mode and plot the computed solution.

Plot and visualize the displacement in the y-direction. The maximum displacement it around *29.6 µm* which agrees well with the theoretically expected solution.

- Press the
**Plot Options***Toolbar*button. - Select
**y-displacement**from the*Predefined surface plot expressions*drop-down menu. - Select the
**Contour Plot**check box. - Select
**x-displacement**from the*Predefined contour plot expressions*drop-down menu. 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 m-script text file, or GUI script (.fes) file.

In this section the model is defined and solved as a MATLAB m-script 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 12e-3; -1e-3 1e-3] );

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.22e-10; % Piezoelectric sn coefficient d33 = -0.3e-10; % Piezoelectric sn coefficient prel = 12; % Relative electrical permittivity pvac = 0.885418781762-11; % Elec. perm. of vacuum % Constitutive relations. constrel = [ Emod/(1-nu^2) nu*Emod/(1-nu^2) 0 ; nu*Emod/(1-nu^2) Emod/(1-nu^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 x-derivative, and 3 y-derivative).

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 y-direction with the following command

postplot( fea, 'surfexpr', 'v', ... 'isoexpr', 'u', ... 'deformexpr', { 'u', 'v' } )

An example m-script file for this model can be found in the ex_piezoelectric1 m-script file.

[1] Tseng CI. *Electromechanical Dynamics of a Coupled Piezo-electric/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.