Рассеивание проблемы

Этот пример показывает, как решить простую проблему рассеивания, где вы вычисляете волны, отраженные квадратным объектом, освещенным инцидентными волнами, которые прибывают слева. Этот пример использует программируемый рабочий процесс. Для прикладного решения PDE Modeler смотрите Рассеивающуюся проблему: Приложение PDE Modeler.

Для этой проблемы примите бесконечную горизонтальную мембрану, подвергнутую маленьким вертикальным смещениям U. Мембрана фиксируется на границе объекта. Носитель является гомогенным, и фазовая скорость (скорость распространения) волны, α, является постоянной. Уравнение волны

$\frac{\partial^{2} U}{\partial t^{2}} - α^{2} ▵ U = 0$

Решение U является суммой инцидентной волны V и отраженной волны R:

$U = V + R$

Когда освещение является гармоническим вовремя, можно вычислить поле путем решения одной устойчивой проблемы. Примите, что инцидентная волна является плоской волной, перемещающейся в-x направлении:

$V (x, y, t) = e^{i (- k x - ω t)} = e^{- ikx} \cdot e^{- i ω t}$

Отраженная волна может быть разложена на компоненты времени и пространственный:

$R (x, y, t) = r (x, y) e^{- i ω t}$

Теперь можно переписать уравнение волны как уравнение Гельмгольца для пространственного компонента отраженной волны с номером волны $k = ω / α$ :

$- Δ r - k^{2} r = 0$

Граничное условие Дирихле для контура объекта является U = 0, или с точки зрения инцидента и отраженных волн, R =-V. Для гармонического временем решения и инцидентной волны, перемещающейся в-x направлении, можно записать это граничное условие можно следующим образом:

$r (x, y) = - e^{- ikx}$

Отраженная волна R перемещения, исходящие от объекта. Условие на внешнем вычислительном контуре должно позволить волнам передавать без отражения. Такие условия обычно называются, не отражаясь. Как $| \vec{x} |$ бесконечность подходов, R приблизительно удовлетворяет одностороннему уравнению волны

$\frac{\partial R}{\partial t} + α \vec{ξ} \cdot \nabla R = 0$

Это уравнение рассматривает только волны, перемещающиеся в положительный ξ-direction. Здесь, ξ является радиальным расстоянием от объекта. С гармоническим временем решением это уравнение превращается в обобщенное Нейманово граничное условие

$\vec{ξ} \cdot \nabla r - ikr = 0$

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

numberOfPDE = 1;
model = createpde(numberOfPDE);

Задайте переменные, которые задают проблему:

g: функция спецификации геометрии. Для получения дополнительной информации смотрите раздел Parametrized Function документации для 2D Создания Геометрии и код для scatterg.m.
k, c, a, f: коэффициенты и неоднородный термин.

g = @scatterg;
k = 60;
c = 1;
a = -k^2;
f = 0;

Преобразуйте геометрию и добавьте ее к модели.

geometryFromEdges(model,g);

Постройте геометрию и отобразите метки ребра для использования в определении граничного условия.

figure; 
pdegplot(model,'EdgeLabels','on'); 
axis equal
title 'Geometry with Edge Labels Displayed';
ylim([0,1])

Примените граничные условия.

bOuter = applyBoundaryCondition(model,'neumann','Edge',(5:8),'g',0,'q',-60i);
innerBCFunc = @(loc,state)-exp(-1i*k*loc.x);
bInner = applyBoundaryCondition(model,'dirichlet','Edge',(1:4),'u',innerBCFunc);

Задайте коэффициенты.

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

Сгенерируйте mesh.

generateMesh(model,'Hmax',0.02);
figure
pdemesh(model); 
axis equal

Решите для комплексной амплитуды. Действительная часть векторного u хранит приближение к решению с действительным знаком уравнения Гельмгольца.

result = solvepde(model);
u = result.NodalSolution;

Постройте решение.

figure
pdeplot(model,'XYData',real(u),'Mesh','off');
colormap(jet)

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

figure
m = 10;
maxu = max(abs(u));
for j = 1:m
    uu = real(exp(-j*2*pi/m*sqrt(-1))*u);
    pdeplot(model,'XYData',uu,'ColorBar','off','Mesh','off');
    colormap(jet)
    caxis([-maxu maxu]);
    axis tight
    ax = gca;
    ax.DataAspectRatio = [1 1 1]; 
    axis off
    M(j) = getframe;
end

Чтобы проигрывать фильм, используйте команду movie(M).

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

Рассеивание проблемы

Документация Partial Differential Equation Toolbox

Поддержка