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

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

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

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

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

-σ=f

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

D=ρ

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

-{σD}={f-ρ}

Для 2D анализа, σ имеет компоненты σ11,σ22, и σ12=σ21, и D имеет компоненты D1 и D2.

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

{σ11σ22σ12D1D2}=[C11C12e11e31C12C22e13e33G12e14e34e11e13e14-E1e31e33e34-E2]{ϵ11ϵ22γ12-E1-E2}

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

{ϵ11ϵ22γ12}={uxvyuy+vx}

Запишите электрическое поле в терминах электрического потенциала ϕ:

{E1E2}=-{ϕxϕy}

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

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

-(cu)+au=f

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

-xk(cijklujxl)+aijuj=fi

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

u={uvϕ}

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

u={uxuyvxvyϕxϕy}

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

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

[c(1)c(2)c(4)c(6)c(11)c(13)c(3)c(5)c(7)c(12)c(14)c(8)c(9)c(15)c(17)c(10)c(16)c(18)c(19)c(20)c(21)]

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

[c1111c1112c1211c1212c1311c1312c1122c1221c1222c1321c1322c2211c2212c2311c2312c2222c2321c2322c3311c3312c3322]{uxuyvxvyϕxϕy}

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

[C11C12e11e31G12G12e14e34G12e14e34C22e13e33-E1-E2]{uxuyvxvyϕxϕy}

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

Создайте модель 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

Ссылки

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

  2. Пифорд, V. "Моделирования конечного элемента Пьезоэлектрических Активных Структур". Доктор философии diss., Брюссельский свободный университет, 2001.

Для просмотра документации необходимо авторизоваться на сайте