Зафиксированная, квадратная изотропная пластина с универсальной загрузкой давления

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

УЧП и граничные условия для тонкой пластины

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

$\nabla^{2} (D \nabla^{2} w) = - p$

где $D$ изгибающаяся жесткость пластины, данной

$D = \frac{E h^{3}}{12 (1 - ν^{2})}$

и $E$ модуль эластичности, $ν$ отношение Пуассона, и $h$ толщина пластины. Поперечное отклонение пластины $w$ и $p$ загрузка давления.

Граничные условия для зафиксированных контуров $w = 0$ и $w^{'} = 0$ где $w^{'}$ производная $w$ в направлении, нормальном к контуру.

Partial Differential Equation Toolbox™ не может непосредственно решить четвертое уравнение пластины порядка, показанное выше, но это может быть преобразовано в следующие два дифференциальных уравнения с частными производными второго порядка.

$\nabla^{2} w = v$

$D \nabla^{2} v = - p$

где $v$ новая зависимая переменная. Однако не очевидно, как задать граничные условия для этой системы второго порядка. Мы не можем непосредственно задать граничные условия для обоих $w$ и $w^{'}$ . Вместо этого мы непосредственно предписываем $w^{'}$ быть нулем и использовать следующий метод, чтобы задать $v^{'}$ таким способом обеспечить это $w$ также равняется нулю на контуре. Жесткие "пружины", которые прикладывают поперечную силу сдвига к ребру пластины, распределяются вдоль контура. Сила сдвига вдоль контура из-за этих пружин может быть записана $n \cdot D \nabla v = - k w$ где $n$ нормальное к контуру и $k$ жесткость пружин. Значение $k$ должно быть достаточно большим это $w$ приблизительно нуль во всех точках на контуре, но не столь большой, что числовые ошибки заканчиваются, потому что матрица жесткости плохо обусловлена. Это выражение является обобщенным Неймановым граничным условием, поддержанным Partial Differential Equation Toolbox™

В определении Partial Differential Equation Toolbox™ для эллиптической системы, $w$ и $v$ зависимые переменные являются you(1), и u (2). Два дифференциальных уравнения с частными производными второго порядка могут быть переписаны как

$- \nabla^{2} u_{1} + u_{2} = 0$

$- D \nabla^{2} u_{2} = p$

который является формой, поддержанной тулбоксом. Вход, соответствующий этой формулировке, показывают в разделах ниже.

Создайте модель УЧП

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

numberOfPDE = 2;
model = createpde(numberOfPDE);

Проблемные параметры

E = 1.0e6; % modulus of elasticity
nu = .3; % Poisson's ratio
thick = .1; % plate thickness
len = 10.0; % side length for the square plate
hmax = len/20; % mesh size parameter
D = E*thick^3/(12*(1 - nu^2));
pres = 2; % external pressure

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

Для одного квадрата геометрия и mesh легко заданы как показано ниже.

gdm = [3 4 0 len len 0 0 0 len len]';
g = decsg(gdm,'S1',('S1')');

Создайте сущность геометрии.

geometryFromEdges(model,g);

Постройте геометрию и отобразите метки ребра для использования в определении граничного условия.

figure; 
pdegplot(model,'EdgeLabels','on');
ylim([-1,11])
axis equal
title 'Geometry With Edge Labels Displayed';

Содействующее определение

Документация относительно коэффициентов УЧП показывает требуемые форматы для a и c матриц. Самая удобная форма для c в этом примере $n_{c} = 3 N$ из таблицы, где $N$ количество дифференциальных уравнений. В этом примере $N = 2$ . $c$ тензор, в форме $N \times N$ матрица $2 \times 2$ подматрицы показывают ниже.

$[\begin{array}{cccc} c (1) & c (2) & \cdot & \cdot \\ \cdot & c (3) & \cdot & \cdot \\ \cdot & \cdot & c (4) & c (5) \\ \cdot & \cdot & \cdot & c (6) \end{array}]$

С шестью строками матрицей одного столбца c задан ниже. Записи в полном $2 \times 2$ матрица и $2 \times 1$ f вектор следуют непосредственно из определения системы 2D уравнения, показанной выше.

c = [1 0 1 D 0 D]';
a = [0 0 1 0]';
f = [0 pres]';
specifyCoefficients(model,'m',0,'d',0,'c',c,'a',a,'f',f);

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

k = 1e7; % spring stiffness

Задайте распределенные пружины на всех четырех ребрах.

bOuter = applyBoundaryCondition(model,'neumann','Edge',(1:4),...
                                     'g',[0 0],'q',[0 0; k 0]);

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

generateMesh(model, 'Hmax', hmax);

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

Решение вычисляется с помощью solvepde функционируйте и поперечное отклонение построено с помощью pdeplot функция. Для сравнения поперечное отклонение в центре пластины также вычисляется с помощью аналитического решения этой проблемы.

res = solvepde(model);
u = res.NodalSolution;
numNodes = size(model.Mesh.Nodes,2);
figure
pdeplot(model,'XYData',u(1:numNodes),'Contour','on');
title 'Transverse Deflection'

numNodes = size(model.Mesh.Nodes,2);
fprintf('Transverse deflection at plate center(PDE Toolbox) = %12.4e\n',...
                                                  min(u(1:numNodes,1)));

Transverse deflection at plate center(PDE Toolbox) =  -2.7632e-01

Вычислите аналитическое решение.

wMax = -.0138*pres*len^4/(E*thick^3);
fprintf('Transverse deflection at plate center(analytical) = %12.4e\n', wMax);

Transverse deflection at plate center(analytical) =  -2.7600e-01

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