scattering problem -凯发k8网页登录

solve a simple scattering problem, where you compute the waves reflected by an object illuminated by incident waves. this example shows how to solve a scattering problem using the electromagnetic workflow. for the general pde workflow, see .

for this problem, assume that the domain is an infinite horizontal membrane that is subjected to small vertical displacements. the membrane is fixed at the object boundary. the medium is homogeneous, and the phase velocity (propagation speed) of a wave, α, is constant. for this problem, assume α=1. in this case, the frequency wave number equals the frequency, k=ω/α=ω.

to solve the scattering problem, first create an electromagnetic model for harmonic analysis.

emagmodel = createpde("electromagnetic","harmonic");

specify the frequency as 4π.

omega = 4*pi;

represent the square surface with a diamond-shaped hole. define a diamond in a square, place them in one matrix, and create a set formula that subtracts the diamond from the square.

square = [3; 4; -5; -5; 5; 5; -5; 5; 5; -5];
diamond = [2; 4; 2.1; 2.4; 2.7; 2.4; 1.5; 1.8; 1.5; 1.2];
gd = [square,diamond];
ns = char('square','diamond')';
sf = 'square - diamond';

create the geometry.

g = decsg(gd,sf,ns);
geometryfromedges(emagmodel,g);

include the geometry in the model and plot it with the edge labels.

figure; 
pdegplot(emagmodel,"edgelabels","on"); 
xlim([-6,6])
ylim([-6,6])

specify the vacuum permittivity and permeability values as 1.

emagmodel.vacuumpermittivity = 1;
emagmodel.vacuumpermeability = 1;

specify the relative permittivity, relative permeability, and conductivity of the material.

electromagneticproperties(emagmodel,"relativepermittivity",1, ...
                                    "relativepermeability",1, ...
                                    "conductivity",0);

apply the absorbing boundary condition on the edges of the square. specify the thickness and attenuation rate for the absorbing region by using the thickness, exponent, and scaling arguments.

electromagneticbc(emagmodel,"edge",[1 2 7 8], ...
                            "farfield","absorbing", ...
                            "thickness",2, ...
                            "exponent",4, ...
                            "scaling",1);

apply the boundary condition on the diamond edges.

innerbcfunc = @(location,~) [-exp(-1i*omega*location.x); ...
                            zeros(1,length(location.x))];
binner = electromagneticbc(emagmodel,"edge",[3 4 5 6], ...
                                     "electricfield",innerbcfunc);

generate a mesh.

generatemesh(emagmodel,"hmax",0.1);

solve the harmonic analysis model for the frequency ω=4π.

result = solve(emagmodel,"frequency",omega);

plot the real part of the x-component of the resulting electric field.

u = result.electricfield;
figure
pdeplot(emagmodel,"xydata",real(u.ex),"mesh","off");
colormap(jet)

interpolate the resulting electric field to a grid covering the portion of the geometry, for x and y from -1 to 4.

v = linspace(-1,4,101);
[x,y] = meshgrid(v);
eintrp = interpolateharmonicfield(result,x,y);

reshape eintrp.ex and plot the x-component of the resulting electric field.

eintrpx = reshape(eintrp.electricfield.ex,size(x));
figure
surf(x,y,real(eintrpx),"linestyle","none");
view(0,90)
colormap(jet)

using the solutions for a vector of frequencies, create an animation showing the corresponding solution to the time-dependent wave equation.

result = solve(emagmodel,"frequency",omega/10:omega);
figure
for m = 1:length(omega/10:omega)
    u = result.electricfield;
    pdeplot(emagmodel,"xydata",real(u.ex(:,m)),"mesh","off");
    colormap(jet)
    m(m) = getframe;
end

to play the animation, use the movie function. for example, to play the animation five times in a loop with 3 frames per second, use the movie(m,5,3) command.