CFD——非结构网格梯度计算(不修正)

发布时间 2023-07-24 17:01:28作者: 摩天仑

image
本案例在计算非结构网格的梯度时,不使用修正方法。将直接使用f'处的∅值

计算流程
image

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]

# C单元与周围单元质心连线,CF矢量
CF_List = np.random.rand(N, 2)
for i in range(N):
  CF_List[i] = P_F_List[i] - P_C

# 计算两直线的交点函数
def intersection_point(point1, point2,point3,point4):
  n1n2 = point2 - point1
  CF1 = point4 - point3
  n1n2_k = n1n2[1]/n1n2[0]
  n1n2_b = point1[1] - n1n2_k * point1[0]
  CF1_k = CF1[1]/CF1[0]
  CF1_b = point3[1] - CF1_k * point3[0]
  x = (n1n2_b - CF1_b) / ( CF1_k - n1n2_k)
  y = n1n2_k * x + n1n2_b
  return np.array([x,y])

# 计算两点之间的距离
def distance_two_points(point1,point2):
  dis = ((point1[0]-point2[0])**2+(point1[0]-point2[0])**2)**0.5
  return dis

# 两直线的交点
P_nn_CF_List = np.random.rand(N, 2)
for i in range(N):
  if i < N-1:
    P_nn_CF_List[i] = intersection_point(P_n_List[i],P_n_List[i+1],P_C,P_F_List[i])
  else:
    P_nn_CF_List[i] = intersection_point(P_n_List[i],P_n_List[0],P_C,P_F_List[i])

# gc比例
gc_CF_List = np.random.rand(N)
for i in range(N):
  gc_CF_List[i] = distance_two_points(P_nn_CF_List[i],P_F_List[i])/distance_two_points(P_C,P_F_List[i])

# 顶点连线与质心连线交点的phi值
phi_f_sk_List = np.random.rand(N)
for i in range(N):
  phi_f_sk_List[i] = gc_CF_List[i] * phi_C + (1 - gc_CF_List[i]) * phi_F_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]

print(Sf_List)
print(phi_f_sk_List)

# 计算单元C的梯度,gra_phi_C
gra_phi_C = np.random.rand(2)
sumi = 0
sumj = 0
for i in range(N):
  sumi += phi_f_sk_List[i] * Sf_List[i][0]
  sumj += phi_f_sk_List[i] * Sf_List[i][1]
gra_phi_C[0] = sumi/V_C
gra_phi_C[1] = sumj/V_C