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

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

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

УЧП для пьезоэлектрического тела

Эластичное поведение тела описано уравнениями равновесия

$- \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-направление (с 3 направлениями), выравнивается с "подпертым шестами" направлением материала. Для 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}}$

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

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

Преобразование уравнений к форме тулбокса УЧП

Тулбокс УЧП требует, чтобы система эллиптических уравнений была выражена в форме

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

Документация для функционального assempde показывает, что удобно просмотреть тензор $c_{i j k l}$ как $N \times N$ матрица $2 \times 2$ подматрицы. Самая удобная форма для $c$ входной параметр для этого симметричного, $N = 3$ система имеет 21 строку в $c$ и описан подробно в документации Тулбокса УЧП. Это повторяется здесь для удобства.

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

Пьезоэлектрическая биморфная модель привода

Теперь, когда мы определили уравнения для 2D пьезоэлектрического материала, мы готовы применить их к определенной модели. Модель является консольным лучом 2D слоя, который был экстенсивно изучен (например, судьи 1 и 2). Это задано как "биморфное", потому что несмотря на то, что оба слоя сделаны из того же Фторида Polyvinylidene (PVDF) материалом в верхнем слое точки направления поляризации вниз (минус направление Y) и в нижнем слое, это подчеркивает. Схематический из консольного луча показан на рисунке ниже.

Этот рисунок не должен масштабироваться; фактическое отношение толщины/длины равняется 100, таким образом, луч является очень тонким. Когда напряжение применяется между более низкими и верхними поверхностями луча, оно отклоняет в направлении Y; один слой сокращается, и другой слой удлиняет. Устройства этого типа могут быть спроектированы, чтобы обеспечить необходимое движение или силу для различных приложений.

Создайте Модель УЧП с тремя зависимыми переменными

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

N = 3;
model = createpde(N);

Создание геометрии

Простая геометрия 2D слоя луча может быть создана путем определения суммы двух прямоугольников.

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

Эти две линии ниже содержат столбцы матрицы описания геометрии (GDM) для двух прямоугольных слоев. GDM является первым входным параметром к decsg и описывает основные геометрические сущности в модели.

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

Создайте сущность геометрии и добавьте к модели PDE.

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

Свойства материала и содействующая спецификация

Материал в обоих слоях луча является Фторидом 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];

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

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

Преобразуйте диэлектрическую матрицу от постоянного напряжения до постоянной деформации

D_const_strain = D_const_stress - pzeD'*pzeE;

Как уже отмечалось выше, удобно просмотреть этот 21 коэффициент, требуемый assempde как 3 x 3 массива 2 x 2 подматрицы. cij матрицы, заданные ниже, являются 2 x 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);

Определение граничного условия

В данном примере главному ребру геометрии (ребро 1) предписали напряжение как 100 вольт. Нижнему ребру геометрии (ребро 2) предписали напряжение как 0 вольт (т.е. основанный). Левое ребро геометрии (ребра 6 и 7) имеет u и v смещения равный нуль (т.е. зафиксированный). Напряжение и заряд являются нулем на правильном ребре геометрии (т.е. q = 0).

V = 100;

Установите напряжение (компонент решения 3) на верхнем краю к V.

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

Установите напряжение (компонент решения 3) на базовом краю обнулять.

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

Установите смещения X и Y (компоненты решения 1 и 2) на левом конце (ребра геометрии 6 и 7) обнулять.

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

Поймайте в сети генерацию

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

hmax = 5e-04;
msh = generateMesh(model,'Hmax',hmax,...
                         'GeometricOrder','quadratic',...
                         'MesherVersion','R2013a');

Решение для конечного элемента

result = solvepde(model);

Извлеките NodalSolution свойство от результата, это имеет x-отклонение в столбце 1, y-отклонение в столбце 2 и электрический потенциал в столбце 3. Найдите минимальное y-отклонение и постройте компоненты решения.

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

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

Аналитическое решение

Простое, аппроксимированное, аналитическое решение было получено для этой проблемы в ссылке 1.

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

Analytical tip deflection is:  -3.3000e-05

Сводные данные

Цветные контурные графики x-отклонения и y-отклонения показывают стандартное поведение классического консольного решения для луча. Линейное распределение напряжения через толщину луча как ожидалось. Существует хорошее соглашение между решением для конечного элемента Тулбокса УЧП и аналитическим решением из ссылки 1.

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

Ссылки

Хван, W. S.; Припаркуйтесь, H. C; Моделирование Конечного элемента Пьезоэлектрических Датчиков и Приводов. Журнал AIAA, 31 (5), стр 930-937, 1993.
Пифорд, V; моделирование конечного элемента пьезоэлектрических активных структур. Диссертация, Брюссельский свободный университет, 2001.

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