CFD——使用扩展法计算网格梯度

求解流程

在高斯格林公式中，需要用到phi_f，以下是求解phi_f1的步骤（这里只给出phi_f1）

求解

代码实现

# 使用扩展法计算

import numpy as np

# ----------------------------初始化条件---------------------------#
N = 6

# 各质心的位置坐标
P_C = np.array([13,11])
P_F_List = np.array([[4.5,9.5],[8,3],[17,3.5],[22,10],[16,20],[7,18]])

# C单元各个顶点坐标
P_n_List = np.array([[9,14],[8,8],[12,5],[17,9],[17.5,14],[12,17]])

# 各单元值
phi_C = 167
phi_F_List =  np.array([56.75,35,80,252,356,151])

# 各单元的梯度
gra_phi_F_List = np.array([[10.5,5.5],[4,9],[4.5,18],[115,23],[21,17],[19,8]])

# C单元体积
V_C = 76

# ---------------------------计算面法向--------------------------#
# C单元顶点连线nn矢量
nn_List = np.random.rand(N, 2)
for i in range(N):
  if i < 5:
    nn_List[i] = P_n_List[i+1] - P_n_List[i]
  else:
    nn_List[i] = P_n_List[0] - P_n_List[i]

# 计算面法向
Sf_List = np.random.rand(N, 2)
for i in range(N):
  Sf_List[i][0] = nn_List[i][1]
  Sf_List[i][1] = -nn_List[i][0]

# -----------------------------使用紧凑法------------------------#
# 计算phi_f1
# dis_dF
dis_n1F6 = ((P_n_List[0][0] - P_F_List[5][0]) ** 2 + (P_n_List[0][1] - P_F_List[5][1]) ** 2) ** 0.5
dis_n1F1 = ((P_n_List[0][0] - P_F_List[0][0]) ** 2 + (P_n_List[0][1] - P_F_List[0][1]) ** 2) ** 0.5
dis_n1C = ((P_n_List[0][0] - P_C[0]) ** 2 + (P_n_List[0][1] - P_C[1]) ** 2) ** 0.5
dis_n2F1 = ((P_n_List[1][0] - P_F_List[0][0]) ** 2 + (P_n_List[1][1] - P_F_List[0][1]) ** 2) ** 0.5
dis_n2F2 = ((P_n_List[1][0] - P_F_List[1][0]) ** 2 + (P_n_List[1][1] - P_F_List[1][1]) ** 2) ** 0.5
dis_n2C = ((P_n_List[1][0] - P_C[0]) ** 2 + (P_n_List[1][1] - P_C[1]) ** 2) ** 0.5
print(dis_n1F6)
print(dis_n1F1)
print(dis_n1C)
print(dis_n2F1)
print(dis_n2F2)
print(dis_n2C)

phi_n1 = (phi_F_List[5]/dis_n1F6 + phi_F_List[0]/dis_n1F1 + phi_C/dis_n1C) / (1/dis_n1F6 + 1/dis_n1F1 + 1/dis_n1C)
print(phi_n1)
phi_n2 = (phi_F_List[0]/dis_n2F1 + phi_F_List[1]/dis_n2F2 + phi_C/dis_n2C) / (1/dis_n2F1 + 1/dis_n2F2 + 1/dis_n2C)
print(phi_n2)
phi_f1 = 0.5 * (phi_n1 + phi_n2)