计算科学 - 使用高阶方法进行边界处理 - 吾爱随笔录

使用高阶方法进行边界处理

计算科学计算物理学泊松

2021-12-02 04:41:43

我写了一个代码来求解二维泊松方程，其中到处都是齐次的狄利克雷 BC，源项为 -1。我正在使用经典的 Jacobi 迭代方法。网格是 $N_x \, \mathrm{x} \, N_y$ , 并且所有的数组都是从 0 开始的。

该代码使用 5 点模板（2 $^{\mathrm{nd}}$ 阶有限差分法）。使用这个模板，代码循环所有内部点（即从 $1$ 到 $N_x-2$ 在里面 $x$ 方向），而边界点永远不会被触及。

但是，当使用 9 点模板时，代码会出现差异，4 $^\mathrm{th}$ 阶中心差分法。同样，仅迭代内部点。由于离散化在每侧使用两个点，因此我使用一阶近似值推断边界点（鬼点），即 $U(-1) = U(0), U(Nx) = U(Nx-1)$ .

我必须在边界附近使用特殊的离散化方法吗？我可以使用非对称有限差分方案，但我希望我的方法（带有鬼点）能够工作。下面是我的代码示例。对于这种情况，网格是等距的（ $\Delta_x = \Delta_y$ ，或者在代码中 $h(1) = h(2)$ ）。

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~               
! Jacobi scheme                                                                                          
if(method_order == second_order) then                                                                    
    ! 2nd order method: 5 points stencil                                                                 
    do jj=min_y,max_y                                                                                    
        do ii=min_x,max_x                                                                                
            Ujacob(ii,jj) = 0.25_dp                                                     &               
                &   * (Uold(ii-1,jj) + Uold(ii+1,jj) + Uold(ii,jj-1) + Uold(ii,jj+1)    &               
                &   - h(1)**2 * Urhs(ii,jj))                                                             
        end do                                                                                           
    end do                                                                                               
else if(method_order == fourth_order) then                                                               

    ! Fill buffer points                                                                                  
    Ujacob(     -1, 0:grid_ny-1) = -Ujacob(        1, 0:grid_ny-1)                                        
    Ujacob(grid_nx, 0:grid_ny-1) = -Ujacob(grid_nx-2, 0:grid_ny-1)                                        

    Ujacob(0:grid_nx-1,      -1) = -Ujacob(0:grid_nx-1,         1)                                        
    Ujacob(0:grid_nx-1, grid_ny) = -Ujacob(0:grid_nx-1, grid_ny-2)  

    ! 4th order method: 9 points stencil                            
    do jj=min_y,max_y                                                                                    
        do ii=min_x,max_x                                                                                
            Ujacob(ii,jj) = 1./60._dp                                                   &               
        &   *(-1._dp * (Uold(ii-2,jj) + Uold(ii+2,jj) + Uold(ii,jj-2) + Uold(ii,jj+2))  &               
        &   + 16._dp * (Uold(ii-1,jj) + Uold(ii+1,jj) + Uold(ii,jj-1) + Uold(ii,jj+1))  &               
        &   - 12._dp * h(1)**2 * Urhs(ii,jj))                                                            
        end do                                                                                           
    end do                                                                                               
end if

编辑：

按照 Wolfgang 的建议，我通过镜像边界旁边发生的事情来填充我的幽灵点。但是，这无济于事，代码仍然存在分歧。

1个回答

CFD 中常见的 Ghost 细胞处理

鬼细胞处理是 CFD 中的一门艺术，我们不能针对所有问题随意选择鬼细胞。让我们假设 $0$ 是边界点并且 $-1$ 是鬼点和 $u$ 是可变的，那么常用的幽灵细胞程序是，u(-1) =u(0)，这是某种意义上的等价 $\frac{du}{dx}=0$ 另一个过程是 u(-1) = -u(0)，如果我们希望壁面的法向速度为零 ( $v_n=0$ ) 常用于欧拉方程。

在这里，我将自己限制为一维热方程。从它背后的物理原理，我们可以将其扩展到其他椭圆方程。所有上述 BC 描述的热方程都可能失败。

如果我们强制 u(-1) =u(0) 相当于隔热墙 BC，但我们不确定墙是否隔热。在没有 BC 的情况下在两个点保持相同的温度是非物理的
如果我们使用 u(-1) = -u(0)，这也违反了物理，假设 u(0)=1000 K，u(-1)=-1000 K。它是物理的吗？-1000 K 不存在。即使存在，我们也绝不允许数值方案出现如此大的跳跃，除非有冲击。由于它是椭圆解，因此预计是平滑的。

稳态一维热方程为

\frac{d^{2} u}{d x^{2}} = 0

$\frac{d^2u}{dx^2} = 0$

2^{n d}

$2^{nd}$ 阶 FD 等价物是

\frac{d^{2} u}{d x^{2}} = (u_{i - 1} - 2 u_{i} + u_{i + 1}) / h^{2} = 0

$\frac{d^2u}{dx^2} = (u_{i-1}-2u_i+u_{i+1})/h^2 =0$ 然后

u_{i - 1} = 2 * u_{i} - u_{i + 1}

$u_{i-1}=2*u_i-u_{i+1}$ 对于 i=0 和 nx 这变成了这里

- 1

$-1$ 和

n x + 1

$nx+1$ 是鬼细胞点

u_{- 1} = 2 * u_{0} - u_{1}

$u_{-1}=2*u_0-u_{1}$

u_{n x + 1} = 2 * u_{n x} - u_{n x - 1}

$u_{nx+1}=2*u_{nx}-u_{nx-1}$

这些方程可用于幽灵单元点，请注意这些方程是二阶精度的，可能对解有影响。因为内部点是四阶方案。这里我附上了一个示例matlab代码，你可以玩一下。

clc
clear all
close all
nx=101;
x=linspace(0,1,nx);
dx=1/(nx-1);
u=0*x; %I.C
u_res=u;
u(2)=10; %bc because x(1) is ghost cell
u(nx-1)=10; %bc because x(nx) is ghost cell
itmax=10000;
for j1=1:1:itmax
for i1=1:nx
    if i1==1;
        u(i1)=u(2)*2-u(i1+2);  % Ghost cell values are framed 
                                    %from GDE in second order
    end

    if i1==nx;
        u(i1)=u(nx-1)*2-u(nx-2);
    end
    if i1==2
        u(i1)=10;
    end
    if i1==nx-1
        u(i1)=100;
    end
    if i1>2 && i1<nx-1
        u(i1)=(-u(i1-2)+16*u(i1-1)+16*u(i1+1)-u(i1+2))/30;       
    end
end
for i1=3:1:nx-2
    u_res(i1)=(-u(i1-2)+16*u(i1-1)+16*u(i1+1)-u(i1+2))/30-u(i1); %residue calc.
end
if max(abs(u_res))<0.0001 %convergence creteria
    break
end

end
plot(x(2:nx-1),u(2:nx-1))

如果你做以下事情，我会很高兴

您可以检查此问题的二阶幽灵单元处理对收敛的影响
你可以设置一个 $4^{th}$ 在代码中订购有偏差的幽灵单元插值方案并检查其效果。

其它你可能感兴趣的问题

上一篇显式和隐式预处理之间的区别下一篇计算科学硬件基准数据库