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

Настройка и решение PDE

Определение функции для нелинейного f коэффициент в синтаксисе, как дали в f Коэффициенте для specifyCoefficients.

function bcMatrix = myfffun(region,state)

bcMatrix = 1+10*region.z.^2+region.y;

Импортируйте геометрию и исследуйте метки граней.

model = createpde;
g = importGeometry(model,'Plate10x10x1.stl');
pdegplot(g,'FaceLabels','on','FaceAlpha',0.5)

Грани пронумерованы с 1 по 6.

Создайте коэффициенты и граничные условия.

c = 1;
a = 0;
d = 1;
f = @myfffun;
specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f);

applyBoundaryCondition(model,'dirichlet','face',1:6,'u',0);

Задайте нулевое начальное условие.

setInitialConditions(model,0);

Создайте сетку со сторонами не более 0,3.

generateMesh(model,'Hmax',0.3);

Установите время от 0 до 0,2 и решите PDE.

tlist = 0:0.02:0.2;
results = solvepde(model,tlist);

Печать фрагментов в решении

Создание сетки из (x,y,z) точки, где x = 5, y находится в диапазоне от 0 до 10, и z находится в диапазоне от 0 до 1. Интерполяция решения к этим точкам сетки и все время.

yy = 0:0.5:10;
zz = 0:0.25:1;
[YY,ZZ] = meshgrid(yy,zz);
XX = 5*ones(size(YY));
uintrp = interpolateSolution(results,XX,YY,ZZ,1:length(tlist));

Матрица решения uintrp имеет 11 столбцов, по одному для каждого раза в tlist. Возьмем интерполированное решение для второго столбца, что соответствует времени 0.02.

usol = uintrp(:,2);

Элементы usol происходит от интерполяции решения к XX, YY, и ZZ матрицы, каждая из которых 5 на 21, соответствующие z-by-y переменные. Измениться usol до того же размера 5 на 21 и составить график поверхности решения. Также создайте графики поверхности, соответствующие временам 0,06, 0,10 и 0,20.

figure
usol = reshape(usol,size(XX));
subplot(2,2,1)
surf(usol)
title('t = 0.02')
zlim([0,1.5])
xlim([1,21])
ylim([1,5])

usol = uintrp(:,4);
usol = reshape(usol,size(XX));
subplot(2,2,2)
surf(usol)
title('t = 0.06')
zlim([0,1.5])
xlim([1,21])
ylim([1,5])

usol = uintrp(:,6);
usol = reshape(usol,size(XX));
subplot(2,2,3)
surf(usol)
title('t = 0.10')
zlim([0,1.5])
xlim([1,21])
ylim([1,5])

usol = uintrp(:,11);
usol = reshape(usol,size(XX));
subplot(2,2,4)
surf(usol)
title('t = 0.20')
zlim([0,1.5])
xlim([1,21])
ylim([1,5])

Настройка и решение PDE

Задача состоит в решении уравнения Пуассона с нулевыми граничными условиями Дирихле для сложной геометрии. Уравнение Пуассона

Дифференциальное уравнение в частных производных Toolbox™ решает уравнения в форме

Таким образом, вы можете представить проблему, установив $$c\mathrm{=}$$ 1 $$и\mathrm{a}$$= 0. Произвольно $$\mathrm{задается}\mathrm{}$$f = 10.

c = 1;
a = 0;
f = 10;

Первым шагом в решении любой 3-D проблемы PDE является создание модели PDE. Это контейнер, содержащий количество уравнений, геометрию, сетку и граничные условия для PDE. Создайте модель, а затем импортируйте 'ForearmLink.stl' файл и просмотр геометрии.

N = 1;
model = createpde(N);
importGeometry(model,'ForearmLink.stl');
pdegplot(model,'FaceAlpha',0.5)
view(-42,24)

Определение коэффициентов PDE

Включить коэффициенты PDE в model.

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

Создайте нулевые граничные условия Дирихле на всех гранях.

applyBoundaryCondition(model,'dirichlet','Face',1:model.Geometry.NumFaces,'u',0);

Создайте сетку и решите PDE.

generateMesh(model);
result = solvepde(model);

Печать решения в виде контурных фрагментов

Поскольку граничные условия равны $$u\mathrm{=}$$ 0 на всех гранях, $$\mathrm{решение}$$ u ненулевое только внутри. Для проверки внутреннего пространства следует использовать прямоугольную сетку, охватывающую геометрию с интервалом в одну единицу в каждом направлении координат.

[X,Y,Z] = meshgrid(0:135,0:35,0:61);

Для печати и анализа создайте PDEResults объект из решения. Интерполировать результат в каждой точке сетки.

V = interpolateSolution(result,X,Y,Z);
V = reshape(V,size(X));

Печать фрагментов контура для различных значений $$z$$.

