有限元法计算二维圆柱绕流问题——Python代码实现

百宝袋里是数学和代码

已于 2023-07-12 10:00:16 修改

阅读量1.5k

点赞数 5

分类专栏：计算流体力学文章标签： python

于 2023-05-11 02:07:33 首次发布

本文链接：https://blog.csdn.net/weixin_52046577/article/details/130612970

版权

计算流体力学专栏收录该内容

4 篇文章 1 订阅

订阅专栏

一、问题描述

选取流函数Ψ为变量，对拉普拉斯方程进行求解（右边界为自然边界条件，其余边界为本质边界条件）；

网格数据文件的生成暂时不在本文中详述。

二、节点和单元的数据读取

import numpy as np
import matplotlib as plt
from mpl_toolkits.mplot3d import Axes3D

# 打开文件
try:
    with open('grid.dat', 'r') as f:
        # 读取节点数NP和单元数NE
        line = f.readline().strip()
        NP, NE = map(int, line.split())

        # 读取节点坐标X[0:NP,0]和X[0:NP,1]
        X = np.zeros((NP, 2),dtype = float)
        for i in range(NP):
            line = f.readline().strip()
            X[i] = np.array(list(map(float, line.split())))

        # 读取单元节点对应关系数组NOD[0:NE,0:3]
        NOD = np.zeros((NE,3),dtype = int)
        for i in range(NE):
            line = f.readline().strip()
            NOD[i] = np.array(list(map(int, line.split())))-1 
            #注意文件里的索引是从1开始的，而节点坐标数组X是从0开始的

except FileNotFoundError:
    print("文件不存在")
except Exception as e:
    print("文件读取失败：", e)

三、设置本质边界条件

#设置本质边界条件
Ψ = np.zeros(NP)
substantial_bound_index = []
unknown_index = []
for i in range(NP):
    if X[i,1] == 0 :
        Ψ[i] = 0
        substantial_bound_index.append(i)
    elif X[i,1] == 2:
        Ψ[i] = 2
        substantial_bound_index.append(i)
    elif X[i,0] == -3.5:
        Ψ[i] = X[i,1]
        substantial_bound_index.append(i)
    elif abs(X[i,0]**2 + X[i,1]**2 - 1)<= 1e-3:#如果网格尺度变化，此处可能需要调整
        Ψ[i] =  0
        substantial_bound_index.append(i)
    else:
        unknown_index.append(i)

四、计算单元方程和整体方程

A_overall = np.zeros((NP,NP)) #整体方程的系数矩阵
f_overall = np.zeros(NP) #整体方程的右端项
#遍历所有单元，求单元方程的系数矩阵，并累加到整体方程的系数矩阵上
for i in range(NE):
    X_i = np.zeros(3)
    Y_i = np.zeros(3)
    node_index = np.zeros(3)
    node_index = NOD[i,:] #是从0开始的
    X_i = X[node_index,0]
    Y_i = X[node_index,1]
    
    A = 0.5*((X_i[1]-X_i[0])*(Y_i[2]-Y_i[0])-(Y_i[1]-Y_i[0])*(X_i[2]-X_i[0]))
    b1 = (Y_i[1]-Y_i[2])/(2*A)
    b2 = (Y_i[2]-Y_i[0])/(2*A)
    b3 = (Y_i[0]-Y_i[1])/(2*A)
    c1 = -(X_i[1]-X_i[2])/(2*A)
    c2 = -(X_i[2]-X_i[0])/(2*A)
    c3 = -(X_i[0]-X_i[1])/(2*A)
    
    A_overall[node_index[0],node_index[0]] += b1*b1 + c1*c1
    A_overall[node_index[1],node_index[1]] += b2*b2 + c2*c2
    A_overall[node_index[2],node_index[2]] += b3*b3 + c3*c3
    
    A_overall[node_index[0],node_index[1]] += b1*b2+c1*c2
    A_overall[node_index[1],node_index[0]] += b1*b2+c1*c2
    
    A_overall[node_index[0],node_index[2]] += b1*b3+c1*c3
    A_overall[node_index[2],node_index[0]] += b1*b3+c1*c3
    
    A_overall[node_index[1],node_index[2]] += b2*b3+c2*c3
    A_overall[node_index[2],node_index[1]] += b2*b3+c2*c3

#代入本质边界条件上的函数值，消元，只剩下待求节点的函数值
#计算消元后待求节点对应的右端项
for i in unknown_index:
    sum = 0
    for j in substantial_bound_index:
        sum += (-1)*A_overall[i,j]* Ψ[j]
    f_overall[i] = sum

#解线性方程组
A_tosolve = np.zeros((len(unknown_index),len(unknown_index)))
A_tosolve = A_overall[np.ix_(unknown_index,unknown_index)]
f_tosolve = f_overall[np.ix_(unknown_index)]
sol = np.linalg.solve(A_tosolve, f_tosolve)

#得到完整的节点Ψ数组
pos = 0
for i in unknown_index:
    Ψ[i] = sol[pos]
    pos+=1

