Finite Difference Boundary Condition

Question

Shaikhah Alkhadhr 2021-10-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1575098-finite-difference-boundary-condition

评论： Shaikhah Alkhadhr 2021-10-31

Hello Everyone,

I have implemented a Finite Difference approach to approximate the wave equation with a pulse as a source function. I need some help from someone who is also familiar with finite difference. As you can see in my code I have reflective and absorbing boundary conditions that I switch between( commented and uncommented). When I run my absorbing boundary conditions I can still see significant reflectance happening at the boundary. Please advise on how to perform full absorption at the boundary, or at least major absorption where reflecting effect is negligible.

Thanks everyone!

tic
c = 1;
Cx = 1;
Cy = 1;
Cx2 = Cx^2;
Cy2 = Cy^2;
T = 3;
Lx = 3;
Ly = 3;
dt = 0.01;
Nt = round(T/dt);
t_input = linspace(0, Nt*dt, Nt);
dx = dt*c/Cx;
dy = dt*c/Cy;
Nx = round(Lx/dx);
Ny = round(Ly/dy);
x_input = linspace(0, Lx, Nx);
y_input= linspace(0, Ly, Ny);
CFL = 0.5
% stability_limit = (1/double(c))*(1/sqrt(1/dx^2 + 1/dy^2));
%  if dt <= 0                % max time step?
%      safety_factor = -dt;    % use negative dt as safety factor
%      dt = safety_factor*stability_limit;
%  elseif dt > stability_limit
%     fprintf('error: dt=%.4f exceeds the stability limit %.6f.\n', dt, stability_limit)
% end
dt2 = dt^2;
u   = zeros(Nx,Ny);   % solution array
u_1 = zeros(Nx,Ny);   % solution at t-dt
u_2 = zeros(Nx,Ny);   % solution at t-2*dt
u_output = zeros(Nx, Ny, Nt);    % full solution to export
Ix = 1:Nx;
Iy = 1:Ny;
It = 1:Nt;
%IC
u_1 = fspecial('gaussian',Nx,5);
% Special formula for first time step
n = 1;
% First step requires a special formula
dt = sqrt(dt2);  % save
Cx2 = CFL*Cx2;  Cy2 = CFL*Cy2; dt2 = CFL*dt2;  % redefine
D1 = 1;  D2 = 0;
for i = 2:Nx-1
    for j = 2:Ny-1
        u_xx = u_1(i-1,j) - 2*u_1(i,j) + u_1(i+1,j);
        u_yy = u_1(i,j-1) - 2*u_1(i,j) + u_1(i,j+1);
        u(i,j) = D1*u_1(i,j) - D2*u_2(i,j) + Cx2*u_xx + Cy2*u_yy + ...
            dt2*f1(x_input(i), y_input(j), t_input(n)) + dt2*f2(x_input(i), y_input(j), t_input(n))+...
            dt2*f3(x_input(i), y_input(j), t_input(n)) + dt2*f4(x_input(i), y_input(j), t_input(n));
        u(i,j) = u(i,j)+ dt*V(x_input(i), y_input(j));
    end
end
u_2 = u_1;
u_1 = u;
u_output(:,:,1) = u;
% Advance in steps
D1 = 2;  D2 = 1;
for n = 2:Nt-1
    % Absorbing boundary condition    
    u(1,:) = u_1(2,:) + ((CFL - 1)/(CFL + 1)*(u(2,:) - u_1(1,:)));
    u(end,:) = u_1(end - 1,:) +  ((CFL - 1)/(CFL + 1)*(u(end-1,:) - u_1(end,:)));
    u(:,1) = u_1(:,2) + ((CFL - 1)/(CFL + 1)*(u(:,2) - u_1(:,1)));
    u(:,end) = u_1(:,end - 1) +  ((CFL - 1)/(CFL + 1)*(u(:,end-1) - u_1(:,end)));
    
    %reflection boundary condition
    %Boundary condition u=0
    %     u_1(1,:) = 0;
    %     u_1(end,:) = 0;
    %     u_1(:,1) = 0;
    %     u_1(:,end) = 0;
    
    for i =2:Nx-1
        for j =2:Ny-1
            u_xx = u_1(i-1,j) - 2*u_1(i,j) + u_1(i+1,j);
            u_yy = u_1(i,j-1) - 2*u_1(i,j) + u_1(i,j+1);
            u(i,j) = D1*u_1(i,j) - D2*u_2(i,j) + Cx2*u_xx + Cy2*u_yy; 
        end
    end
    
    u_2 = u_1;
    u_1 = u;
    u_output(:,:,n) = u; 
end
toc
%Visualize
Zlimits = [-0.1, 0.1];
for frame = 1:10:Nt
    s = surf(x_input,y_input,u_output(:,:,frame),'linewidth',2);
    s.EdgeColor='none';
    grid on;
    axis([min(x_input) max(x_input) min(y_input) max(y_input) Zlimits(1) Zlimits(2)]);
    xlabel('X axis','fontSize',14);
    ylabel('Y axis','fontSize',14); 
    zlabel('Wave Amplitude','fontSize',14);
    %zlim=[-0.3 0.3];
    titlestring = [' TIME = ',num2str(t_input(frame)), 'second'];
    title(titlestring ,'fontsize',14);                            
    h=gca; 
    get(h,'FontSize');
    set(h,'FontSize',14);
    fh = figure(6);
    set(fh, 'color', 'white'); 
    F=getframe(fh);
end
movie(fh,F,1);
s = surf(x_input,y_input,u_output(:,:,299),'linewidth',2);
axis([min(x_input) max(x_input) min(y_input) max(y_input) Zlimits(1) Zlimits(2)]);
s.EdgeColor='none';