【CFD算法】伪势LBM模拟壁面接触角（python）

伪势模型（LBGK-EDM）模拟壁面接触角（Python ）

最近学习 LBM，把大佬给的 MATLAB 代码用 Python 重写了一遍。可以实现不同润湿性表面的接触角模拟以及粗糙表面接触角的模拟计算，采用的是伪势模型。

以下是实现效果：

代码实现

# -*- coding:UTF-8 -*-
# Author:DABAIMENG
# Time:2023/7/14 8:00
from numpy import *
from matplotlib import cm
import matplotlib.pyplot as plt

np.set_printoptions(threshold=np.inf)

基本参数设置

# LB 基本参数
G = -5.0
rho_liq = 1.95    # 液体密度
rho_gas = 0.15    # 气体密度
nx = 200          # 定义计算尺寸
ny = 200
noput = 1

# 壁面密度，可调整 0.80 模拟不同亲疏水性
rho_boundary = rho_gas + 0.80 * (rho_liq - rho_gas)

# 粒子分布函数、平衡分布函数、速度分量分配内存
f = zeros((9, nx, ny))
fx = zeros((9, nx, ny))
feq = zeros((9, nx, ny))
feq1 = zeros((9, nx, ny))
u = zeros((nx, ny))
v = zeros((nx, ny))

# 修正速度分配内存
uf = zeros((nx, ny))
vf = zeros((nx, ny))
psi = zeros((nx, ny))       # 势函数分配内存
forcx = zeros((nx, ny))     # 外力项分配内存
forcy = zeros((nx, ny))
rho = rho_gas * ones((nx, ny))  # 初始化密度为液体密度
pressure = zeros((nx, ny))
x = zeros(nx)
y = zeros(ny)

D2Q9 模型参数设置

# D2Q9 模型参数设置
w = [1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36, 4/9]
opp = [2, 3, 0, 1, 6, 7, 4, 5, 8]
cx = [1, 0, -1, 0, 1, -1, -1, 1, 0]
cy = [0, 1, 0, -1, 1, 1, -1, -1, 0]
c2 = 1/3
omega = 1

# 液滴半径
Rd = 42
err = 10
counter = 1

壁面障碍物设定

# 壁面障碍物设定
height_barrier = 10
obst = zeros((nx, ny))
obstk = zeros((9, nx, ny))
a = 1

for i in range(nx):
    if mod(i-5, 23) < 3:
        obst[i-1, 0:height_barrier] = a
        obst[i-1, ny-height_barrier:ny] = a

obst[0:10, 0:height_barrier] = a
obst[0:10, ny-height_barrier:ny] = a
obst[nx-12:nx, 0:height_barrier] = a
obst[nx-12:nx, ny-height_barrier:ny] = a

for i in range(nx):
    for j in range(ny):
        for k in range(8):
            if j + cy[k] < 0 or j + cy[k] > ny-1:
                obstk[k, i, j] = 1
            elif i + cx[k] < 0 or i + cx[k] > nx-1:
                obstk[k, i, j] = 2
            elif obst[i + cx[k], j + cy[k]] == 1:
                obstk[k, i, j] = 1

初始化

# 初始化液滴
for i in range(nx):
    for j in range(ny):
        if (i - nx/2)*(i - nx/2) + (j - Rd - 7)*(j - Rd - 7) < Rd*Rd:
            rho[i, j] = rho_liq
        else:
            rho[i, j] = rho_gas

# 计算初始压力
for j in range(ny):
    for i in range(nx):
        psi[i, j] = 1 - exp(-rho[i, j])
pressure = rho/3 + G * psi * psi/6

碰撞函数

def collision(nx, ny, u, v, cx, cy, omega, f, rho, w, forcx, forcy):
    for i in range(nx):
        for j in range(ny):
            t1 = u[i, j] * u[i, j] + v[i, j] * v[i, j]
            t3 = (u[i, j] + omega * forcx[i, j] / rho[i, j]) * (u[i, j] + omega * forcx[i, j] / rho[i, j]) + (
                  v[i, j] + omega * forcy[i, j] / rho[i, j]) * (v[i, j] + omega * forcy[i, j] / rho[i, j])
            for k in range(9):
                t2 = u[i, j] * cx[k] + v[i, j] * cy[k]
                t4 = (u[i, j] + omega * forcx[i, j] / rho[i, j]) * cx[k] + (
                      v[i, j] + omega * forcy[i, j] / rho[i, j]) * cy[k]
                feq[k, i, j] = rho[i, j] * w[k] * (1.0 + 3.0 * t2 + 4.5 * t2 - 1.5 * t1)
                feq1[k, i, j] = rho[i, j] * w[k] * (1.0 + 3.0 * t4 + 4.5 * t4 - 1.5 * t3)
                f[k, i, j] = (1 - omega) * f[k, i, j] - (1 - omega) * feq[k, i, j] + feq1[k, i, j]
    return f