figure
colormap jet
contourslice(X,Y,Z,V,[],[],0:5:60)
xlabel('x')
ylabel('y')
zlabel('z')
colorbar
view(-11,14)
axis equal

Печать фрагментов контура для различных значений $$y$$.

figure
colormap jet
contourslice(X,Y,Z,V,[],1:6:31,[])
xlabel('x')
ylabel('y')
zlabel('z')
colorbar
view(-62,34)
axis equal

Сохранение памяти путем оценки по мере необходимости

Для больших проблем при создании тонкой сетки 3-D может не хватать памяти. Кроме того, оценка решения по полной сетке может занять много времени. Чтобы сэкономить память и время, вычисляйте только в тех точках, которые вы выводите на график. Этот метод можно также использовать для интерполяции к наклонным сеткам или к другим поверхностям.

Например, интерполировать решение в сетку на наклонной плоскости $$\mathrm{}\mathrm{}\mathrm{}\mathrm{0\le x\le 135}$$, $$\mathrm{}\mathrm{}\mathrm{0\le y\le 35}$$и $$z=\mathrm{}\mathrm{}\mathrm{x/10}\mathrm{+}$$y/2. Печать данных горизонталей и цветных поверхностей. Используйте тонкую сетку с интервалом 0,2.

[X,Y] = meshgrid(0:0.2:135,0:0.2:35);
Z = X/10 + Y/2;
V = interpolateSolution(result,X,Y,Z);
V = reshape(V,size(X));
figure
subplot(2,1,1)
contour(X,Y,V);
axis equal
title('Contour Plot on Tilted Plane')
xlabel('x')
ylabel('y')
colorbar
subplot(2,1,2)
surf(X,Y,V,'LineStyle','none');
axis equal
view(0,90)
title('Colored Plot on Tilted Plane')
xlabel('x')
ylabel('y')
colorbar

Импорт и просмотр геометрии

model = createpde;
importGeometry(model,'Block.stl');
pdegplot(model,'FaceLabels','on','FaceAlpha',0.5)
view(-42,24)

Задать граничные условия

Задайте грани 1, 2 и 5 в качестве мест, где частица может сбежать. На этих гранях решение $$u\mathrm{=}$$0. Сохранить граничные условия по умолчанию для граней 3, 4 и 6.

applyBoundaryCondition(model,'dirichlet','Face',[1,2,5],'u',0);

Создание коэффициентов PDE

PDE:

В синтаксисе Toolbox™ дифференциального уравнения в частных производных

Это уравнение переводится в коэффициенты c = 1, a = 0, и f = 2. Введите коэффициенты.

c = 1;
a = 0;
f = 2;
specifyCoefficients(model,'m',0,'d',0,'c',c','a',a,'f',f);

Создание сетки и решение PDE

Инициализируйте сетку.

generateMesh(model);

Решите проблему PDE.

results = solvepde(model);

Анализ решения на графике контурного фрагмента

Создайте сетку и интерполируйте решение в сетку.

[X,Y,Z] = meshgrid(0:135,0:35,0:61);
V = interpolateSolution(results,X,Y,Z);
V = reshape(V,size(X));

Создание графика горизонталей для пяти фиксированных значений y-координата.

figure
colormap jet
contourslice(X,Y,Z,V,[],0:4:16,[])
xlabel('x')
ylabel('y')
zlabel('z')
xlim([0,100])
ylim([0,20])
zlim([0,50])
axis equal
view(-50,22)
colorbar

Частица имеет наибольшее среднее время выхода вблизи точки $$(x,y,z\left)\mathrm{}\mathrm{=}\mathrm{}\right(\mathrm{}\mathrm{}$$100,0,0).

Использование градиентов для скрученных и упорядоченных графиков

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

[X,Y,Z] = meshgrid(1:9:99,1:3:20,1:6:50);
[gradx,grady,gradz] = evaluateGradient(results,X,Y,Z);

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

gradx = reshape(gradx,size(X));
grady = reshape(grady,size(Y));
gradz = reshape(gradz,size(Z));

figure
quiver3(X,Y,Z,gradx,grady,gradz)
axis equal
xlabel 'x'
ylabel 'y'
zlabel 'z'
title('Quiver Plot of Estimated Gradient of Solution')

Постройте график аппроксимации градиента. Начните обтекание из более разреженного набора начальных точек.

hold on
[sx,sy,sz] = meshgrid([1,46],1:6:20,1:12:50);
streamline(X,Y,Z,gradx,grady,gradz,sx,sy,sz)
title('Quiver Plot with Streamlines')
hold off

Оптимизация показывает, что небольшие значения y и z дать большее среднее время выхода. Они также показывают, что x-координата оказывает значительное влияние на u когда x мал, но когда x больше 40, большие значения мало влияют на u. Аналогично, когда z меньше 20, его значения мало влияют на u.