Julia CFD|12 二維泊松方程

Poisson方程的表達式為：

與Laplace方程不同，Poisson方程帶有源項。

1 方程離散

Poisson的離散方式與Laplace方程類似：

改成迭代式形式為：

2 初始值與邊界值

假設計算區域初始條件下

對於源項

3 計算代碼

using PyPlot
matplotlib.use("TkAgg")
 
nx = 50
ny = 50
nt = 100
xmin = 0
xmax = 2
ymin = 0
ymax = 1
 
dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)
 
 
# 初始化
p = zeros((ny, nx));
b = zeros((ny, nx));
x = range(xmin, xmax, length = nx);
y = range(xmin, xmax, length = ny);
 
# 源項
b[floor(Int, ny / 4), floor(Int, nx / 4)] = 100;
b[floor(Int, 3 * ny / 4), floor(Int, 3 * nx / 4)] = -100;
 
function plot2D(x, y, p)
    fig = PyPlot.figure(figsize=(11,7), dpi=100)
    ss = PyPlot.surf(x,y,p, rstride=1, cstride=1, cmap=ColorMap("coolwarm"),
            linewidth=0)
    xlim(0,2)
    ylim(0,1)
    PyPlot.show()
end
查看初始值分布。
plot2D(x,y,p)
初始分布如下圖所示。
迭代計算代碼如下。
for it in 1:nt
    pd = copy(p)
    p[2:end-1,2:end-1] = ((pd[2:end-1,3:end] + pd[2:end-1,1:end-2])*dy^2 +
        (pd[3:end,2:end-1]+pd[1:end-2,2:end-1])*dx^2 -b[2:end-1,2:end-1]*dx^2*dy^2)/(2*(dx^2+dy^2))
 
    p[1,:] .= 0
    p[ny,:] .= 0
    p[:,1] .= 0
    p[:,nx] .= 0
end
計算結果如下圖所示。
完整代碼如下。
using PyPlot
matplotlib.use("TkAgg")
 
nx = 50
ny = 50
nt = 100
xmin = 0
xmax = 2
ymin = 0
ymax = 1
 
dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)
 
# 初始化
p = zeros((ny, nx));
b = zeros((ny, nx));
x = range(xmin, xmax, length = nx);
y = range(xmin, xmax, length = ny);
 
# 源項
b[floor(Int, ny / 4), floor(Int, nx / 4)] = 100;
b[floor(Int, 3 * ny / 4), floor(Int, 3 * nx / 4)] = -100;
plot2D(x, y, p)
 
for it in 1:nt
    pd = copy(p)
    p[2:end-1,2:end-1] = ((pd[2:end-1,3:end] + pd[2:end-1,1:end-2])*dy^2 +
        (pd[3:end,2:end-1]+pd[1:end-2,2:end-1])*dx^2 -b[2:end-1,2:end-1]*dx^2*dy^2)/(2*(dx^2+dy^2))
 
    p[1,:] .= 0
    p[ny,:] .= 0
    p[:,1] .= 0
    p[:,nx] .= 0
end
 
function plot2D(x, y, p)
    fig = PyPlot.figure(figsize=(11,7), dpi=100)
    ss = PyPlot.surf(x,y,p, rstride=1, cstride=1, cmap=ColorMap("coolwarm"),
            linewidth=0)
    xlim(0,2)
    ylim(0,1)
    PyPlot.show()
end
 
plot2D(x, y, p)