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

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

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

2Ut2-α2U=0

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

U=V+R

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

V(x,y,t)=ei(-kx-ωt)=e-ikxe-iωt

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

R(x,y,t)=r(x,y)e-iωt

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

-Δr-k2r=0

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

r(x,y)=-e-ikx

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

Rt+αξR=0

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

ξr-ikr=0

Чтобы решить рассеивающуюся задачу с помощью программируемого рабочего процесса, сначала создайте модель 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])

Figure contains an axes object. The axes object with title Geometry with Edge Labels Displayed contains 9 objects of type line, text.

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

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

Figure contains an axes object. The axes object contains 2 objects of type line.

Решите для комплексной амплитуды. Действительная часть векторного 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 contains an axes object. The axes object with title Real Value Solution of Helmholtz Equation contains an object of type patch.

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

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