0. Introduction

计算流体力学基础
大作业说明文档

学号:

姓名:

Tel:

xxxx年xx月

[TOC]

0. Introduction

本文档为计算流体力学基础大大作业的说明文档

使用语言：Python 3.7版本；
本作业中用到的额外库：numpy (用于矩阵运算)， pandas (读取csv文件用)，matplotlib （用于绘制3D图片）
Python额外库安装方法python -m pip install --user numpy matplotlib pandas
可选库：Tensorflow(C语言级别的矩阵计算速度，为机器学习、科学计算优化，比numpy的矩阵运算快5倍，矩阵较大时还可使用GPU计算)，可取代numpy。
CPU配置：Intel(R) Xeon(R) CPU E5-2637 v4 @ 3.50GHZ 3.50GHZ (2 处理器)，计算时间供参考。

本次大作业分为三个，其文件夹格式为：

O网格生成
1. 输入：input.txt
2. 运行文件：ogrid.py
3. 输出：Result/Cp-x.dat, Result/flowfield.dat, Result/o-grid.dat
一维Burgers方程数值求解
1. 输入：input.txt
2. 运行文件：burgers.py
3. 输出：2D文件夹中的mu=0.0.dat, mu=0.01.dat, mu=0.05.dat
4. 可选输出：3D文件夹中的三个三维图，mu=0.0.png, mu=0.01.png, mu=0.05.png
准一维Laval流动
1. 输入：input.txt
2. 运行文件：laval.py
3. 输出：Result/Cx=0.05, 0.2, 0.15.dat

1. O网格生成与数值求解

1.1. 输入文档说明

输入文档格式如下：

Points,Radius,iterlimit,errlimit,Velocity 
101,20,20000,1e-05,20.0

Points表示翼型取点个数N，其中后缘点算作两个点；
Radius表示远场边界条件半径r；
iterlimit表示最大容许的迭代次数；
errlimit表示最大误差限，即迭代前后误差小于errlimit时停止迭代；
Velocity表示最大来流速度$V_\infty$。

1.2. 运行文档ogrid.py说明

文档定义两个类：

class Ogrid():
    # 输入参数：迭代次数N, 半径radius, 迭代和误差限
    def __init__(self, N, radius, iterlimit, errlimit)
    # 初始化网格函数
    def initialGrid(self)
    # 网格迭代生成函数
    def gridGenerate(self)
    # 网格导出为dat文件函数
    def gridExport(self, x, y)
class GridSolver():
    # 输入参数：网格点坐标x, y, 自由来流速度Vf, 迭代和误差限
    def __init__(self, x, y, Vf, errlimit, iterlimit)
    # 流函数迭代求解器
    def phiSolver(self)
    # 流场参数求解器
    def flowSolver(self)
    # 流畅参数x, y, u, v, Cp导出函数
    def flowExport(self, x, y, u, v, Cp)

迭代中差分格式数学表达式 $$ u_i^j=u_{i+1}^j-u_{i-1}^j $$ 差分格式python代码实现

u[:, 1:-1]=u[:, 2:] - u[:, 0:-2]

main函数实现

if __name__ == '__main__':
	input = pd.read_csv('input.txt', sep=',', usecols=['Points', 'Radius', \
        'iterlimit','errlimit', 'Velocity'])
	N, r, iterlimit, errlimit, Vf = input['Points'][0], input['Radius'][0], \
    input['iterlimit'][0], input['errlimit'][0], input['Velocity'][0]

	ogrid = Ogrid(N, r, iterlimit, errlimit)
	x, y = ogrid.gridGenerate()
	ogrid.gridExport(x, y)
	flow = GridSolver(x, y, Vf, errlimit, iterlimit)
	u, v, Cp, phi = flow.flowSolver()
	flow.flowExport(x, y, u, v, Cp)

参考运行时间：18.0 s

1.3. 运行结果展示

1.3.1. o网格生成图片

Fig 1. 翼型后缘网格和全局网格

1.3.2. 流场显示

Fig 2. 压力系数云图与流线图

Fig 3. x, y方向速度云图

1.3.3. 翼型表面压力系数曲线

Fig 4. 翼型表面压力系数曲线

