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

Настройка и решение УЧП

Задайте функцию для нелинейного f коэффициент в синтаксисе, как задано в f Коэффициент для specifyCofficients.

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);

Создать mesh со сторонами не более 0,3.

generateMesh(model,'Hmax',0.3);

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

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])

Настройка и решение УЧП

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

Partial Differential Equation Toolbox™ решает уравнения в виде

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

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

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

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

Задайте коэффициенты УЧП

Включите коэффициенты УЧП в model.

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

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

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

Создайте mesh и решите УЧП.

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

Постройте график решения как срезы контура

Потому что граничные условия $$u=0$$ на всех лицах, решение $$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 сетки. Кроме того, может потребоваться много времени для оценки решения на полной сетке. Чтобы сэкономить память и время, оценивайте только в точках, которые вы строите. Можно также использовать этот метод для интерполяции с наклонными сетками или с другими поверхностями.

Например, интерполируйте решение сетки на наклонной плоскости $$0\le x\le 135$$, $$0\le y\le 35$$, и $$z=x/10+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=0$$. Сохраните граничные условия по умолчанию на гранях 3, 4 и 6.

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

Создание коэффициентов УЧП

УЧП:

В Partial Differential Equation 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);

Создание Mesh и решение УЧП

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

generateMesh(model);

Решить УЧП.

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)=(100,0,0)$$.

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

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

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

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

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.