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

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

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

В этом примере модель является консольным лучом 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

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

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

Ссылки

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

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