有限元法计算二维圆柱绕流问题——无粘有环量的情况

百宝袋里是数学和代码

已于 2023-07-12 10:01:22 修改

阅读量2.6k

点赞数 5

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

于 2023-06-29 05:12:46 首次发布

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

版权

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

4 篇文章

订阅专栏

文章详细介绍了使用有限元法解决二维有环量圆柱绕流问题的步骤，包括网格划分、边界条件设置、求解流程及误差分析。通过改变环量大小和外边界距离，探讨了这些因素对计算结果的影响，发现误差主要来源于边界条件和网格尺度。此外，还展示了计算得到的流线图、流速和压强分布，并对比了数值解与解析解的差异。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

写在前面

该系列文章是笔者的学习记录，旨在通过分享促进交流加深理解，文中难免有许多纰漏与理解不当之处，若能批评指正必定感激不尽。

侧重分享经典CFD方法的具体实现结果，详细的理论推导请参考教材。

一、问题描述

与上一篇有限元法计算二维圆柱绕流问题——Python代码实现不同的是，本文考虑有环量圆柱绕流问题的求解，这同样是一个有精确解的问题。加入了平面点涡流动以后，圆柱绕流问题表现出更贴近现实的升力的现象。

取一个充分大的区间，在外边界提无穷远边界条件。下面将逐步讲解用有限元法计算的过程，并针对数值结果考察外边界距离以及网格尺度对误差的影响。

回顾-有限元法简介：

有限元法是求解计算力学问题的重要方法。首先需要将求解区域进行离散化，划分为许多几何形状简单规则的单元，如二维的三角形或四边形，三维的四面体或六面体；然后，在每个单元内，用一个比较简单的解析函数来逼近微分方程的解，此函数在单元内用一组选定的单元基函数的线性组合表示，而其中的系数通常是节点上的物理量，它是待求解的。这样，每个单元只要有适当数量的节点物理量的值，就可以满足对插值函数的光滑性和精度的要求；接下来，在满足微分方程和相应的初边值条件下，对全部子域进行积分：对每个单元分别进行积分，形成“单元方程”，通过总体合成，得到总体有限元方程组；最后，用适当的方法解方程组，可得节点物理量的值，进而可表示出各单元内的近似解。

二、网格划分

计算区域：由于有环量的圆柱绕流问题左右对称（阻力为0）而上下不对称（升力不为0），因此计算圆柱左边的绕流区域即可。为了便于施加边界条件，体现环量的作用，将左上的1/4区域和左下的1/4区域分别划分网格计算（网格的生成直接利用对称性）。

本文编写了一份朴素的代码将待求解区域划分为三角形网格。总体来讲，首先分区并在每个子区域内均匀地取网格点，然后借助Header only constrained delaunay tessellator的C++代码中的三角剖分方法生成网格文件。示例网格如下图所示（分2层）

但由于在加大网格密度时计算效率不高，简便起见，利用Gmsh来生成误差分析所需的均匀三角网格。所有用到的网格已放在文件中。

三、有限元求解流程（总体与上一次相似）

1. 节点和单元的数据读取

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

如果是读取Gmsh生成的msh文件：

import meshio

# 读取 MSH 文件
mesh = meshio.read("D:/area_24.0_14.4_gird2.msh")

# 提取网格数据
points = mesh.points  # 点坐标
cells = mesh.cells    # 单元数据
# 还可以提取其他的网格属性，如单元类型、单元标签等

NP = points.shape[0]
NE = [cell_data for cell_data in cells if cell_data.type == 'triangle'][0].data.shape[0]

# 读取节点坐标X[0:NP,0]和X[0:NP,1]
X = np.zeros((NP, 2),dtype = float)
X = points[:,0:2]
X2 = np.zeros((NP, 2),dtype = float)
X2[:,0] = points[:,0]
X2[:,1] = - points[:,1]
# 读取单元节点对应关系数组NOD[0:NE,0:3]
NOD = np.zeros((NE,3),dtype = int)
NOD = [cell_data for cell_data in cells if cell_data.type == 'triangle'][0].data
#注意文件里的索引是从1开始的，而节点坐标数组X是从0开始的

如果是读取自己的网格程序生成的dat文件：

# 打开文件
try:
    with open('cylinder.dat', 'r') as f:
        # 读取节点数NP和单元数NE（自己生成的网格文件有些问题 实际的单元数比文件开头写的数量少）
        line = f.readline().strip()
        NP, NE_max = 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_max,3),dtype = int)
        NE=0 #记录一下真实存在的单元数量
        for i in range(NE_max):
            line = f.readline().strip()
 
            if line:#读到文件末尾就结束，针对不知道单元确切数量的情况
                NOD[i] = np.array(list(map(int, line.split())))-1 
                NE+=1
            else:
                break
            #注意文件里的索引是从1开始的，而节点坐标数组X是从0开始的

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

NOD = NOD[:NE, :]

2. 设置边界条件

#几何和流动参数
gamma = 5 #环量大小
r = 1
Rx= X[:,0].max()- X[:,0].min()
Ry = X[:,1].max()- X[:,1].min()

我们通过施加不同的边界条件来体现出有环量与无环量的不同。有环量的情况相当于绕圆柱额外施加一个环量为Gamma的平面点涡流动，易得其流函数为 $-\frac{\Gamma}{2\pi}ln(R)$ 。

