本发明涉及柔性植被消浪的数值计算方法，尤其涉及一种柔性植被消浪的参数化计算方法。
背景技术：
1、基于计算流体力学模型（cfd）的柔性植被消浪的数值计算方法主要可分为直接计算法和参数化计算方法。近年来，部分学者通过将计算流体力学方法和有限元方法结合实现了柔性植被与流体相互作用的流固耦合计算。上述直接计算方法由于需要刻画植被细节特征或者需要计算流固耦合过程，导致计算量过大，难以实现高密度植被消浪计算，无法满足工程需求。而参数化计算方法则是将植被区域概化为一个阻力区域，通过在控制方程中将纳入植被的阻力源项来实现植被消浪的计算。在该方法中无需刻画每一根植物，因此具有较高的计算效率。参数化计算方法的控制方程主要为ns方程（纳维-斯托克斯方程、navier-stokes equations），同时结合vof方程（volume of fluid方程）计算自由液面的波动。
2、目前的参数化计算方法在cfd计算中一般是先划定好植被的区域（计算前先指定哪些网格属于植被区域），在计算过程中对处于植被区域内的网格求解中施加阻力项和惯性力项。在计算过程中植被区域的范围是不可变的（如图 1，参数化的植被区域（红色区域）以及波面线（黑线）），因此在cfd中的植被消浪参数化模型主要用于计算刚性植被的消浪。此外在整个植被区域内拖曳力系数同样保持不变，该值主要通过一些经验性公式或者通过和实验数据率定来确定。
3、一些学者也通过引入等效高度的概念将参数化计算方法应用于柔性植被的计算中。具体而言，考虑波浪作用下柔性植被的弯曲而导致的植被高度降低的现象，将弯曲后的柔性植被视为一定高度下的刚性植被（该高度通常小于静水条件下的柔性植被自身高度），使得柔性植被和在该高度下的刚性植被具有相同的消浪能力，该高度成为等效高度。
4、但是，上述参数化计算方法在计算柔性植被时仍然存在如下的问题：1、目前的等效高度方法中，利用提前计算好的等效高度来确定消浪计算过程中的植被区域，并且所确定的柔性植被区域在计算过程中仍然是保持不变的，这个一点与实际情况并不相符，实际中柔性植被会随着波浪运动不断起伏。此外，计算等效高度的公式中，需要用到波浪最大水平流速以及波浪最大水质点位移，现有的方法中多是选取植被区域的前缘的波高、周期来确定波浪的流速和水质点位移，计算时将通过该数据计算出的等效高度应用与整个植被区域。但是在实际情况中，如果植被区域较长，波浪在植被区域内是不断衰减的，因此植被的等效高度同样应该是不断变化的，现有的模型无法准确的反映柔性植被瞬时等效高度的变化。2、植被区域拖曳力系数的确定也存在和等效高度类似的问题，即的值在植被区域内是保持不变的，但是由阻力源项的计算公式可知的值与植被瞬时受力以及瞬时局部流速相关（如图2所示，来自文献hu z, suzuki t, zitman t, et al. laboratorystudy on wave dissipation by vegetation in combined current–wave flow[j].coastal engineering, 2014, 88: 131-142.），恒定的难以准确的反映波浪在植被区域内部所受到的阻力变化特征，因此难以准确的计算植被区域内部流场变化特征（尤其是对于植被区域内存在泥沙输移的情况）；此外，多数模型中为了获取更好的消浪波高变化计算结果，的值都是通过和实验数据对比来率定出来的，但是对于缺乏实验数据的计算中（如在对原型尺度的计算），的确定较为困难。
5、鉴于此，特提出此发明。
技术实现思路
1、本发明的目的是为了解决现有技术中存在的缺点，而提出的一种柔性植被消浪的参数化计算方法，在计算的过程中，可以根据瞬时垂向平均流速和波面位移计算柔性植被的瞬时等效高度和拖曳力系数，经验证本发明的方法能够较好的预测柔性植被区域的消浪过程。
2、为了实现上述目的，本发明采用了如下技术方案：
3、一种柔性植被消浪的参数化计算方法，其特征在于，包括如下步骤：
4、s1：利用openfoam建立计算区域网格；
5、s2：根据静水条件下植被总高度，通过setfields和toposet工具确定植被区域内网格的集合cellzone；
6、s3：通过toposet和boxtoface工具生成网格面集合facezone；
7、s4：通过前处理程序cellfaceproject遍历cellzone投影后的网格中心点，获取体标量矩阵cellfaceaddressing，并将cellfaceaddressing写入初始时刻文件夹；
8、s5：读取字典文件中的参数，包括：植被茎宽度、植被茎高度、单位面积内植被茎的密度、植被叶片宽度、植被叶片厚度、植被叶片长度、每个茎上叶片的根数、植被茎和叶片的杨氏模量、静水深度、采样线长度、采样插值方法、拖曳力计算公式中的系数 a、b和c；
9、s6：根据cellfaceaddressing中各元素值构建插值采样线；
10、s7：计算各插值采样线上的瞬时垂向平均水流流速和瞬时波面位移，并计算各插值采样线处的植被的等效高度；
11、s8：根据各插值采样线处的植被的等效高度获取体标量矩阵vegetationcell；
12、s9：遍历cellzone，通过vegetationcell中值为1的网格计算拖曳力系数，并根据计算得到的拖曳力系数求解ns方程和vof方程。
13、进一步，所述s4包括如下步骤：
14、s41：将facezone中每一个面的顶点坐标投影到z=0的水平面，得到若干投影后的面；
15、s42：将cellzone中每一个网格中心点坐标投影到z=0的水平面，得到若干投影后的网格中心点；
16、s43：构建体标量矩阵cellfaceaddressing，其中各元素与s42中cellzone投影后的网格中心点一一对应，其中各元素初始值为-1；
17、s44：遍历s42中cellzone投影后的网格中心点，通过射线法判断其是否位于s41中facezone投影后的任何一个面，如果是则将cellfaceaddressing中与该投影后网格中心点对应的元素赋值为与其对应的facezone投影后的面的编号值；
18、s45：遍历s42中cellzone投影后的各网格中心点后，通过零梯度边界条件更新cellfaceaddressing边界上的值，获得最终的cellfaceaddressing，并将其写入初始时刻文件夹。
19、进一步，所述s6包括如下步骤：
20、s61：通过c++标准库中的std::map函数建立cellfaceaddressing中元素与计算域底边界网格面心坐标之间的映射；
21、s62：计算具有相同的映射值的底边界网格面心坐标的坐标平均值，并将该坐标平均值作为插值采样线的底坐标；
22、s63：沿垂向向上平移预设范围内的距离作为插值采样线的顶坐标；
23、s64：通过midpointandfaceset类构建各网格簇的插值采样线。
24、进一步，在所述s61前还包括如下步骤：
25、s60：读取cellfaceaddressing中各元素的值以及计算区域网格底边界上所有网格面心坐标。
26、进一步，在所述s61前还包括如下步骤：
27、s601：通过pstream::gatherlist函数获取每一个处理器上的cellfaceaddressing中各元素的值以及计算区域网格底边界上所有网格面心坐标；
28、在所述s61中，在主处理器上通过c++标准库中的std::map函数建立cellfaceaddressing中元素与网格面心坐标之间的映射；
29、在所述s62中，将主处理器中的插值采样线的底坐标的集合发送给各其余处理器；
30、在所述s64后还包括如下步骤：
31、s65：通过globalindex类中的gatheronly函数，使主处理器上获取存储于各其余处理器中的插值采样线的采样点坐标集合，并利用coordset中的gathersort函数在主处理器上建立插值采样线上各采样点顶坐标和底坐标的排序映射。
32、进一步，所述s7包括如下步骤：
33、s71：在每一个时间步长，通过设定的采样插值法在各插值采样线上对流速u和液体体积分数α进行插值；
34、s72：通过梯形积分法对u和α进行积分，分别获取各插值采样线上的瞬时垂向平均流速u(t)和瞬时水深，将并瞬时水深与静水深的差值的绝对值作为瞬时波面位移a(t)；
35、s72：根据各插值采样线上的瞬时垂向平均流速和瞬时波面位移计算各插值采样线处的植被的等效高度。
36、进一步，所述s72通过如下步骤计算植被的等效高度：
37、s721：计算植被茎柯西数和植被叶片柯西数，计算公式如下：
38、，
39、，
40、，
41、，
42、公式中，为水的密度，为植被茎的宽度，为植被叶片的宽度，为植被茎的高度，为植被叶片的高度，为植被茎的杨氏模量，为植被叶片的杨氏模量，为植被叶片的厚度；
43、s722：计算植被茎部的等效高度和植被叶片的等效高度，计算公式如下：
44、，
45、，
46、，
47、，
48、公式中，k为常数，一般可取0.94；
49、s373：计算植被的等效高度，计算公式如下：
50、，
51、公式中，为单个茎上的叶片数。
52、进一步，所述s71包括如下步骤：
53、s711：在每一个时间步长，通过设定的采样插值法在各插值采样线上对流速u和液体体积分数α进行插值；
54、s712：在主处理器上通过globalindex类中的gatherinplace函数获取各其余处理器上的插值结果和对应的采样线上插值点坐标；
55、s713：在主处理器上通过coordset类中的gathersort函数中的排序映射对获取各其余处理器上的插值结果和对应的插值点坐标进行排序；
56、s72包括如下步骤：
57、s721：在主处理器上通过梯形积分法对u和α进行积分，分别获取各插值采样线上的瞬时垂向平均流速u(t)和瞬时水深，将并瞬时水深与静水深的差值的绝对值作为瞬时波面位移a(t)；
58、s722：将主处理器上计算得到的u(t)和a(t)发送至各其余处理器。
59、进一步，所述s8包括如下步骤：
60、s81：构建体标量矩阵vegetationcell，其中各元素与cellfaceaddressing中各元素一一对应，其中各元素初始值为0；
61、s82：判断cellfaceaddressing中各元素对应的网格簇内网格中心高度是否小于与其对应的插值采样线处的植被的等效高度，如果是则将vegetationcell对应元素置1。
62、进一步，在s9中通过如下公式计算拖曳力系数：
63、
64、，
65、公式中，a、b和c为常数，μ为水的动力粘度值。
66、本发明与现有技术相比，其有益效果为：
67、1、在计算的过程中，在cellzone中根据瞬时垂向平均流速和波面位移计算柔性植被的瞬时等效高度，根据瞬时等效高度可确定计算区域网格中起到消浪作用的柔性植被区域；此外，在cellzone中根据瞬时垂向平均流速计算拖曳力系数，可实现高效率、高精度的柔性植被计算；
68、2、经验证计算所得的有效波高和拖曳力系数符合实际中的植被运动情况，能够较为准确的计算植被区域内的沿程波高衰减变化趋势，能够较好的预测柔性植被区域的消浪过程。
1.一种柔性植被消浪的参数化计算方法，其特征在于，包括如下步骤：
2.根据权利要求1所述的一种柔性植被消浪的参数化计算方法，其特征在于，所述s4包括如下步骤：
3.根据权利要求2所述的一种柔性植被消浪的参数化计算方法，其特征在于，所述s6包括如下步骤：
4.根据权利要求3所述的一种柔性植被消浪的参数化计算方法，其特征在于，在所述s61前还包括如下步骤：
5.根据权利要求3所述的一种柔性植被消浪的参数化计算方法，其特征在于，在所述s61前还包括如下步骤：
6.根据权利要求3所述的一种柔性植被消浪的参数化计算方法，其特征在于，所述s7包括如下步骤：
7.根据权利要求6所述的一种柔性植被消浪的参数化计算方法，其特征在于，所述s72通过如下步骤计算植被的等效高度：
8.根据权利要求6所述的一种柔性植被消浪的参数化计算方法，其特征在于，所述s71包括如下步骤：
9.根据权利要求6所述的一种柔性植被消浪的参数化计算方法，其特征在于，所述s8包括如下步骤：
10.根据权利要求9所述的一种柔性植被消浪的参数化计算方法，其特征在于，在s9中通过如下公式计算拖曳力系数：
