Решение PDE с неконстантными граничными условиями

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

Геометрия

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

% Rectangle is code 3, 4 sides, followed by x-coordinates and then y-coordinates
R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]';
% Circle is code 1, center (.5,0), radius .2
C1 = [1,.5,0,.2]';
% Pad C1 with zeros to enable concatenation with R1
C1 = [C1;zeros(length(R1)-length(C1),1)];
geom = [R1,C1];

% Names for the two geometric objects
ns = (char('R1','C1'))';

% Set formula
sf = 'R1-C1';

% Create geometry
g = decsg(geom,sf,ns);

% Create geometry model
model = createpde;

% Include the geometry in the model and view the geometry
geometryFromEdges(model,g);
pdegplot(model,'EdgeLabels','on')
xlim([-1.1 1.1])
axis equal

Figure contains an axes. The axes contains 9 objects of type line, text.

Скалярная задача

Ребро 3 имеет условия Дирихле со значением 32.
Ребро 1 имеет условия Дирихле со значением 72.
Ребра 2 и 4 имеют условия Дирихле, которые линейно интерполируются между ребрами 1 и 3.
Круговые ребра (с 5 по 8) имеют условия Неймана с q = 0, g = -1.

applyBoundaryCondition(model,'dirichlet','Edge',3,'u',32);
applyBoundaryCondition(model,'dirichlet','Edge',1,'u',72);
applyBoundaryCondition(model,'neumann','Edge',5:8,'g',-1); % q = 0 by default

Ребрам 2 и 4 нужны функции, которые выполняют линейную интерполяцию. Каждое ребро может использовать ту же функцию, которая возвращает значение $u (x, y) = 52 + 20 x$ .

Реализовать эту простую интерполяцию можно в анонимной функции.

myufunction = @(location,state)52 + 20*location.x;

Включите функцию для ребер 2 и 4. Чтобы помочь ускорить решатель, разрешите векторизованную оценку.

applyBoundaryCondition(model,'dirichlet','Edge',[2,4],...
                                         'u',myufunction,...
                                         'Vectorized','on');

Решить эллиптический УЧП с этими граничными условиями, используя параметры c = 1, a = 0, и | f = 10|. Поскольку более короткая прямоугольная сторона имеет длину 0,8, чтобы гарантировать, что mesh не слишком груба, выберите максимальный размер сетки Hmax = 0.1.

specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',10);
generateMesh(model,'Hmax',0.1);
results = solvepde(model);
u = results.NodalSolution;
pdeplot(model,'XYData',u)

Figure contains an axes. The axes contains an object of type patch.

Система PDE

Предположим, что система имеет N = 2.

Ребро 3 имеет условия Дирихле со значениями [32,72].
Ребро 1 имеет условия Дирихле со значениями [72,32].
Ребра 2 и 4 имеют условия Дирихле, которые интерполируются между условиями на ребрах 1 и 3 и включают синусоидальное изменение.
Круговые ребра (ребра с 5 по 8) имеют q = 0 и g = -10.

model = createpde(2);
geometryFromEdges(model,g);

applyBoundaryCondition(model,'dirichlet','Edge',3,'u',[32,72]);
applyBoundaryCondition(model,'dirichlet','Edge',1,'u',[72,32]);
applyBoundaryCondition(model,'neumann','Edge',5:8,'g',[-10,-10]);

Первый компонент ребер 2 и 4 удовлетворяет уравнению $u_{1} (x) = 52 + 20 x + 10 \sin (π x^{3})$ .

Второй компонент удовлетворяет $u_{2} (x) = 52 - 20 x - 10 \sin (π x^{3})$ .

Запись файла функции myufun.m который включает эти уравнения в синтаксис, описанный в неконстантных граничных условиях.

function bcMatrix = myufun(location,state)
bcMatrix = [52 + 20*location.x + 10*sin(pi*(location.x.^3));
    52 - 20*location.x - 10*sin(pi*(location.x.^3))]; % OK to vectorize
end

Включите эту функцию в граничное условие ребра 2 и ребра 4.

applyBoundaryCondition(model,'dirichlet','Edge',[2,4],...
                                         'u',@myufun,...
                                         'Vectorized','on');

Решить эллиптический УЧП с этими граничными условиями, с параметрами c = 1, a = 0, и f = (10,-10). Поскольку более короткая прямоугольная сторона имеет длину 0,8, чтобы гарантировать, что mesh не слишком груба, выберите максимальный размер сетки Hmax = 0.1.

specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',[10;-10]);
generateMesh(model,'Hmax',0.1);
results = solvepde(model);
u = results.NodalSolution;

subplot(1,2,1)
pdeplot(model,'XYData',u(:,1),'ZData',u(:,1),'ColorBar','off')
view(-9,24)
subplot(1,2,2)
pdeplot(model,'XYData',u(:,2),'ZData',u(:,2),'ColorBar','off')
view(-9,24)

Figure contains 2 axes. Axes 1 contains an object of type patch. Axes 2 contains an object of type patch.

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

Решение PDE с неконстантными граничными условиями

Геометрия

Скалярная задача

Система PDE

Partial Differential Equation Toolbox производными

Поддержка