SDT-piezo Contents   Functions      PDF Index |
In this very simple example, the electric field and the strains are all constant, so that the electric potential and the displacement field are linear. The example is treated analytically in section section 2.4. It is therefore possible to obtain an exact solution using a single volumic 8-node finite element (with linear shape functions, the nodal unknowns being the displacements in x,y and z and the electric potential φ). Consider a piezoelectric patch whose dimensions and material properties are given in Table 3.1. The material properties correspond to the material SONOX_P502_iso in m_piezo.
Property Value b 10 mm w 10 mm h 2 mm E 54 GPa ν 0.41 d31=d32 -185 10−12 pC/N (or m/V) d33 440 10−12 pC/N (or m/V) d15=d24 560 10−12 pC/N (or m/V) є33T=є22T=є11T 1850 є0 є0 8.854 10−12 Fm−1
We first produce the mesh, associate the material properties and define the electrodes with
d_piezo('TutoPatch-s1')
. The default material is SONOX_P502_iso. The number of elements in the x, y and z directions are given by nx,ny and nz.
% See full example as MATLAB code in d_piezo('ScriptTutoPatch') t_avc('DefineStyles'); %% Step 1 Build mesh - Define electrodes % Meshing script can be viewed with sdtweb d_piezo('MeshPatch') model=d_piezo('MeshPatch lx=1e-2 ly=1e-2 h=2e-3 nx=1 ny=1 nz=1'); % Define electrodes model=p_piezo('ElectrodeMPC Top -ground',model,'z==2e-3'); model=p_piezo('ElectrodeMPC Bottom -Input "Free patch"',model,'z==0');
The information about the nodes associated to each electrode can be obtained through the following call:
p_piezo('TabInfo',model)
The material can be changed for example to PIC_255 with the following call, and the full set of mechanical, piezoelectric and permittivity matrices can be obtained in order to check consistency with the datasheet (d_piezo('TutoPatch-s2') ):
%% Step 2 Define material properties model.pl=m_piezo('dbval 1 -elas 2 PIC_255'); p_piezo('TabDD',model) % creates the table with full set of matrices
The next step consists in defining the boundary conditions and load case using
d_piezo('TutoPatch-s3')
. We consider here two cases, the first one where the patch is free to expand, and the second one where it is mechanically constrained (all mechanical degrees of freedom are equal to 0).
%% Step 3 Compute static response % to avoid rigid body mode model=stack_set(model,'info','Freq',10); def=fe_simul('dfrf',model); def.lab={'Free patch, axial'}; def.fun=[0 1]; def=feutil('rmfield',def,'data','LabFcn'); % Append mechanically constrained structure % can't call fe_simul because no free DOF % see code with sdtweb d_piezo('scriptFullConstrain') def=d_piezo('scriptFullConstrain',model,def); def.lab{2}='Constrained patch, axial';
We can look at the deformed shape, and plot the electric field for both cases.
(d_piezo('TutoPatch-s4')
%% Step 4 Visualize deformed shape cf=feplot(model,def); % Electric field representation p_piezo('viewElec EltSel "matid1" DefLen 20e-4 reset',cf); fecom('colormap',[1 0 0]);fecom('undef line');iimouse('resetview'); cf.mdl.name='E-field'; % Sets figure name t_avc('SetStyle',cf); feplot(cf);
For the free patch deformed shape, we compute the mean strains from which d31, d32 and d33 are deduced. The values are found to be equal to the analytical values used in the model. Note that the parameters of the constitutive equations can be recovered using (d_piezo('TutoPatch-s5') ):
%% Step 5 : check constitutive law % Decompose constitutive law CC=p_piezo('viewdd -struct',cf); %
where the fields of CC are self-explanatory. The parameters which are not directly defined are computed from the equations presented in Section section 2.1.
% Display and compute mean strains a=p_piezo('viewstrain -curve -mean',cf); % Strain S fprintf('Relation between mean strain on free structure and d_3i\n'); E3=a.Y(9,1); disp({'E3 mean' a.Y(9,1) 1/2e-3 'E3 analytic'}) disp([a.X{1}(1:3) num2cell([a.Y(1:3,1)/E3 CC.d(3,1:3)']) ... {'d_31';'d_32';'d_33'}])
For the constrained patch, we compute the mean stress from which we can compute the e31, e32 and e33 values which are found to be equal to the analytical values used in the model:
% Display and compute mean stresses b=p_piezo('viewstress -curve -mean',cf); % Stress T fprintf('Relation between mean stress on pure electric and e_3i\n'); disp([b.X{1}(1:3) num2cell([b.Y(1:3,2)/-E3 CC.e(3,1:3)']) ... {'e_31';'e_32';'e_33'}]) % Mean stress/strain disp([b.X{1} num2cell(b.Y(:,2)) num2cell(a.Y(:,1)) a.X{1}])
We can also compute the charge and the charge density (in pC/m2) accumulated on the electrodes, and compare with the analytical values (d_piezo('TutoPatch-s6') ):
%% Step 6 Check capacitance values % Theoretical values of Capacitance and charge density - free patch CT=CC.epst_r(3,3)*8.854e-12*1e-2*1e-2/2e-3; %% Capacitance - free patch CdensT=CC.epst_r(3,3)*8.854e-12/2e-3*1e12; %% charge density - free patch % Theoretical values of Capacitance and charge density - constrained patch CS=CC.epss_r(3,3)*8.854e-12*1e-2*1e-2/2e-3; %% Capacitance - free patch CdensS=CC.epss_r(3,3)*8.854e-12/2e-3*1e12; %% charge density - free patch % Represent charge density (C/S) value on the electrodes % - compare with analytical values cut=p_piezo('electrodeviewcharge',cf,struct('EltSel','matid 1')); b=fe_caseg('stressobserve',cut,cf.def);b=reshape(b.Y,[],2); disp([{'','CdensT','CdensS'};{'Numeric';'Theoretical'} ... num2cell([mean(abs(b));CdensT CdensS])]) iimouse('zoom reset'); % Compute the value of the total charge (from reaction at electrical dof) % Compare with analytical values p_piezo('electrodeTotal',cf) % disp('Theoretical values of capacitance') disp([{'CT';'CS'} num2cell([CT;CS])]) cf.mdl.name='Charge'; % Sets figure name t_avc('SetStyle',cf); feplot(cf); fecom('ch 1')
The results clearly show the very large difference of charge density between the two cases (free patch or constrained patch).
For this simple static example, a finer mesh can be used, but it does not lead to more accurate results (this can be done by changing the values in the call to d_piezo('mesh') for example:
% Build mesh with refinement model=d_piezo('MeshPatch lx=1e-2 ly=1e-2 h=2e-3 nx=5 ny=5 nz=2'); % Now a model with quadratic elements model=d_piezo('MeshPatch lx=1e-2 ly=1e-2 h=2e-3 Quad');
As for the patch in extension, as the fields are also uniform (see section 2.5), the problem can be modelled with a single 8-node element. The patch is meshed and then the poling is aligned with the −y axis by performing a rotation of 90o around the x-axis (d_piezo('TutoPatchShear-s1') ).
t_avc('DefineStyles'); % See full example as MATLAB code in d_piezo('ScriptTutoShearPatch') %% Step 1 Build mesh and define electrodes %Meshing script can be viewed with sdtweb d_piezo('MeshPatch') model=d_piezo('MeshPatch lx=1e-2 ly=1e-2 h=2e-3 nx=1 ny=1 nz=1'); % Define electrodes model=p_piezo('ElectrodeMPC Top -ground',model,'z==2e-3'); model=p_piezo('ElectrodeMPC Bottom -Input "Free patch"',model,'z==0'); % Rotate basis to align poling direction with y (-90 deg around x) model.bas=basis('rotate',[],'rx=-90',1); %create local basis with id=1 model=feutil('setpro 1 COORDM=1',model); % assign basis with id=1 to pro=1
Then the response is computed both for the free case and the fully constrained case. The deformed shape for the free case is shown in Figure 3.6 together with the applied electric field:
(d_piezo('TutoShearPatch-s2')
)
%% Step 2 Compute static response % to avoid rigid body mode model=stack_set(model,'info','Freq',10); def=fe_simul('dfrf',model); def.lab={'Free patch, shear'}; def.fun=[0 1]; def=feutil('rmfield',def,'data','LabFcn'); % Append mechanically constrained structure % can't call fe_simul because no free DOF % see code with sdtweb d_piezo('scriptFullConstrain') def=d_piezo('scriptFullConstrain',model,def); def.lab{2}='Constrained patch, shear';
(d_piezo('TutoShearPatch-s3') )
%% Step 3 Vizualise deformed shape cf=feplot(model,def); fecom('undef line'); % Electric field representation p_piezo('viewElec EltSel "matid1" DefLen 20e-4 reset',cf); cf.mdl.name='E-field'; % Sets figure name t_avc('SetStyle',cf); feplot(cf);
The mean of shear strain and stress is evaluated and compared to the d24 piezo coefficient.
Note that the mean values are computed in the global yz axis for which a negative strain corresponds to a positive strain in the local 23 axis.
Finally the capacitance is evaluated and compared to the theoretical values, showing a perfect agreement, and demonstrating the difference with the extension case for CS.
(d_piezo('TutoShearPatch-s4')
)
%% Step 4 : Check constitutive law % Decompose constitutive law CC=p_piezo('viewdd -struct',cf); % % Display and compute mean strains a=p_piezo('viewstrain -curve -mean',cf); % Strain S fprintf('Relation between mean strain on free structure and d_24\n'); E3=a.Y(9,1); disp({'E3 mean' a.Y(9,1) 1/2e-3 'E3 analytic'}) disp([a.X{1}(4) num2cell([a.Y(4,1)/E3 CC.d(2,4)']) ... {'d_24'}]) % Display and compute mean stresses b=p_piezo('viewstress -curve -mean',cf); % Stress T fprintf('Relation between mean stress on pure electric and e_24 \n'); disp([b.X{1}(4) num2cell([b.Y(4,2)/-E3 CC.e(2,4)']) ... {'e_24'}])
% Mean stress/strain disp([b.X{1} num2cell(b.Y(:,2)) num2cell(a.Y(:,1)) a.X{1}]) % Theoretical values of Capacitance and charge density - free patch CT=CC.epst_r(2,2)*8.854e-12*1e-2*1e-2/2e-3; %% Capacitance - free patch CdensT=CC.epst_r(2,2)*8.854e-12/2e-3*1e12; %% charge density - free patch % Theoretical values of Capacitance and charge density - constrained patch CS=CC.epss_r(2,2)*8.854e-12*1e-2*1e-2/2e-3; %% Capacitance - constrained patch CdensS=CC.epss_r(2,2)*8.854e-12/2e-3*1e12; %% charge density - constrained patch % Represent charge density (C/S) value on the electrodes % - compare with analytical values cut=p_piezo('electrodeviewcharge',cf,struct('EltSel','matid 1')); b=fe_caseg('stressobserve',cut,cf.def);b=reshape(b.Y,[],2); disp([{'','CdensT','CdensS'};{'Numeric';'Theoretical'} ... num2cell([mean(abs(b));CdensT CdensS])]) iimouse('zoom reset');
(d_piezo('TutoShearPatch-s5') )
%% Step 5 Check capacitance % Compute the value of the total charge (from reaction at electrical dof) % Ccompare with analytical values p_piezo('electrodeTotal',cf) % disp('Theoretical values of capacitance') disp([{'CT';'CS'} num2cell([CT;CS])])