Проблема рассеяния

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

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

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

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

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

Теперь можно переписать волновое уравнение как уравнение Гельмгольца для пространственной составляющей отраженной волны с волновым числом $\mathit{k}=$λ/α:

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

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

В этом уравнении учитываются только волны, движущиеся в положительном направлении. В данном случае, start- радиальное расстояние от объекта. При гармоническом решении по времени это уравнение превращается в обобщённое граничное условие Неймана

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

numberOfPDE = 1;
model = createpde(numberOfPDE);

Укажите переменные, определяющие проблему:

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

Создайте сетку.

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) команда.