|
FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
|
|
FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
EX_PIEZOELECTRIC1 Piezoelectric bending of a beam example.
[ FEA, OUT ] = EX_PIEZOELECTRIC1( VARARGIN ) Example for piezoelectric bending of a beam.
[1] V. Pierfort, Finite Element Modeling of Piezoelectric Active Structures, ULB, Faculty of Applied Sciences, 2000.
[2] W-S. Hwang, H.C. Park, Finite Element Modeling of Piezoelectric Sensors and Actuators, AIAA Journal, Vol. 31, No. 5, May 1993.
[3] C-I. Tseng, Electromechanical Dynamics of a Coupled Piezo-electric/Mechanical System Applied to Vibration Control and Distributed Sensing, Univ. of Kentucky, Lexington, July 1989.
Accepts the following property/value pairs.
Input Value/{Default} Description
-----------------------------------------------------------------------------------
l scalar {0.12} Length of beam
h scalar {2e-3} Height (thickness) of beam
delV scalar {200} Potential difference
igrid scalar 1/{0} Cell type (0=quadrilaterals, 1=triangles)
hmax scalar {20} Cell resolution
sfund string {sflag1} Shape function for displacements
sfunp string {sflag1} Shape function for potential
iplot scalar 0/{1} Plot solution (=1)
.
Output Value/(Size) Description
-----------------------------------------------------------------------------------
fea struct Problem definition struct
out struct Output struct
cOptDef = { ...
'l', 0.12; ...
'h', 2e-3; ...
'delV', 100; ...
'igrid', 0; ...
'hmax', 20; ...
'sfund', 'sflag2'; ...
'sfunp', 'sflag1'; ...
'iplot', 1; ...
'fid', 1 };
[got,opt] = parseopt( cOptDef, varargin{:} );
fid = opt.fid;
i_impl = 1;
% Geometry definition.
fea.sdim = { 'x' 'y' }; % Names of space coordinates.
gobj1 = gobj_rectangle( 0, opt.l, 0, opt.h/2, 'R1' );
gobj2 = gobj_rectangle( 0, opt.l, -opt.h/2, 0, 'R2' );
fea.geom.objects = { gobj1 gobj2 };
% Grid generation.
switch opt.igrid
case -1
fea.grid = rectgrid( opt.hmax, opt.hmax, [ 0 opt.l; -opt.h/2 opt.h/2] );
fea.grid = quad2tri( fea.grid );
case 0
fea.grid = rectgrid( opt.hmax, opt.hmax, [ 0 opt.l; -opt.h/2 opt.h/2] );
fea.grid.s( selcells( fea, 'y<=0' ) ) = 2;
case 1
fea.grid = gridgen( fea, 'hmax', opt.h/(opt.hmax/4.6), 'fid', fid );
fea.grid.s(:) = 1;
end
ind_c = selcells( fea, 'y<=0' );
fea.grid.s( ind_c ) = 2;
ix = ismember( fea.grid.b(1,:), ind_c );
fea.grid.b(3,ix) = fea.grid.b(3,ix) + 4;
[~,~,ib] = unique(fea.grid.b(3,:));
fea.grid.b(3,:) = ib;
% 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.885418781762e-11; % Electrical permittivity 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) };
if( i_impl )
for i=1:4
c{1,3}{i} = [num2str(c{1,3}{i}),'*(2*(y<0)-1)'];
end
else
for i=1:4
fea.expr{i,1} = ['piezoel_sn1',num2str(i)];
fea.expr{i,2} = { -c{1,3}{i} c{1,3}{i} };
c{1,3}{i} = ['piezoel_sn1',num2str(i)];
end
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) };
if( i_impl )
for i=1:4
c{2,3}{i} = [num2str(c{2,3}{i}),'*(2*(y<0)-1)'];
end
else
for i=1:4
fea.expr{i+4,1} = ['piezoel_sn2',num2str(i)];
fea.expr{i+4,2} = { -c{2,3}{i} c{2,3}{i} };
c{2,3}{i} = ['piezoel_sn2',num2str(i)];
end
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) };
% Dependent variable names.
fea.dvar = { 'u' 'v' 'V' };
n_dvar = length(fea.dvar);
% Finite element shape functions.
fea.sfun = { opt.sfund opt.sfund opt.sfunp };
% Define equations.
bilinear_form = [ 2 2 3 3 ;
2 3 2 3 ];
for i=1:n_dvar
for j=1:n_dvar
fea.eqn.a.form{i,j} = bilinear_form;
fea.eqn.a.coef{i,j} = c{i,j}(:)';
end
end
% Source terms (set to zero).
fea.eqn.f.form = { 1 1 1 };
fea.eqn.f.coef = { 0 0 0 };
% Boundary conditions.
n_bdr = max(fea.grid.b(3,:)); % Number of boundaries.
fea.bdr.d = cell(n_dvar,n_bdr);
fea.bdr.n = cell(n_dvar,n_bdr);
i_top = findbdr( fea, ['y>=',num2str(opt.h/2-sqrt(eps))] );
i_bottom = findbdr( fea, ['y<=',num2str(-opt.h/2+sqrt(eps))] );
i_left = findbdr( fea, ['x<=',num2str(sqrt(eps))] );
[fea.bdr.d{3,i_top}] = deal( opt.delV ); % Set potential to dV on top boundary.
[fea.bdr.d{3,i_bottom}] = deal( 0 ); % Set potential to 0V on bottom boundary.
[fea.bdr.d{1:2,i_left}] = deal( 0 ); % Set displacements to 0 on left boundary.
% Check and parse problem struct.
fea = parseprob( fea );
% Call to stationary solver.
fea.sol.u = solvestat( fea, 'fid', fid );
% Postprocessing.
if( opt.iplot )
YSCALE = 3;
axlim = [0, opt.l, -YSCALE*opt.h/2, YSCALE*opt.h/2];
DSCALE = 20;
subplot(2,2,1)
postplot( fea, 'surfexpr', 'u', 'isoexpr', 'u', 'setaxes', 'off' )
axis(axlim);
title('x-displacement')
subplot(2,2,2)
postplot( fea, 'surfexpr', 'v', 'isoexpr', 'v', 'setaxes', 'off' )
axis(axlim);
title('y-displacement')
subplot(2,2,3)
postplot( fea, 'surfexpr', 'V', 'isoexpr', 'V', 'setaxes', 'off' )
axis(axlim);
title('Electric potential')
subplot(2,2,4)
plotsubd( fea, 'labels', 'off', 'setaxes', 'off' )
axis(axlim);
title(['Displacement plot (at ',num2str(DSCALE),' times scale)'])
ind1 = sub2ind( size(fea.grid.c), fea.grid.b(2,:), fea.grid.b(1,:) );
p1b = fea.grid.p( :, fea.grid.c(ind1) );
ind2 = sub2ind( size(fea.grid.c), mod(fea.grid.b(2,:),size(fea.grid.c,1))+1, fea.grid.b(1,:) );
p2b = fea.grid.p( :, fea.grid.c(ind2) );
up1 = DSCALE*evalexpr( 'u', p1b, fea );
vp1 = DSCALE*evalexpr( 'v', p1b, fea );
up2 = DSCALE*evalexpr( 'u', p2b, fea );
vp2 = DSCALE*evalexpr( 'v', p2b, fea );
hold on
for i=1:size(p1b,2)
plot( [p1b(1,i)+up1(i) p2b(1,i)+up2(i)], [p1b(2,i)+vp1(i) p2b(2,i)+vp2(i)], '-b', 'linewidth', 2 )
end
set( gca, 'ytick', [] )
end
% Error checking.
ind_dof_v = fea.eqn.dofm{2}(:) + fea.eqn.ndof(1);
out.v_max = max(abs( fea.sol.u(ind_dof_v) ));
v_max_ref = abs(-3/2*d31*opt.delV*(opt.l/opt.h)^2);
out.err = abs(v_max_ref - out.v_max)/v_max_ref;
out.pass = out.err <= 0.1;
if ( nargout==0 )
clear fea out
end