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

В этом примере показано, как решить двойную задачу электростатики эластичности.

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

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

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

$- \nabla \cdot σ = f$

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

$\nabla \cdot D = ρ$

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

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

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

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

${\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-направление (третье направление) выравнивается с "подпертым шестами" направлением материала. Для 2D анализа выровняйте "подпертое шестами" направление с осью 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}}$

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

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}$

где повторные индексы подразумевают суммирование. Для 2D пьезоэлектрической системы в этом примере, системном векторе $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}}$

Излучите геометрию

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

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 object. The axes object contains 10 objects of type line, text.

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

Задайте свойства материала слоев луча. Материал в обоих слоях является polyvinylidene фторидом (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 object. The axes object with title X-Deflection, meters contains 12 objects of type patch, line.

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

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

Ссылки

Хван, Добейтесь-Seok, и парк Hyun Chul. "Моделирование Конечного элемента Пьезоэлектрических Датчиков и Приводов". Журнал 31, № 5 AIAA (май 1993): 930-937.
Пифорд, V. "Моделирования конечного элемента Пьезоэлектрических Активных Структур". Доктор философии diss., Брюссельский свободный университет, 2001.

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