В этом примере показано, как вычислить отклонение структурной пластины при загрузке давления.
Дифференциальное уравнение с частными производными для тонкой изотропной пластины с загрузкой давления
где изгибающаяся жесткость пластины, данной
и модуль эластичности, отношение Пуассона, толщина пластины, поперечное отклонение пластины, и загрузка давления.
Граничные условия для зафиксированных контуров и , где производная в направлении, нормальном к контуру.
Partial Differential Equation Toolbox™ не может непосредственно решить это уравнение пластины четвертого порядка. Преобразуйте уравнение четвертого порядка в эти два дифференциальных уравнения с частными производными второго порядка, где новая зависимая переменная.
Вы не можете непосредственно задать граничные условия для обоих и в этой системе второго порядка. Вместо этого задайте это 0, и задать так, чтобы также равняется 0 на контуре. Чтобы задать эти условия, используйте жесткие "пружины", распределенные вдоль контура. Пружины прикладывают поперечную силу сдвига к ребру пластины. Задайте силу сдвига вдоль контура из-за этих пружин как , где нормальное к контуру, и жесткость пружин. Это выражение является обобщенным Неймановым граничным условием, поддержанным тулбоксом. Значение должно быть достаточно большим так, чтобы приблизительно 0 во всех точках на контуре. Это также должно быть малым достаточно, чтобы избежать числовых ошибок из-за плохо обусловленной матрицы жесткости.
Тулбокс использует зависимые переменные и вместо и . Перепишите две переменные использования дифференциальных уравнений с частными производными второго порядка и :
Создайте модель PDE для системы двух уравнений.
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) следуют непосредственно из определения системы 2D уравнения.
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
Сравните результат с deflecion в центре пластины, вычисленном аналитически.
wMax = -.0138*pres*len^4/(E*thick^3)
wMax = -0.2760