根据圆柱绕流问题的定义，圆柱表面为流线，不妨设bc边界上 $\psi_{bc} = -\frac{\Gamma}{2\pi}ln(r)$ ；
ba边界上流函数不再是常数，设置ba上 $\psi_{ba} = -\frac{\Gamma}{2\pi}ln(R)$ ，R是到圆心的距离；
根据无穷远均匀来流（速度为1）的边界条件，ae边界与流线垂直， $\psi_{ae} = -\frac{\Gamma}{2\pi}ln(Rx) + y$ ；
de边界沿板仍为流线（外边界距离足够远时近似满足）， $\psi_{ed}= \psi_{e} = -\frac{\Gamma}{2\pi}ln(Rx) + Ry$ ；
增加的平面点涡流动仍然与cd边界垂直，无切向速度，仍然满足 $\frac{\partial \psi_{cd}}{\partial \mathbf{n}} =0$ 的自然边界条件。

#设置本质边界条件
Ψ = np.zeros(NP)
substantial_bound_index = []
unknown_index = []
for i in range(NP):
    if X[i,1] == 0 :
        Ψ[i] = - gamma*np.log(abs(X[i,0]))/(2*math.pi)
        substantial_bound_index.append(i)
    elif X[i,1] == Ry:
        Ψ[i] = Ry - gamma*np.log(abs(Rx))/(2*math.pi)
        substantial_bound_index.append(i)
    elif X[i,0] == -Rx:
        Ψ[i] = X[i,1] - gamma*np.log(Rx)/(2*math.pi)
        substantial_bound_index.append(i)
    elif abs(X[i,0]**2 + X[i,1]**2 - r*r)<= 1e-3:#如果网格尺度变化可能需要调整判断标准
        Ψ[i] = - gamma*np.log(r)/(2*math.pi)
        substantial_bound_index.append(i)
    else:
        unknown_index.append(i)

#验证本质边界条件的设置是否正确（主要是圆弧上的）
for i in substantial_bound_index:
    if(X[i,1]<=3 and X[i,0]>=-4):
        plt.scatter(X[i,0],X[i,1],alpha = 0.5,s=2)
plt.show()

3. 计算单元方程和整体方程

为了在网格加密时保持可接受的计算速度和存储效率，改用稀疏矩阵的存储方式。可处理的节点数量能达到数十万，解方程组的速度明显加快。

import scipy.sparse as sp
import scipy.sparse.linalg as splinalg

A_overall = sp.lil_matrix((NP,NP), dtype=float)
f_overall = np.zeros(NP) #整体方程的右端项，由于切向速度为0所以本来为0

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]
    # print(X[2,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]))
    # print(i,A)
    #由于自己的网格文件有些问题，存在三点共线的情况，直接忽略掉
    （有限元最终的方程组是基于节点而非单元的，没有影响）
    if A==0:
        continue
    
    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_f = 0
    for j in substantial_bound_index:
        sum_f += (-1)*A_overall[i,j]* Ψ[j]
    f_overall[i] = sum_f
#注：移项前的右端项为0，因为自然边界上切向速度为0

#首先得到待求解unknown_index节点变量所对应的A_tosolve系数矩阵（这里使用了一些稀疏矩阵的访问技巧）
import bisect
A_overall_coo = A_overall.tocoo()
A_tosolve = sp.lil_matrix((len(unknown_index), len(unknown_index)), dtype=float)

for i, j, val in zip(A_overall_coo.row, A_overall_coo.col, A_overall_coo.data): 
    #由于unknown_index是从小到大有序的，使用二分查找来加快速度
    idx_i = bisect.bisect_left(unknown_index, i)
    idx_j = bisect.bisect_left(unknown_index, j)

    if idx_i < len(unknown_index) and unknown_index[idx_i] == i and idx_j < len(unknown_index) and unknown_index[idx_j] == j:
        A_tosolve[idx_i, idx_j] = val


f_tosolve = f_overall[np.ix_(unknown_index)]
csc_A_tosolve = A_tosolve.tocsc()
# 在 CSC 矩阵上执行线性方程组求解，效率更高
sol = splinalg.spsolve(csc_A_tosolve, f_tosolve)

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

4. 计算单元和节点的流速和压强分布

首先利用前面计算出的节点流函数值，插值得到每个单元的流速；
然后，每个节点的流速则用相邻单元流速的加权平均（以单元面积为权重）得到；
最后，通过伯努利方程计算出节点的压强。

#计算单元的速度
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]))
    if A==0: #vx vy数组的中间会有一些单元是空的没赋值，但是没关系
        continue
    
    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)

5. 绘制流线图、流速和压强

先利用我们自己的数据来创建绘图网格（由于自己的网格有点问题，必须筛选出 valid_triangles）

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


def are_points_collinear(p1, p2, p3):
    """
    检查三个点是否共线
    """
    return abs(p1[0] * (p2[1] - p3[1]) + p2[0] * (p3[1] - p1[1]) + p3[0] * (p1[1] - p2[1])) <= 1e-16

# 用于存储有效的三角形索引
valid_triangles = []

# 检查三角形的拓扑关系数组 triangles
for triangle in NOD:
    p1, p2, p3 = X[triangle[0]], X[triangle[1]], X[triangle[2]]
    if not are_points_collinear(p1, p2, p3):
        # 顶点不共线，将三角形索引添加到有效的三角形数组中
        valid_triangles.append(triangle)


tri = Triangulation(X[:, 0], X[:, 1],triangles=valid_triangles)

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

#解析解的流线图
for i in range(NP):
    Ψ_true[i] = X[i, 1]*(1-1/(X[i,0]**2 + X[i,1]**2)) - gamma*np.log(X[i,0]**2 + X[i,1]**2)/(2*np.pi)
#其余都同上