Отклонение пьезоэлектрического привода

В этом примере показано, как решить проблему связанной упругости и электростатики.

Пьезоэлектрические материалы деформируются под приложенным напряжением. Наоборот, деформация пьезоэлектрического материала создает напряжение. Поэтому анализ пьезоэлектрической части требует решения набора связанных дифференциальных уравнений в частных производных с отклонениями и электрическим потенциалом в качестве зависимых переменных.

В этом примере модель представляет собой двухслойную консольную балку с обоими слоями, выполненными из одного и того же материала поливинилиденфторида (PVDF). Направление поляризации указывает вниз (отрицательное направление y) в верхнем слое и вверх в нижнем слое. Типичное отношение длины к толщине составляет 100. При подаче напряжения между нижней и верхней поверхностями балки она отклоняется в направлении y, поскольку один слой укорачивается, а другой - удлиняется.

Уравнения равновесия описывают упругое поведение твердого тела:

В данном случае $$\lambda$$ - тензор напряжения, а $$f$$ - вектор силы тела. Закон Гаусса описывает электростатическое поведение твердого вещества:

$$D$$ является электрическим смещением, а $$\mathrm{(«»)}$$ - распределенным бесплатным зарядом. Объедините эти две системы PDE в одну систему:

Для 2-го анализа $$\mathrm{у\; \sigma }$$ есть компоненты $${}_{\mathrm{}\mathrm{\sigma 11}},{}_{\mathrm{}\mathrm{\sigma 22}},$$ и $${}_{\mathrm{}\mathrm{\sigma 12}}{}_{\mathrm{}\mathrm{=\sigma 21}}$$, и $$\mathrm{у\; D}$$ есть компоненты $${}_{\mathrm{D1}}$$ и $${}_{\mathrm{D2}}$$.

Конститутивные уравнения для материала определяют тензор напряжения и вектор электрического смещения в терминах тензора деформации и электрического поля. Для 2-D анализа ортотропного пьезоэлектрического материала в условиях плоского напряжения можно записать эти уравнения как

$${}_{\mathrm{Cij}}$$ - коэффициенты упругости, $${\mathcal{}}_{\mathrm{Ei}}$$ - электрические проницаемости, и $${}_{\mathrm{eij}}$$ - коэффициенты пьезоэлектрических напряжений. Коэффициенты пьезоэлектрических напряжений в этих уравнениях соответствуют обычным обозначениям в пьезоэлектрических материалах, где направление z (третье направление) выравнивается с «полированным» направлением материала. Для анализа 2-D выровняйте направление «полюса» по оси Y. Запишите вектор деформации в терминах x-смещения $$u$$ и y-смещения $$v$$:

Запишите электрическое поле с точки зрения электрического потенциала$$$$

Можно подставить уравнения деформаций-смещений и уравнения электрического поля в конститутивные уравнения и получить систему уравнений для напряжений и электрических смещений в терминах смещения и производных электрического потенциала. Подстановка результирующих уравнений в уравнения системы PDE дает систему уравнений, которые включают в себя расхождение производных смещения и электрического потенциала. На следующем шаге расположите эти уравнения в соответствии с формой, требуемой для панели инструментов.

Дифференциальное уравнение в частных производных Toolbox™ требует, чтобы система эллиптических уравнений была выражена в векторной форме:

или в тензорной форме:

где повторные индексы подразумевают суммирование. Для 2-й пьезоэлектрической системы в этом примере системный вектор $$u$$

Это система $$N\mathrm{=}$$ 3. Градиент $$$$u равен

Дополнительные сведения об указании коэффициентов в формате, требуемом панелью инструментов, см. в разделе:

Коэффициент c в этом примере является тензором. Ее можно представить в виде матрицы 3 на 3 из блоков 2 на 2:

Для отображения членов конститутивных уравнений в форму, требуемую панелью инструментов, запишите тензор c и градиент решения в этой форме:

Из этого уравнения можно сопоставить традиционные конститутивные коэффициенты с формой, требуемой для матрицы c. Знак минус в уравнениях для электрического поля включается в матрицу c для соответствия условию панели инструментов.

model = createpde(3);

L = 100e-3; % Beam length in meters
H = 1e-3; % Overall height of the beam
H2 = H/2; % Height of each layer in meters

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']');

geometryFromEdges(model,g);

figure
pdegplot(model,'EdgeLabels','on','FaceLabels','on')
xlabel('X-coordinate, meters')
ylabel('Y-coordinate, meters')
axis([-.1*L,1.1*L,-4*H2,4*H2])
axis square

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;

relPermittivity = 12;

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];

pzeE = c2d*pzeD;
D_const_stress = [relPermittivity 0; 0 relPermittivity]*permittivityFreeSpace;

D_const_strain = D_const_stress - pzeD'*pzeE;

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)];
ctop = [c11(:); c12(:); c22(:); -c13(:); -c23(:); -c33(:)];
cbot = [c11(:); c12(:); c22(:);  c13(:);  c23(:); -c33(:)];

f = [0 0 0]';
specifyCoefficients(model,'m',0,'d',0,'c',ctop,'a',0,'f',f,'Face',2);
specifyCoefficients(model,'m',0,'d',0,'c',cbot,'a',0,'f',f,'Face',1);

voltTop = applyBoundaryCondition(model,'mixed','Edge',1,...
                                               'u',100,...
                                               'EquationIndex',3);

voltBot = applyBoundaryCondition(model,'mixed','Edge',2,...
                                               'u',0,...
                                               'EquationIndex',3);

clampLeft = applyBoundaryCondition(model,'mixed','Edge',6:7,...
                                                 'u',[0 0],...
                                                 'EquationIndex',1:2);

msh = generateMesh(model,'Hmax',5e-4);
result = solvepde(model)

result = 
  StationaryResults with properties:

    NodalSolution: [3605x3 double]
       XGradients: [3605x3 double]
       YGradients: [3605x3 double]
       ZGradients: [0x3 double]
             Mesh: [1x1 FEMesh]

rs = result.NodalSolution;

feTipDeflection = min(rs(:,2));
fprintf('Finite element tip deflection is: %12.4e\n',feTipDeflection);

Finite element tip deflection is:  -3.2900e-05

tipDeflection = -3*d31*100*L^2/(8*H2^2);
fprintf('Analytical tip deflection is: %12.4e\n',tipDeflection);

Analytical tip deflection is:  -3.3000e-05

varsToPlot = char('X-Deflection, meters', ...
                  'Y-Deflection, meters', ...
                  'Electrical Potential, Volts');
for i = 1:size(varsToPlot,1)
  figure;
  pdeplot(model,'XYData',rs(:,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')
  axis square
end