【CFD小工坊】尝试完成一个简单的溃坝流算例(1)

2024-03-04 18:12

本文主要是介绍【CFD小工坊】尝试完成一个简单的溃坝流算例(1),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

【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模拟时长
A10.0m0m(或Hwet0.036.0s
B10.0m5.0m0.07.2s

网格生成

计算网格可通过Triangle库(Triangle库的使用方式可见《学习使用Triangle库》)、SMS或DHI-MIKE等软件辅助生成。生成后,需要通过处理得到两类数据:

  1. 网格节点位置;包括节点总数、编号及坐标值。(用于制作points.txt文件)
  2. 网格几何关系;包括网格编号、构成各个三角形网格的节点序号。(用于制作cells.txt文件)

上述数据文件的格式详见《【CFD小工坊】模型网格(三角形网格)》。(该博文已于2024年3月更新)

我制作的网格如下图所示。网格节点总数3412,三角形网格总数6537。
在这里插入图片描述
之后编辑得到points.txt文件和cells.txt文件,它们的格式如下:

  1. points.txt(第一行是节点总数,之后依次是1~3412号节点的x、y坐标,一行数据对应一个节点)
3412
0	0
0	200
60	95
60	170
95	0
95	95
95	170
95	200
...
  1. 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 meshgridReadGrid();// Construct the gridConstructGrid();// Read the parameters and other dataReadInput();// Initialization//Initial();// output at the initial timeoutput_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 parameterstime_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 elevationzbp = (double *)malloc(sizeof(double) * Np);zbc = (double *)malloc(sizeof(double) * Nc);bot_type = get_int_val("BotType");if (bot_type == 0) {// zb = constantbot0 = 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 fieldseta = (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 elevationini_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 velocitiesini_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 parametersvis0 = get_double_val("Viscosity");n_fric = get_double_val("Manning");// 6. numerical parametersDepDry = get_double_val("DepDry");DepWet = get_double_val("DepWet");theta_fric = get_double_val("ThetaFric");// 7. output setupsout_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. timeFILE* 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。
在这里插入图片描述
至此,我们初步实现了网格数据、模型参数读取和结果数据输出的功能!


  1. 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 ↩︎

这篇关于【CFD小工坊】尝试完成一个简单的溃坝流算例(1)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/774021

相关文章

csu 1446 Problem J Modified LCS (扩展欧几里得算法的简单应用)

这是一道扩展欧几里得算法的简单应用题,这题是在湖南多校训练赛中队友ac的一道题,在比赛之后请教了队友,然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意:把题意转变下就变成了:求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点:求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个

hdu2289(简单二分)

虽说是简单二分,但是我还是wa死了  题意:已知圆台的体积,求高度 首先要知道圆台体积怎么求:设上下底的半径分别为r1,r2,高为h,V = PI*(r1*r1+r1*r2+r2*r2)*h/3 然后以h进行二分 代码如下: #include<iostream>#include<algorithm>#include<cstring>#include<stack>#includ

usaco 1.3 Prime Cryptarithm(简单哈希表暴搜剪枝)

思路: 1. 用一个 hash[ ] 数组存放输入的数字,令 hash[ tmp ]=1 。 2. 一个自定义函数 check( ) ,检查各位是否为输入的数字。 3. 暴搜。第一行数从 100到999,第二行数从 10到99。 4. 剪枝。 代码: /*ID: who jayLANG: C++TASK: crypt1*/#include<stdio.h>bool h

uva 10387 Billiard(简单几何)

题意是一个球从矩形的中点出发,告诉你小球与矩形两条边的碰撞次数与小球回到原点的时间,求小球出发时的角度和小球的速度。 简单的几何问题,小球每与竖边碰撞一次,向右扩展一个相同的矩形;每与横边碰撞一次,向上扩展一个相同的矩形。 可以发现,扩展矩形的路径和在当前矩形中的每一段路径相同,当小球回到出发点时,一条直线的路径刚好经过最后一个扩展矩形的中心点。 最后扩展的路径和横边竖边恰好组成一个直

poj 1113 凸包+简单几何计算

题意: 给N个平面上的点,现在要在离点外L米处建城墙,使得城墙把所有点都包含进去且城墙的长度最短。 解析: 韬哥出的某次训练赛上A出的第一道计算几何,算是大水题吧。 用convexhull算法把凸包求出来,然后加加减减就A了。 计算见下图: 好久没玩画图了啊好开心。 代码: #include <iostream>#include <cstdio>#inclu

uva 10130 简单背包

题意: 背包和 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#include <stack>#include <vector>#include <queue>#include <map>

Java 后端接口入参 - 联合前端VUE 使用AES完成入参出参加密解密

加密效果: 解密后的数据就是正常数据: 后端:使用的是spring-cloud框架,在gateway模块进行操作 <dependency><groupId>com.google.guava</groupId><artifactId>guava</artifactId><version>30.0-jre</version></dependency> 编写一个AES加密

JAVA用最简单的方法来构建一个高可用的服务端,提升系统可用性

一、什么是提升系统的高可用性 JAVA服务端,顾名思义就是23体验网为用户提供服务的。停工时间,就是不能向用户提供服务的时间。高可用,就是系统具有高度可用性,尽量减少停工时间。如何用最简单的方法来搭建一个高效率可用的服务端JAVA呢? 停工的原因一般有: 服务器故障。例如服务器宕机,服务器网络出现问题,机房或者机架出现问题等;访问量急剧上升,导致服务器压力过大导致访问量急剧上升的原因;时间和

简单的角色响应鼠标而移动

actor类 //处理移动距离,核心是找到角色坐标在世界坐标的向量的投影(x,y,z),然后在世界坐标中合成,此CC是在地面行走,所以Y轴投影始终置为0; using UnityEngine; using System.Collections; public class actor : MonoBehaviour { public float speed=0.1f; CharacterCo

docker-compose安装和简单使用

本文介绍docker-compose的安装和使用 新版docker已经默认安装了docker-compose 可以使用docker-compose -v 查看docker-compose版本 如果没有的话可以使用以下命令直接安装 sudo curl -L https://github.com/docker/compose/releases/download/1.16.1/docker-c