外力计算函数

def force(nx, ny, cx, cy, rho, rho_boundary, rho_liq, rho_gas, w, G, obstk):
    for j in range(nx):
        for i in range(ny):
            forcxs = 0.0
            forcys = 0.0
            for k in range(8):
                newx = i + cx[k]
                newy = j + cy[k]
                if obstk[k, i, j] == 1:
                    if rho[i, j] > rho_liq / 8:
                        psi[i, j] = 1 - exp(-rho_boundary)
                        forcxs = forcxs - w[k] * G * psi[i, j] * cx[k]
                        forcys = forcys - w[k] * G * psi[i, j] * cy[k]
                    else:
                        psi[i, j] = 1 - exp(-rho_gas)
                        forcxs = forcxs - w[k] * G * psi[i, j] * cx[k]
                        forcys = forcys - w[k] * G * psi[i, j] * cy[k]
                elif obstk[k, i, j] == 2:
                    psi[i, j] = 1 - exp(-rho[i, j])
                    forcxs = forcxs - w[k] * G * psi[i, j] * cx[k]
                    forcys = forcys - w[k] * G * psi[i, j] * cy[k]
                else:
                    psi[i, j] = 1 - exp(-rho[newx, newy])
                    forcxs = forcxs - w[k] * G * psi[i, j] * cx[k]
                    forcys = forcys - w[k] * G * psi[i, j] * cy[k]
            forcx[i, j] = (1 - exp(-rho[i, j])) * forcxs
            forcy[i, j] = (1 - exp(-rho[i, j])) * forcys
    return forcx, forcy, psi

宏观参数计算函数

def ruv(u, v, f, psi, rho, obst):
    for i in range(nx):
        for j in range(ny):
            if obst[i, j] == 0:
                rho[i, j] = sum(f[:, i, j], axis=0)
                u[i, j] = (sum(f[[0, 4, 7], i, j], axis=0) - sum(f[[2, 5, 6], i, j], axis=0)) / rho[i, j]
                v[i, j] = (sum(f[[1, 4, 5], i, j], axis=0) - sum(f[[3, 6, 7], i, j], axis=0)) / rho[i, j]
                pressure = rho[i, j] / 3 + 1.5 * G * psi[i, j] * psi[i, j]
            else:
                rho[i, j] = nan
                u[i, j] = nan
                v[i, j] = nan
    return rho, u, v, pressure

边界条件

def boundary(nx, ny, f, f_old, obstk, opp, obst):
    fx = np.array(f, copy=True)
    for i in range(nx):
        for j in range(ny):
            for k in range(8):
                if obstk[k, i, j] == 1:
                    fx[opp[k], i, j] = f_old[k, i, j]
    fx1 = np.array(fx, copy=True)
    fx2 = np.array(fx, copy=True)
    for j in range(ny):
        if obst[1, j] == 0:
            fx2[0, 0, j] = fx1[0, nx-1, j]
            fx2[4, 0, j] = fx1[4, nx-1, j]
            fx2[7, 0, j] = fx1[7, nx-1, j]
            fx2[2, nx-1, j] = fx1[2, 0, j]
            fx2[5, nx-1, j] = fx1[5, 0, j]
            fx2[6, nx-1, j] = fx1[6, 0, j]
    return fx2, fx1

迁移项

def stream(f):
    f_old = np.array(f, copy=True)
    for i in range(9):
        f[i, :, :] = roll(roll(f[i, :, :], cx[i], axis=0), cy[i], axis=1)
    return f, f_old

def uvf(u, v, forcx, forcy, rho):
    for i in range(nx):
        for j in range(ny):
            uf[i, j] = u[i, j] + omega * forcx[i, j] / (2 * rho[i, j])
            vf[i, j] = v[i, j] + omega * forcy[i, j] / (2 * rho[i, j])
    return uf, vf

主循环

counter = 1
while counter < 1000:
    err = 0
    forcx, forcy, psi = force(nx, ny, cx, cy, rho, rho_boundary, rho_liq, rho_gas, w, G, obstk)
    print('碰撞%s' % counter)
    f = collision(nx, ny, u, v, cx, cy, omega, f, rho, w, forcx, forcy)
    print('迁移%s' % counter)
    f, f_old = stream(f)
    f, fill = boundary(nx, ny, f, f_old, obstk, opp, obst)
    print('宏观量计算%s' % counter)
    rho, u, v, pressure = ruv(u, v, f, psi, rho, obst)
    print('速度修正%s' % counter)
    uf, vf = uvf(u, v, forcx, forcy, rho)
    print('***********')
    rho_tot = 0
    for i in range(nx):
        for j in range(ny):
            if rho[i, j] != '':
                rho_tot = rho_tot + rho[i, j]
    print(rho_tot)
    plt.imshow(rho, cmap=cm.bwr)
    plt.savefig(r"./计算结果/" + str(counter) + ".png")
    plt.clf()
    counter = counter + 1