Создайте символьный УЧП, который представляет Лапласиан переменной u(x,y).

syms u(x,y) f
pdeeq = laplacian(u,[x y]) == f

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)=f$$

Извлеките коэффициенты УЧП как символьные выражения и отобразите их значения.

symCoeffs = pdeCoefficients(pdeeq,u,'Symbolic',true);
structfun(@disp,symCoeffs)

$$0$$

$$0$$

$$\left(\begin{array}{cc}-1& 0\\ 0& -1\end{array}\right)$$

$$f$$

$$0$$

pdeCoefficients преобразует символьный УЧП в скалярное уравнение PDE формы

$$m\frac{{\partial }^{2}u}{\partial t^2}+d\frac{\partial u}{\partial t}-\nabla \cdot (c\nabla u)+au=f$$,

и извлеките коэффициенты aCMD, и f в структуру symCoeffs.

Выберите значение для f. Начиная с symCoeffs символьные объекты, используют pdeCoefficientsToDouble преобразовывать коэффициенты в double тип данных. Коэффициенты с double тип данных является допустимыми входными параметрами для specifyCoefficients функция в Тулбоксе УЧП.

symCoeffs = subs(symCoeffs,f,-3);
coeffs = pdeCoefficientsToDouble(symCoeffs)

coeffs = struct with fields:
    a: 0
    c: [4x1 double]
    m: 0
    d: 0
    f: -3

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

Система УЧП представляет отклонение зафиксированной структурной пластины при универсальной загрузке давления. Система УЧП с зависимыми переменными $$u_1$$ и $$u_2$$ дают

$$-{\nabla }^2{u}_1+u_2=0$$,

$$-D{\nabla }^2{u}_2=p$$,

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

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

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

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

model = createpde(2);

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

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

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

E = 1.0e6;
h_thick = 0.1;
nu = 0.3;
D = E*h_thick^3/(12*(1 - nu^2));
syms pres

Объявите систему УЧП как систему символьные уравнения. Извлеките коэффициенты УЧП и возвратите их в символьной форме.

syms u1(x,y) u2(x,y)
pdeeq = [-laplacian(u1) + u2; -D*laplacian(u2) - pres];
symCoeffs = pdeCoefficients(pdeeq,[u1 u2],'Symbolic',true)

symCoeffs = struct with fields:
    m: [1x1 sym]
    a: [2x2 sym]
    c: [4x4 sym]
    f: [2x1 sym]
    d: [1x1 sym]

Отобразите коэффициенты mACF, и d.

structfun(@disp,symCoeffs)

$$0$$

$$\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$$

$$\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& \frac{25000}{273}& 0\\ 0& 0& 0& \frac{25000}{273}\end{array}\right)$$

$$\left(\begin{array}{c}0\\ \mathrm{pres}\end{array}\right)$$

$$0$$

Замените значением pres использование subs функция. Начиная с выходных параметров subs символьные объекты, используют pdeCoefficientsToDouble функционируйте, чтобы преобразовать коэффициенты в double тип данных, который делает их допустимыми входными параметрами для Тулбокса УЧП.

symCoeffs = subs(symCoeffs,pres,2);
coeffs = pdeCoefficientsToDouble(symCoeffs)

coeffs = struct with fields:
    a: [4x1 double]
    c: [16x1 double]
    m: 0
    d: 0
    f: [2x1 double]

Задайте коэффициенты УЧП для модели PDE.

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

Задайте пружинную жесткость. Задайте граничные условия путем определения распределенных пружин на всех четырех ребрах.

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

Задайте размер mesh геометрии и сгенерируйте mesh для модели PDE.

hmax = len/20;
generateMesh(model,'Hmax',hmax);

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

res = solvepde(model);

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

u = res.NodalSolution;

Постройте поперечное отклонение пластины.

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

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

wMax = min(u(1:numNodes,1))

wMax = -0.2763

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

pres = 2;
wMax = -.0138*pres*len^4/(E*h_thick^3)

wMax = -0.2760