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

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

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

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

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

Для просмотра документации необходимо авторизоваться на сайте