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

Этот пример показывает, как решить связанную задачу эластичности-электростатики.

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

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

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

$- \nabla \cdot σ = f$

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

$\nabla \cdot D = ρ$

$D$ является электрическим перемещением, и $ρ$ распространяется бесплатно. Объедините эти две системы PDE в одну систему:

$- \nabla \cdot {\begin{array}{c} σ \\ D \end{array}} = {\begin{array}{c} f \\ - ρ \end{array}}$

Для анализа 2-D, $σ$ имеет компоненты $σ_{11}, σ_{22},$ и $σ_{12} = σ_{21}$ , и $D$ имеет компоненты $D_{1}$ и $D_{2}$ .

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

${\begin{array}{c} σ_{11} \\ σ_{22} \\ σ_{12} \\ D_{1} \\ D_{2} \end{array}} = [\begin{array}{ccccc} C_{11} & C_{12} & e_{11} & e_{31} \\ C_{12} & C_{22} & e_{13} & e_{33} \\ G_{12} & e_{14} & e_{34} \\ e_{11} & e_{13} & e_{14} & - E_{1} \\ e_{31} & e_{33} & e_{34} & - E_{2} \end{array}] {\begin{array}{c} ϵ_{11} \\ ϵ_{22} \\ γ_{12} \\ - E_{1} \\ - E_{2} \end{array}}$

$C_{i j}$ являются упругими коэффициентами, $E_{i}$ являются электрическими проницаемостями, и $e_{i j}$ являются пьезоэлектрическими коэффициентами напряжения. Пьезоэлектрические коэффициенты напряжения в этих уравнениях соответствуют обычному обозначению в пьезоэлектрических материалах, где z-направление (третье направление) совпадает с «полюсным» направлением материала. Для анализа 2-D выровняйте «полярное» направление по оси Y. Напишите вектор деформации в терминах перемещения X $u$ и y-перемещение $v$ :

${\begin{array}{c} ϵ_{11} \\ ϵ_{22} \\ γ_{12} \end{array}} = {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{array}}$

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

${\begin{array}{c} E_{1} \\ E_{2} \end{array}} = - {\begin{array}{c} \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

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

Partial Differential Equation Toolbox™ требует, чтобы система эллиптических уравнений была выражена в вектор форме:

$- \nabla \cdot (c \otimes \nabla u) + a u = f$

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

$- \frac{\partial}{\partial x_{k}} (c_{i j k l} \frac{\partial u_{j}}{\partial x_{l}}) + a_{i j} u_{j} = f_{i}$

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

$u = {\begin{array}{c} u \\ v \\ ϕ \end{array}}$

Это $N = 3$ система. Градиент $u$ является

$\nabla u = {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

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

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

$[\begin{array}{cccccc} c (1) & c (2) & c (4) & c (6) & c (11) & c (13) \\ \cdot & c (3) & c (5) & c (7) & c (12) & c (14) \\ \cdot & \cdot & c (8) & c (9) & c (15) & c (17) \\ \cdot & \cdot & \cdot & c (10) & c (16) & c (18) \\ \cdot & \cdot & \cdot & \cdot & c (19) & c (20) \\ \cdot & \cdot & \cdot & \cdot & \cdot & c (21) \end{array}]$

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

$[\begin{array}{cccccc} c_{1111} & c_{1112} & c_{1211} & c_{1212} & c_{1311} & c_{1312} \\ \cdot & c_{1122} & c_{1221} & c_{1222} & c_{1321} & c_{1322} \\ \cdot & \cdot & c_{2211} & c_{2212} & c_{2311} & c_{2312} \\ \cdot & \cdot & \cdot & c_{2222} & c_{2321} & c_{2322} \\ \cdot & \cdot & \cdot & \cdot & c_{3311} & c_{3312} \\ \cdot & \cdot & \cdot & \cdot & \cdot & c_{3322} \end{array}] {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

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

$[\begin{array}{cccccc} C_{11} & \cdot & \cdot & C_{12} & e_{11} & e_{31} \\ \cdot & G_{12} & G_{12} & \cdot & e_{14} & e_{34} \\ \cdot & \cdot & G_{12} & \cdot & e_{14} & e_{34} \\ \cdot & \cdot & \cdot & C_{22} & e_{13} & e_{33} \\ \cdot & \cdot & \cdot & \cdot & - E_{1} & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & - E_{2} \end{array}] {\begin{array}{c} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial ϕ}{\partial x} \\ \frac{\partial ϕ}{\partial y} \end{array}}$

Геометрия балки

Создайте модель УЧП. Уравнения линейной упругости имеют три компонента, поэтому модель должна иметь три уравнения.

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

Figure contains an axes. The axes contains 10 objects of type line, text.

Свойства материала

Задайте свойства материала слоев балки. Материал в обоих слоях представляет собой поливинилиденфторид (PVDF), термопластичный полимер с пьезоэлектрическим поведением.

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;

Можно просмотреть 21 коэффициент как матрицу 3 на 3 блоков 2 на 2. Матрицы cij являются блоками 2 на 2 в верхнем треугольнике этой матрицы.

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

Граничные условия

Установите напряжение (компонент решения 3) в верхней части луча (ребро 1) на 100 вольт.

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

Задайте, что нижняя часть балки (ребро 2) заземлена, установив напряжение равным 0.

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

Укажите, что левая сторона (ребра 6 и 7) зажата, установив смещения X и Y (решения раствора 1 и 2) равными 0.

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

Напряжение и заряд на правой стороне балки равны нулю. Соответственно, используйте граничное условие по умолчанию для ребер 3 и 4.

Конечные элементные и аналитические решения

Сгенерируйте mesh и решите модель.

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]

Доступ к решению в узловых местоположениях. Первый столбец содержит x-отклонение. Второй столбец содержит отклонение y. Третий столбец содержит электрический потенциал.

rs = result.NodalSolution;

Найдите минимальное отклонение Y.

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

Figure contains an axes. The axes with title X-Deflection, meters contains 12 objects of type patch, line.

Figure contains an axes. The axes with title Y-Deflection, meters contains 12 objects of type patch, line.

Figure contains an axes. The axes with title Electrical Potential, Volts contains 12 objects of type patch, line.

Ссылки

Хван, У-Сок и Хён Чул Парк. Конечный элемент Моделирования пьезоэлектрических датчиков и Приводы. Журнал AIAA 31, № 5 (май 1993): 930-937.
Пьфорд, В. «Конечный Элемент моделирование пьезоэлектрических активных структур». PhD diss., Universite Libre De Bruxelles, 2001.

Документация