五、计算单元和节点的流速和压强分布

首先利用前面计算出的节点流函数值，插值得到每个单元的流速；

然后，每个节点的流速则用相邻单元流速的加权平均（以单元面积为权重）得到；

最后，通过伯努利方程计算出节点的压强。

#计算单元的速度
vx = np.zeros(NE)
vy = np.zeros(NE)

#计算节点的速度——每个节点的速度采用相邻单元速度的面积加权平均
node_sum_area = np.zeros(NP) #每个节点相邻的累积面积
vx_node = np.zeros(NP)
vy_node = np.zeros(NP)

    #遍历所有单元
for i in range(NE):
    X_i = np.zeros(3)
    Y_i = np.zeros(3)
    node_index = np.zeros(3)
    node_index = NOD[i,:] #是从0开始
    X_i = X[node_index,0]
    Y_i = X[node_index,1]
    
    A = 0.5*((X_i[1]-X_i[0])*(Y_i[2]-Y_i[0])-(Y_i[1]-Y_i[0])*(X_i[2]-X_i[0]))
    b1 = (Y_i[1]-Y_i[2])/(2*A)
    b2 = (Y_i[2]-Y_i[0])/(2*A)
    b3 = (Y_i[0]-Y_i[1])/(2*A)
    c1 = -(X_i[1]-X_i[2])/(2*A)
    c2 = -(X_i[2]-X_i[0])/(2*A)
    c3 = -(X_i[0]-X_i[1])/(2*A)
    
    vx[i] = c1*Ψ[NOD[i,0]]+c2*Ψ[NOD[i,1]]+c3*Ψ[NOD[i,2]]
    vy[i] = -b1*Ψ[NOD[i,0]]-b2*Ψ[NOD[i,1]]-b3*Ψ[NOD[i,2]]

    #更新节点的速度均值
    for j in range(3):
        s0 = node_sum_area[NOD[i,j]]
        vx0 = vx_node[NOD[i,j]]
        vx_node[NOD[i,j]] = (s0*vx0+A*vx[i])/(s0+A)
        vy0 = vy_node[NOD[i,j]]
        vy_node[NOD[i,j]] = (s0*vy0+A*vy[i])/(s0+A)
        node_sum_area[NOD[i,j]]+= A
        
#计算节点的压力
p_node = np.zeros(NP)
for i in range(NP):
    p_node[i] = 0.5*(1-vx_node[i]**2 - vy_node[i]**2)

六、计算结果

一、流函数数值解的3D图

from mpl_toolkits.mplot3d import Axes3D
from matplotlib.tri import Triangulation

tri = Triangulation(X[:, 0], X[:, 1])
#一、绘制流函数数值解的3D图
fig = plt.pyplot.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(X[:, 0], X[:, 1], Ψ, triangles=tri.triangles, cmap='viridis',alpha = 0.9)
ax.set_title('Ψ_num')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Ψ')
plt.pyplot.show()

二、数值解的流线图

#二、绘制数值解的流线图
z = Ψ
# 绘制等值线
levels = np.linspace(z.min(), z.max(), 30)
plt.pyplot.tricontour(tri, z, levels=levels, colors='k')
# 添加等值线标签
plt.pyplot.tricontourf(tri, z, levels=levels, cmap='viridis')
plt.pyplot.colorbar()
plt.pyplot.title('streamline _num')
plt.pyplot.show()

对比：解析解的流线图

# 三、绘制解析解的流线图
#解析解
Ψ_true= np.zeros(NP)
for i in range(NP):
    Ψ_true[i] = X[i, 1]*(1-1/(X[i,0]**2 + X[i,1]**2))
z = Ψ_true
# 绘制等值线
levels = np.linspace(z.min(), z.max(), 30)
plt.pyplot.tricontour(tri, z, levels=levels, colors='k')
# 添加等值线标签
plt.pyplot.tricontourf(tri, z, levels=levels, cmap='viridis')
plt.pyplot.colorbar()
plt.pyplot.title('streamline _true')
plt.pyplot.show()

可以观察到数值解与解析解的流线形状是比较相近的。

三、节点流速和压强的数值结果

#四、绘制节点压强的3D图
fig = plt.pyplot.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(X[:, 0], X[:, 1], p_node, triangles=tri.triangles, cmap='viridis',alpha = 0.9)
ax.set_title('pressure')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('p')
plt.pyplot.show()

#五、绘制节点流速的3D图
fig = plt.pyplot.figure()

ax = fig.add_subplot(121, projection='3d')
ax.plot_trisurf(X[:, 0], X[:, 1], vx_node, triangles=tri.triangles, cmap='viridis',alpha = 0.9)
ax.set_title('vx')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('vx')

ax = fig.add_subplot(122, projection='3d')
ax.plot_trisurf(X[:, 0], X[:, 1], vy_node, triangles=tri.triangles, cmap='viridis',alpha = 0.9)
ax.set_title('vy')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('vy')

plt.pyplot.show()

数值解在圆柱壁面附近出现较大的波动。