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

Скрипт Open Live Script

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

Для этой проблемы примите бесконечную горизонтальную мембрану, подвергнутую маленьким вертикальным смещениям 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.
kCAF: Коэффициенты и неоднородный термин.

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)
xlabel 'x'
ylabel 'y'
title('Real Value Solution of Helmholtz Equation')

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

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

Поддержка