2. 一维Burgers方程数值求解

2.1. 输入文档说明

输入文档格式如下：

time,mu
0.5,0.0
1.2,0.01
2.0,0.05

第一列是时间，第二列是粘性系数，可以以$n\times n$的形式输入，即输入$n$个时间，需对应输入$n$个粘性系数。

2.2. 运行文档burgers.py说明

该运行文档定义一个类：

class Burgers:
    # 输入粘性系数，最后的计算停止时间，网格空间划分尺度
    def __init__(self, visco, time, N)
    # 对空间和时间进行离散化
    def discretization(self)
    # 使用MacCormack格式进行迭代计算
    def macCormack(self, u)

main函数实现：

if __name__ == '__main__':
	N = 200
	seper = os.sep 		# Use this seperator to satisfy both Windows and Linux usage

	# Read input file of time and viscosity
	time_visco = pd.read_csv('input.txt', sep = ',', usecols=['time', 'mu'])
	time = max(time_visco['time'])
	visco = time_visco['mu']

	# Compute Burgers equation at different viscosity
	for i in range(0, len(visco)):
		u = np.zeros(N + 1)
		x = np.zeros(N + 1)

		for j in range(0, N + 1):
			u[j] = -0.5 * (-1 + j * 2 / N)
			x[j] = -1 + j * 2 / N

		burgers = Burgers(visco[i], time, N)
		uall = burgers.macCormack(u)

		# Save 2-D data to files
		output = pd.DataFrame({'x': x, 't='+str(time_visco['time'][0]): uall[int(time_visco['time'][0] * 10000)], \
			't='+str(time_visco['time'][1]): uall[int(time_visco['time'][1] * 10000)], \
			't='+str(time_visco['time'][2]): uall[int(time_visco['time'][2] * 10000)]})
		output.to_csv('2D' + seper + 'mu=' + str(visco[i]) + '.dat', sep='\t', index=False)

此外，由于在每次时间步长的计算中记录了计算结果，所以可以将计算结果输出为$x, t, u$格式的三维形式，由于数据量较大，尝试过后发现一个粘性系数的三维数据的输出有160 MB，故直接采用python的matplotlib进行三维图像绘制，仅供参考分析使用。

参考运行时间：

2D数据输出：6.0 s
2D数据和3D图像输出：20.0 s

2.3. 运行结果展示

Fig 5. $\mu=0.0$ 的数值模拟结果

Fig 6. $\mu=0.01$ 的数值模拟结果

Fig 7. $\mu=0.05$ 的数值模拟结果

2.4. 结果分析

随着时间的推进，速度沿空间的发展越来越陡；
无粘情况下 ($\mu=0.0$)速度到最后会发生骤变，梯度有个猛然的变化，并且在速度骤变区域附近会产生振荡；
$\mu=0.01$时, 速度梯度变化没有无粘情况下那么剧烈，速度变化曲线更为平缓，激波的非物理振荡被消除；
$\mu=0.05$时，速度梯度变化更小，且在时间推进过程中，速度曲线间的差异越来越小；
总结：随着粘性的增大，粘性带来的耗散作用消除了速度的梯度的剧烈变化，并消除了解的振荡；而粘性趋近于零时，解越容易出现激波形式的不连续性。

3. 准一维Laval管流动数值求解

3.1. 输入文档说明

输入文档格式如下：

Mach Number,xleft,xright,Pe,Cx,CFL,Errlimit,Iterlimit,xStep
1.16,0.1,1.0,0.8785,0.05,0.5,0.000001,20000,200

Mach Number表示入口马赫数；
xleft表示入口坐标；
xright表示出口坐标；
Pe表示出口压力；
Cx表示人工粘性系数；
CFL表示Courant数；
Errlimit表示最大误差限；
Iterlimit表示最大迭代次数；
xStep表示空间离散的空间步。

3.2. 运行文档laval.py说明

采用守恒型控制方程进行求解，所以需要使用耦合和解耦

该文档定义一个类：

