Зажатый квадратный изотропный диск с равномерным давлением

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

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

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

и $$E$$ - модуль упругости, $$\nu$$ - коэффициент Пуассона, $$h$$ - толщина пластины, $$w$$ - поперечное отклонение пластины, и $$p$$ - нагрузка под давлением.

Граничные условия для зажатых контуров: $$w=0$$ и $$w^{\prime }=0$$, где $$w^{\prime }$$ является производной от $$w$$ в направлении, перпендикулярном контуру.

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

Вы не можете непосредственно задать граничные условия для обоих $$w$$ и $$w^{\prime }$$ в этой системе второго порядка. Вместо этого укажите, что $$w^{\prime }$$ равен 0 и задает $$v^{\prime }$$ так что $$w$$ также равен 0 на контуре. Чтобы указать эти условия, используйте жесткие «пружины», распределенные вдоль контура. Пружины прикладывают поперечную силу сдвига к ребру диска. Определите силу сдвига вдоль контура от этих пружин как $$n\cdot D\nabla v=-kw$$, где $$n$$ является нормалью к контуру, и $$k$$ - жесткость пружин. Это выражение является обобщенным граничным условием Неймана, поддерживаемым тулбоксом. Значение $$k$$ должен быть достаточно большим, чтобы $$w$$ равен приблизительно 0 во всех точках на контуре. Он также должен быть достаточно маленьким, чтобы избежать числовых ошибок из-за плохо обусловленной матрицы жесткости.

Тулбокс использует зависимые переменные ${\mathit{u}}_1$ и ${\mathit{u}}_2$ вместо $$w$$ и $$v$$. Перепишите два дифференциальных уравнений с частными производными производными второго порядка с помощью переменных ${\mathit{u}}_1$ и ${\mathit{u}}_2$:

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

model = createpde(2);

Создайте квадратную геометрию и включите ее в модель.

len = 10;
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'

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

Коэффициент c в этом примере является тензором, который может быть представлен как матрица 2 на 2 блоков 2 на 2:

Эта матрица дополнительно уплощена в вектор-столбец из шести элементов. Значения в полной матрице 2 на 2 (определяющие коэффициент a) и векторе 2 на 1 (определяющем коэффициент f) следуют непосредственно из определения системы двух уравнений.

E = 1.0e6; % Modulus of elasticity
nu = 0.3; % Poisson's ratio
thick = 0.1; % Plate thickness
pres = 2; % External pressure

D = E*thick^3/(12*(1 - nu^2));

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;

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

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

Сгенерируйте mesh.

generateMesh(model);

Решить модель.

res = solvepde(model);

Доступ к решению в узловых местоположениях.

u = res.NodalSolution;

Постройте график поперечного отклонения.

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

Найдите поперечное отклонение в центре диска.

numNodes = size(model.Mesh.Nodes,2);
wMax = min(u(1:numNodes,1))

wMax = -0.2763

Сравните результат с отклонением в центре диска, вычисленным аналитически.

wMax = -.0138*pres*len^4/(E*thick^3)

wMax = -0.2760