【CFD小工坊】尝试完成一个简单的溃坝流算例（1）
- 前言
- 算例简介
- 网格生成
- 数据的读入与输出
- 模型参数的读入
- 网格数据及结果数据的输出
前言
我们从一个简单的算例开始，从实际建模过程中学习和做代码。我选择的算例是一个矩形区域内的溃坝流，其效果如下图所示。图片来自于Caleffi等人1的数值模拟研究。
如图所示，模型的计算区域为一个正方形，中间处有一个有缺口的大坝。在初始时刻，大坝的一侧水位很高，另一侧水位很低（甚至可能没有水）；之后水流从坝口泄出，形成了溃坝流（dam break）。这个算例边界简单，常被用于测试模型对急变流、干湿边界处理等方面的模拟性能，是一个典型的测试算例。
我们的第一任务是实现模型的输入和输出。让模型能够读入我们的网格数据，并能按照一定格式输出模型的结果数据。
算例简介
本算例的计算区域是一个边长为200m的正方形区域，区域内有厚度为10m的大坝（位于
95
m
<
x
<
105
m
95m<x<105m
95m<x<105m处），如下图所示。计算域为平底，底高程为0m，其四周及大坝所对应的壁面处均采用壁面边界条件（closed boundary）。
在初始时刻，
x
<
100
m
x<100m
x<100m区域的水位为10.0m，其余模型参数见下表。
| 工况编号 | 上游水深 | 下游水深 | 曼宁系数 n n n | 模拟时长 |
|---|---|---|---|---|
| A | 10.0m | 0m（或Hwet） | 0.03 | 6.0s |
| B | 10.0m | 5.0m | 0.0 | 7.2s |
网格生成
计算网格可通过Triangle库（Triangle库的使用方式可见《学习使用Triangle库》）、SMS或DHI-MIKE等软件辅助生成。生成后，需要通过处理得到两类数据：
- 网格节点位置；包括节点总数、编号及坐标值。（用于制作points.txt文件）
- 网格几何关系；包括网格编号、构成各个三角形网格的节点序号。（用于制作cells.txt文件）
上述数据文件的格式详见《【CFD小工坊】模型网格（三角形网格）》。（该博文已于2024年3月更新）
我制作的网格如下图所示。网格节点总数3412，三角形网格总数6537。
之后编辑得到points.txt文件和cells.txt文件，它们的格式如下：
- points.txt（第一行是节点总数，之后依次是1~3412号节点的x、y坐标，一行数据对应一个节点）
3412
0 0
0 200
60 95
60 170
95 0
95 95
95 170
95 200
...
- cells.txt（第一行是三角形网格总数，之后依次是每个网格所对应三个节点的节点号；一行数据对应一个三角形网格）
6537
218 991 992
866 867 728
646 2171 2172
1810 1744 1809
143 2379 140
2306 2303 2304
3113 421 2481
4 1123 2775
2254 302 3334
...
数据的读入与输出
首先，我们确定一个主程序代码的框架：
int main() {
int i, j;
// Read the meshgrid
ReadGrid();
// Construct the grid
ConstructGrid();
// Read the parameters and other data
ReadInput();
// Initialization
//Initial();
// output at the initial time
output_plt(0);
// Main Loop
// Main Loop (not finished...)
//while (time <= time_end) {
//
//}
system("pause"); // just for test...
}
首先，调用函数ReadGrid()和ConstructGrid()来读入和构建网格。这两个程序的详细代码可见《【CFD小工坊】模型网格（三角形网格）》。关于网格的读入，不在此赘述了。
随后是通过函数ReadInput()读入一些模型参数，例如糙率、模拟时间等。
之后边是初始化程序并启动主循环while (time <= time_end) ，这一部分我们之后再实现。
此外，在运行前，我们试着通过函数output_plt( )输出一下读入的网格和初始场数据。当然，这个函数在后面输出模型结果的时候都会用到。
模型参数的读入
在模型中，参数文件是input.txt，即程序将读取我们写在input.txt中的数字，作为各个参数。
/*
* Function: ReadInput
* Usage: Read the data from the input files and other data files
* Called by: main
*/
void ReadInput() {
int i, j;
int bot_type; // bot_type= 0: zb=constant; bot_type=1: read zb at grid points; bot_type=2: read zb at cell centers.
int ini_eta; // ini_eta = 0: eta=constant (eta00) initially; ini_eta = 1: read eta from a file (eta0.txt).
int ini_vel; // ini_vel = 0: u,v=zero initially; ini_vel = 1: read u and v from a file (vel0.txt).
double bot0, eta00;
// 1. time parameters
time_str = get_double_val("TimeStr");
time_end = get_double_val("TimeEnd");
CFL = get_double_val("CFL");
dt_min = get_double_val("DtMin");
dt_max = get_double_val("DtMax");
// 2. bottom elevation
zbp = (double *)malloc(sizeof(double) * Np);
zbc = (double *)malloc(sizeof(double) * Nc);
bot_type = get_int_val("BotType");
if (bot_type == 0) {
// zb = constant
bot0 = get_double_val("Bot0");
for (i = 0; i < Np; i++) {
zbp[i] = bot0;
}
for (i = 0; i < Nc; i++) {
zbc[i] = bot0;
}
}
else if (bot_type == 1) {
// read zbp from zb.txt
// ...
}
}
else if (bot_type == 2) {
// read zbc from zb.txt
// ...
}
else {
fprintf(stderr, "Please verify the value of BotType in input.txt. \n");
exit(EXIT_FAILURE);
}
// 3. initial fields
eta = (double *)malloc(sizeof(double) * Nc);
u = (double *)malloc(sizeof(double) * Nc);
v = (double *)malloc(sizeof(double) * Nc);
char* row2 = (char*)malloc(sizeof(char) * N_str);
// 3.1 initial surface elevation
ini_eta = get_int_val("IniEta");
if (ini_eta == 0) {
eta00 = get_double_val("Eta0");
for (i = 0; i < Nc; i++) {
eta[i] = eta00;
}
}
else if (ini_eta == 1) {
// read water elevation from eta0.txt
// ...
}
else {
fprintf(stderr, "Please verify the value of Initial_Eta in input.txt. \n");
exit(EXIT_FAILURE);
}
// 3.2 initial velocities
ini_vel = get_int_val("IniVel");
if (ini_vel == 0) {
for (i = 0; i < Nc; i++) {
u[i] = 0.0;
v[i] = 0.0;
}
}
else if (ini_vel == 1) {
// read velocities from vel0.txt
// ...
}
else {
fprintf(stderr, "Please verify the value of Initial_Vel in input.txt. \n");
exit(EXIT_FAILURE);
}
free(row2);
// 4. open boundary
// ...
// 5. physical parameters
vis0 = get_double_val("Viscosity");
n_fric = get_double_val("Manning");
// 6. numerical parameters
DepDry = get_double_val("DepDry");
DepWet = get_double_val("DepWet");
theta_fric = get_double_val("ThetaFric");
// 7. output setups
out_str = get_double_val("OutStr");
out_int = get_double_val("OutInt");
fprintf(stderr, "The input.txt has been read... \n");
}
针对所需读入参数的类型，我将该函数段的代码分了七个部分。每个部分中用的读入方式有两种，第一种是通过get_val系列函数（这个系列函数会在之后的部分中讲解），它的功能就是从input.txt中读取所需名称的参数；第二种是读入一个初始场，如上述代码种的zb.txt，eta0.txt、vel0.txt。此外，我们还在这里设置了若干标识符，例如bot_type、ini_eta等，前者的意思是我们初始地形的设置方式，后者则是设置初始水位场的方式。
目前，我们暂不启用这些功能，它们有待我们的后续开发。目前，我们先定义我们的初始水位、流速都为0；其余参数暂时不读入。
网格数据及结果数据的输出
我们暂时采用tecplot作为后处理工具，故我们需要输出tecplot格式的结果数据文件和网格文件。代码如下：
/*
* Function: output_plt(count_p)
* Usage: output the results (water depth and velocities) as a tecplot file
* Called by: main
*/
void output_plt(int count_p) {
int i;
char filename[100];
// 1. time
FILE* fp1 = fopen(".\\result\\time.txt", "at");
fprintf(fp1, "%d\t%lf \n", count_p, time);
fclose(fp1);
// 2. create a tecplot file (.plt)
sprintf(filename, ".\\result\\%d.plt", count_p);
FILE* fp2 = fopen(filename, "at");
fprintf(fp2, "TITLE=\"Test\"\n");
fprintf(fp2, "VARIABLES=\"X\",\"Y\",\"Zbot\",\"elev\",\"u\",\"v\"\n");
fprintf(fp2, "ZONE T=\"ZONE1\",N=%d,E=%d,DATAPACKING=BLOCK,ZONETYPE=FETRIANGLE\n", Np, Nc);
fprintf(fp2, "VARLOCATION=([3-6]=CELLCENTERED)\n");
for (i = 0; i < Np; i++) {
fprintf(fp2, "%f\n", xp[i]);
}
for (i = 0; i < Np; i++) {
fprintf(fp2, "%f\n", yp[i]);
}
for (i = 0; i < Nc; i++) {
fprintf(fp2, "%f\n", zbc[i]);
}
for (i = 0; i < Nc; i++) {
fprintf(fp2, "%f\n", eta[i]);
}
for (i = 0; i < Nc; i++) {
fprintf(fp2, "%f\n", u[i]);
}
for (i = 0; i < Nc; i++) {
fprintf(fp2, "%f\n", v[i]);
}
for (i = 0; i < Nc; i++) {
fprintf(fp2, "%d\t%d\t%d\n", node[0][i] + 1, node[1][i] + 1, node[2][i] + 1);
}
fclose(fp2);
fprintf(stderr, "Output %d (Tecplot file) \n", count_p);
}
函数output_plt(count_p)自带参数count_p表示输出结果的次数，若是第一次输出结果，则count_p = 1；若是输出初始场，count_p = 0。之后我们需要在程序.exe同一目录下新建一个result文件夹，用于存放所有的结果文件。
第一类结果文件是时间序列，它记录了每一次输出结果是，对应的模拟时刻。结果存储在result文件夹下的time.txt里。第二类数据是tecplot数据.plt，它需要按照一定的格式，依次输入节点坐标、网格中心的水位、流速等等，最后还要附上网格的几何数据。
通过运行整个模型，我们得到的输出结果如下。我们在初始化时，预设了0的水位场和流速场，因此云图显示全局数据为0。
至此，我们初步实现了网格数据、模型参数读取和结果数据输出的功能！
Valerio Caleffi , Alessandro Valiani & Andrea Zanni (2003) Finite volume method for simulating extreme flood events in natural channels, Journal of Hydraulic Research, 41:2, 167-177, DOI: 10.1080/00221680309499959 ↩︎