class Laval():
  # 输入人工粘性，Courant数，空间步长，离散后的空间坐标点，离散后空间截面面积
  def __init__(self, Cx, CFL, dx, x, Ax)
  # 耦合u，p，rho为U，F，H
  def Couples(self, u, p, rho)
  # 将U解耦为u，p，rho
  def deCouples(self, U)
  # 使用MacCormack格式求解
  def macCormack(self, u, p, rho, errlimit, iterlimit)

main函数实现：

if __name__ == '__main__':
  sep = os.sep
  input = pd.read_csv('input.txt', sep = ',', usecols=['Mach Number', 'xleft', 'xright', \
  'Pe', 'Cx', 'CFL', 'Errlimit', 'Iterlimit', 'xStep']) # Read csv file as input. Attention! inputs are lists.
  Ma, xl, xr, pe, Cx, CFL, errlimit, iterlimit, N = (input['Mach Number'][0], input['xleft'][0], \
    input['xright'][0], input['Pe'][0], input['Cx'][0], input['CFL'][0], input['Errlimit'][0], \
    input['Iterlimit'][0], input['xStep'][0])

  u, p, rho, x, Ax, dx = initCondition(xl, xr, N, pe, Ma)
  laval = Laval(Cx, CFL, dx, x, Ax)
  u, p, rho = laval.macCormack(u, p, rho, errlimit, iterlimit)
  Ma = u / (np.sqrt(1.4 * p / rho))
  flux = rho * u * Ax

  output = pd.DataFrame({'X/L':x, 'Mach': Ma, 'Velocity':u, 'Pressure':p, \
    'Density':rho, 'Flux':flux})
  output.to_csv('Result' + sep + 'Cx=' + str(Cx) + '.dat', sep='\t', index=False)

参考运行时间：5.0 s

3.3. 运行结果展示

Fig 8. $C_x=0.05$ 的数值模拟结果

Fig 9. $C_x=0.15$ 的数值模拟结果

Fig 10. $C_x=0.20$ 的数值模拟结果

## 3.4. 结果分析 ### 3.4.1. 人工粘性系数的影响从上面的图中可以看到：

人工粘性系数为0.05的时候，激波附近出现高频非物理振荡，沿管内流量不变；
人工粘性系数为0.15时，激波附近的高频非物理振荡得到抑制，沿管道内流量不变；
人工粘性系数为0.20时，激波附近的高频非物理振荡几乎得到抑制，沿管道内流量不变；
总结：随着人工粘性系数的增大，激波间断点处的梯度有所减小，斜率略有改变和倾斜，解的精度有所降低。

3.4.2. 与解析解比较

Fig 11. 数值模拟结果和解析解的比较

可以看出，随着人工粘性的增大，数值解有所发散，和解析解的差别有所增大，表明解的精度确实有所降低，人工粘性确实有数值耗散的作用。

孙北大 / NWPU-CFD

0. Introduction

1. O网格生成与数值求解

1.1. 输入文档说明

1.2. 运行文档ogrid.py说明

1.3. 运行结果展示

1.3.1. o网格生成图片

1.3.2. 流场显示

1.3.3. 翼型表面压力系数曲线

2. 一维Burgers方程数值求解

2.1. 输入文档说明

2.2. 运行文档burgers.py说明

2.3. 运行结果展示

2.4. 结果分析

3. 准一维Laval管流动数值求解

3.1. 输入文档说明

3.2. 运行文档laval.py说明

3.3. 运行结果展示

3.4.2. 与解析解比较

简介

发行版

贡献者

近期动态

孙北大 / NWPU-CFD .gitee-modal { width: 500px !important; }

0. Introduction

1. O网格生成与数值求解

1.1. 输入文档说明

1.2. 运行文档ogrid.py说明

1.3. 运行结果展示

1.3.1. o网格生成图片

1.3.2. 流场显示

1.3.3. 翼型表面压力系数曲线

2. 一维Burgers方程数值求解

2.1. 输入文档说明

2.2. 运行文档burgers.py说明

2.3. 运行结果展示

2.4. 结果分析

3. 准一维Laval管流动数值求解

3.1. 输入文档说明

3.2. 运行文档laval.py说明

3.3. 运行结果展示

3.4.2. 与解析解比较

简介

发行版

贡献者

近期动态

搜索帮助

孙北大 / NWPU-CFD