Задача рассеяния

Открыть 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$

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

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

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

numberOfPDE = 1;
model = createpde(numberOfPDE);

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

g: Функция спецификации геометрии. Для получения дополнительной информации смотрите раздел документации Параметризованная функция для создания 2-D геометрии и код для 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])

Figure contains an axes. The axes 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. The axes 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. The axes 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) команда.

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

Задача рассеяния

Partial Differential Equation Toolbox производными

Поддержка