溃漫堤洪水多维耦联数值模型及应用

苑希民，王亚东，田福昌



溃漫堤洪水多维耦联数值模型及应用

苑希民，王亚东，田福昌

(天津大学水利工程仿真与安全国家重点实验室，天津 300350)

针对山洪沟防洪标准偏低、遭遇暴雨洪水时溃漫堤风险较大且较难准确预测的问题，开展洪水多维耦联数值模型研究．基于圣维南方程与VOF法的标准-双方程原理，融合Abbott六点隐格式法、非结构化网格Roe格式的单元中心显式有限体积法与结构化网格PISO算法优势，系统构建溃漫堤洪水多维耦联数值仿真模型．以贺兰山汝箕沟为研究对象，模拟其遭遇百年一遇洪水时，沟道一维、溃口三维与淹没区域二维的多维水流动态耦合演进过程及影响情况．结果表明，溃漫堤洪水多维运动过程与流态模拟效果较好，所建多维洪水耦联数值模型计算精度较高，对多流态洪水动态耦合精确计算与淹没风险准确评估具有重要意义．

溃漫堤洪水；多维耦联；动态耦合；数值模型；汝箕沟

随着自然环境日趋恶劣，极端天气常有发生，致使超常洪水较以往更加频繁．堤防在防洪减灾方面起到了至关重要的作用，但当遭遇超标准洪水时，漫溃堤风险加大．水流数值模拟技术具有计算精度高、水力特征信息提取便利等明显优势，可为防汛指挥决策提供可靠的技术支撑．

近年来，各国学者对洪水演进数值模拟技术进行了广泛深入的研究．在一、二维水动力模型及其耦合研究方面，Blade等[1]基于数值通量耦合方法，分析模拟了天然河道水流运动情况；张靖雨等[2]采用一、二维水动力耦合模型对溃堤洪水与内涝洪水进行叠加分析，模型能够合理反映保护区洪涝风险分布特征，结果较为可靠合理；苑希民等[3]将全二维气相色谱理论与全二维水动力模型结合，为洪水演进模拟提供了新思路，并将该模型成功应用于黄河内蒙段河道与左右岸灌区洪水风险动态耦合模拟；施露等[4]运用Mike Flood建立了水文学模型、一维和二维水动力耦合模型，模拟分析了中小河流洪涝风险情况，研究结果具有重要价值．在二、三维水动力模型及其耦合研究方面，练继建[5]基于VOF法建立了溃坝水流运动三维数值模型，并对河道丁坝附近溃坝波传播的特性及其内部结构进行了模拟研究；王桂芬[6]提出了二、三维嵌套数学模型概念，并应用于天津新港抛泥区的潮流计算；Zounemat等[7]也对二、三维耦合水动力模型进行了研究，其耦合方式为水体表面采用二维浅水模型计算，而深水水体采用三维水动力模型计算．综上所述，国内外学者对溃堤或漫堤洪水一、二维耦合模拟研究较多，少数学者对二、三维水动力耦合模型进行了理论与应用研究，而对于溃漫堤洪水多维耦联计算模型的研究却少见报导．实际堤防溃决洪水分流是极其复杂的三维非恒定流运动过程，堤防及溃口边界对水流流态与分流过程影响很大，通常所建一、二维水动力耦合模型仅能计算河道与溃口断面平均水力要素，难以精确获取纵向水流变化规律与分流 特征．

因此，本文拟融合不同维度水动力数值模型计算优势，基于河道实测断面与DEM数字高程，采用Abbott六点隐格式法、Roe格式离散法与PISO算法，研究建立河道一维、平面区域二维与溃口分流三维的洪水多维耦联数值仿真模型，较为准确地计算沟道一维洪水水面线、淹没区二维洪水淹没特征与溃口分流三维洪水运动过程，并将其应用于贺兰山汝箕沟遭遇百年一遇超标准洪水时溃漫堤洪水演进动态变化的数值模拟．

1 模型原理

1.1 河道一维水动力模型

河道一维水动力学模型采用描述明渠非恒定流的圣维南方程，如式(1)和式(2)所示[8]．

 (2)

式中：为过水面积；为区间流量；为断面过水流量；为重力加速度；为时间；为沿程距离；,为动量修正系数；为河道断面水位；为流量模数．采用稳定性与精准度较高的Abbott六点隐式格式离散上述方程．

1.2 淹没区二维水动力模型

二维平面淹没区域，垂向洪水流速等水力要素的变化远小于水平方向，故采用二维浅水方程对其进行描述计算．具体如式(3)～式(5)所示[3]．

 (3)

 (4)

 (5)

式中：为水深；为水位；为源汇项流量；与分别为、方向上的垂向平均单宽流量；、分别为垂向平均流速在、方向的分量；为曼宁系数；为重力加速度．方程不考虑科氏力和风力影响．

在保证计算区内水量和动量守恒的基础上，采用Roe格式的单元中心显式有限体积法对控制方程进行离散求解，并通过干湿边界处理技术、线状地物与网格耦联方法，优化淹没区二维水动力模型．

1)干湿边界处理优化

为提高模型计算效率及其稳定性，采用干湿边界理论判断网格淹没状态，通过设置干水深及湿水深，确定洪水移动边界．当网格单元计算水深值小于干水深时，单元通量为0；当计算水深值大于干水深而小于湿水深时，网格单元只计算质量通量，无动量通量；而当其大于湿水深时，网格单元质量通量与动量通量均进行计算处理．具体关系如图1所示．

图1 网格水深与通量计算关系

2)线状地物与网格耦联方法

为在保护区内准确定位线状地物，将实测地理参考点线性连接，以其为内边界创建网格，可使堤防、道路等特殊边界与非结构化网格密切耦联．该方式可避免线状地物与网格界面相交，使得参与计算折线与实际线状地物重合，利于模型稳定运行并提高计算精度．图2(a)为常规线状地物耦联方式，图2(b)为改进后的线状地物耦联示意图．

图2 线状地物与网格耦联示意

1.3 溃口分流三维水动力模型

洪水经由河道通过溃口进入堤外保护区域时，水流方向发生急剧改变，堤防受到溃堤水流的强烈冲击，水流流态复杂，三维特征十分明显．一维水动力模型计算溃口分流洪水时，只能得到纵向上无梯度分布的河道及溃口断面平均水力要素，与实际情况差异较大．故本文采用具有较强界面捕捉能力的VOF方法的标准双方程模型对溃口洪水分流过程进行模拟计算[9]．基本控制方程如式(6)和式(7)所示．

 (6)

 (7)

流体自由表面处理方程如式(8)～式(10)所示．

输运方程为

 (8)

式中：为流体体积函数，＝1表示控制体内充满流体，＝0表示无流体；u为x方向上的流速．

紊流方程为

 (9)

 (10)

模型常数

1z＝1.44，2z＝1.92，k＝1.0，s＝1.3

式中：为分子动力黏性系数；t为紊流黏性系数；为紊动动能；为紊动能生成率；为紊动耗散率；k、e分别为、的紊流普朗特数；1z、2z为经验常数．在时间与空间上分别采用Crank-Nicholson格式和LUD格式离散求解控制方程．为使其能更好地满足动量方程和连续性方程，采用PISO算法对方程进行修正．具体见文献[9]．

1.4 多维洪水动态耦联模型

河流溃漫堤洪水耦合计算关系如图3所示．

图3 溃漫堤洪水耦合计算示意

1.4.1 一、二维水动力耦合模型

当河道水位超过堤顶高程时，发生洪水漫溢现象，可将漫溢洪水作为宽顶堰流处理[8]，建立一、二维水动力耦合模型，采用侧向连接方式实现河道与堤外保护区域的动态耦联，利用Villemonte堰流公式计算洪水漫溢流量(如式(11)所示)，模拟漫溢洪水在堤外平面区域的动态演进过程及淹没特征．

 (11)

1.4.2 一、三维水动力耦合模型

为减小河道下游水位顶托影响，采用河床具有正向坡的断面作为一、三维水动力模型的耦合衔接面，边界条件设定如下．

1)上游入流边界条件

以流量过程作为河道上游入流边界条件．为保证一、三维水动力耦合模型计算结果的准确性，考虑将衔接面河道水位作为模型动态耦合衔接条件，如式(12)、式(13)所示．由经验公式确定和，如式(14)和式(15)所示．

 (12)

 (13)

 (14)

 (15)

式中：1为衔接面上游流量；3为衔接面下游流量；1是一维河道水位；3为三维河道水位；3是三维河道入口水流流速；为一、三维耦合衔接面长度．

2)下游出流边界条件

以河道末端断面水位-流量关系作为下游出流边界条件．为确保模型计算稳定，连同溃口上、下游一定范围内的河道及淹没区作为三维水动力模型建模范围．在此可选择溃口下游地形平坦、水流流态较为平稳的淹没区域作为二、三维水动力耦合模型的过渡衔接面，如图4所示．

图4 二、三维耦合模型动态衔接示意

1.4.3 二、三维水动力耦合模型

在计算二、三维水动力耦合模型衔接通量时，三维模型沿水深方向具有多层计算节点，而二维模型仅计算水平方向速度分布，两者相差较大．计算二维模型对三维模型的反馈影响，需获得三维洪水纵向水力要素的变化规律，以此来补充其垂向分层边界条件．为确保流速和水位等变量信息在二、三维模型耦合衔接面上的交互传递[10](如图5所示)，在压强及流量相等的基础上，对三维洪水纵向速度分布进行平均处理，如式(16)和式(17)所示．

 (16)

 (17)

对二、三维耦合模型过渡衔接处回水反馈影响，按式(18)和式(19)对平均流速进行垂线化分布处理.

 (18)

 (19)

式中：为流速分布指数；为水深；为距离河床高程；、分别为、方向基于水深平均的流速；()、()分别为沿河床高程各层节点的、方向流速，水流流经二、三维耦合过渡衔接区域后，垂向无速度梯度．

图5 二、三维耦合模型衔接面流速示意

2 模型验证

为验证多维耦联模型的准确性，本文选用Soares等[11]于2002年做的物理试验对耦联模型进行验证．试验模型由水库段和渠道段组成，其中渠道宽0.495,m，渠底高于水库底部0.33,m，距溃口下游3.92,m处有一90°直角弯，渠道的曼宁糙率系数为0.011，末端为自由出流；水库宽2.44,m、长2.39,m，初始水位0.58,m，高出渠道底部0.25,m．该试验主要研究水库坝体瞬间全溃时水流在渠道中的演进过程.

在此分别建立多维耦联模型和一、二维耦合模型，如图6所示．首先，从上游至下游依次建立二维(A区)、三维(C区)和一维(D区)模型．当水流流经直角弯时，水体流态空间分布较为复杂，故C区域进行三维模拟．其次，因该耦联模型只有二、三维耦合和一、三维耦合，缺少对多维耦联模型中的一、二维耦合模型部分的验证，在此单独建立一、二维耦合模型，如图6(b)所示．模型处理：①采用非结构化三角形网格剖分二维模拟区，最大网格面积0.009,m2； ②二、三维过渡段区(B区)采用0.02,m×0.02,m×0.1,m，C区采用0.02,m×0.02,m×0.01,m结构化正交网格进行剖分；③一维模拟计算区，渠道断面距离0.02～0.04,m，在直角弯进行插值加密．

模型上边界条件为给定的初始水深0.58,m，下边界为出口端水位-流量关系．为使得水库水面水平，在渠道处设一控制闸门，当其上游水位等于0.58,m时，闸门突然提起；模型计算步长为0.01,s，二、三维模型最小步长为0.1,μs，模拟时长＝10,s．图7为水深验证结果，提取模型在＝5,s和＝7,s时的模拟水深值，将其与实测数据进行对比分析．由图可知，二维模拟计算值与实测值最为接近，因二维模拟区水流受直角弯影响较小，流态较为平稳，故该段模拟结果较为准确；而三维模拟水深总体较小，主要原因有：①涡黏性系数各向同性假定，即模型自身限制；②因模型尺寸总体较小，水体表面张力影响较大，而数学模型并未考虑表面张力作用，致使模型计算结果误差较大；对于一维渠道，由模拟结果可知，水深变化总体趋势基本一致，其中该段上游模拟值与实测值较为接近．结果表明，本文所建多维耦联数值计算模型模拟结果较为准确可靠．

A—夺维计算区；B—耦联过渡区；C—三维计算区；D—一维计算区；——耦合线

图7 水深验证结果

3 应用案例

3.1 研究区域概况

本文选择贺兰山汝箕沟出山口至长胜墩拦洪库段作为研究对象(如图8所示)，该河段地处贺兰山东麓中段地区石嘴山市，为宁夏两大暴雨中心之一．区域地势总体为西南高东北低，洪水自出山口流出后，经行洪沟道汇入长胜墩拦洪库．行洪沟道长约11.01,km，左右岸防洪堤防顶宽约5.00,m，沟道平均比降为3.29%，防洪标准为20年一遇．研究区域内主要线状地物有G110、姚汝公路和石银高速公路等，因姚汝公路桥上游转弯处河道弯度较大、水流对右侧凹岸淘刷作用较强，溃堤风险偏高，且堤防一旦溃决，将对溃口下游居民及大量农田造成严重影响，故将此处设定为溃口位置．

图8 区域概况

3.2 多维耦联模型的建立

根据河道洪水、溃口分流洪水、堤防漫溢洪水与淹没区洪水流动各自表现出的不同维度特征，建立溃漫堤多维洪水动态耦联数值计算模型．为减小河道一维与溃口三维水流动态交互影响，于溃口断面上游200,m处设置一、三维耦合面域，并通过二、三维耦合过渡衔接区，实现溃口三维水流与淹没区二维洪水的动态耦合，模拟河道洪水水面线、溃口三维洪水分流过程、溃漫堤洪水运动特征及其淹没风险．研究区网格剖分及模型衔接关系，如图9所示．

图9 二维区域网格剖分及模型衔接关系

汝箕沟沟道一维水动力模型计算河段长度11.01,km，设置控制断面192个．根据《石嘴山市汝箕沟大武口段治理工程初步设计报告》内容论述，设定汝箕沟沟道综合糙率为0.034．为提高模型计算稳定性及运算效率，设定汝箕沟一维水动力模型计算步长为1,s．选择2009年汝箕沟水文站汛期实测洪水过程作为典型洪水，并采用同倍比放大法求得其百年一遇设计洪水流量过程(如图10所示)，作为一、三维水动力耦合模型的上游入流边界条件．

采用非结构化三角形网格剖分堤外平面保护区域，区域总面积318.68,km2，网格面积控制在100～5,000,m2之间，剖分网格数量17.50×104个，计算节点数8.92×104个．计算区域多为农田耕地，房屋村庄面积较小且零散分布(见图11)，因缺乏实测糙率资料，在此依据以往类似区域模型的计算经验[12]，可采用不同地理条件的糙率值与面积权重之积，确定计算区域综合糙率值为0.06．为提高模型计算精度与运算效率，设置干水深值为0.005,m，湿水深值为0.1,m，在满足临界值CFL≤1的条件下，设置计算时间步长为0.01～10,s，模型可根据计算效率与稳定性大小，实现计算步长的自我动态调整．

图10 洪水流量过程

图11 土地利用图

针对姚汝公路桥上游转弯处溃口上、下游200,m河道及溃口外延50,m淹没区，采用矩形网格剖分技术构建研究区域三维地形，网格边长为0.60～1.50,m，剖分网格数量74.47×104个，整体建立溃口分流过程三维仿真计算模型．当溃口处河道断面水位等于设计水位1,106.10,m时，堤防瞬间溃决，并参考《洪水风险图编制技术细则(试行)附录》中的溃口宽度经验公式(如式(20))，溃口处河宽77.00,m，溃口宽度计算值为59.98,m，模型设定溃口宽度为60,m，溃口底高程为1,105.87,m．通过模拟溃口三维洪水分流及二、三维过渡衔接区域水流运动过程，实现多维洪水流态的动态耦合与过渡，准确计算溃漫堤分流洪水淹没风险分布特征．

 Bb＝1.9(lg,B)4.8＋20(20)

式中：b为溃口宽度，m；为河宽，m．

3.3 计算结果与影响分析

根据汝箕沟溃漫堤多维洪水耦联数值模型计算结果分析可知：当汝箕沟遭遇百年一遇洪水时，出山口下游850～2,100,m之间沟段发生洪水漫溢现象，漫溢洪水最高水位超过堤顶高程0.51,m．漫溢洪水流量过程如图12所示，漫溢分流呈现2个峰值，最大分流量为132,m3/s，漫溢总水量为17.70×104,m3，上游来水洪峰过后，河道流量迅速减小，堤防漫溢洪水流量急剧降低，且由于G110及下游沟道堤防阻水壅水作用明显，使得发生漫溢洪水倒灌回流现象，即漫溢流量出现短暂负值，随后水流平稳漫溢流量为零．

图12 漫溢洪水流量过程

汝箕沟溃口位于出山口下游3,150,m处，溃口分流流量过程如图13所示．堤防瞬间溃决，溃口分流流量急剧增加，起溃时分流流量为59,m3/s，分流洪峰流量为139,m3/s，溃堤分流总水量为68.35×104m3．汝箕沟洪水峰值时刻溃口处河道横断面与区域三维流速分布分别如图14和15所示．由图14分析可知：溃口附近流速较大，溃口外侧流速逐渐减小，水面流速低于底部流速；溃口出流水面收缩处流速最大，下游发生水跃现象，水面波动剧烈，流速减小，水深增加，水面逐渐趋于平稳．由图15分析可知：堤防对洪水破堤分流干扰作用明显，表现为侧向出流，且在溃口边缘部位产生漩涡，溃口处流速急剧变化且明显高于其余部位．

图13 溃口流量过程

图14 洪峰流量时溃口处河道横断面流速分布

图15 溃口三维分流流速分布

图16为研究区地形，汝箕沟溃漫堤洪水演进不同时刻淹没水深分布情况如图17所示，与图16对比分析可知：洪水淹没范围不断向地势低洼处扩展，其中以石银高速公路为例，因涵洞位置地形较低，路面较高，使得石银高速靠近溃口一侧发生积水现象，该处水深值较大，阻水作用明显．当洪水演进2,h，漫溢洪水淹没风险较大，主要影响区域为北庄子；当洪水演进6,h，主影响区域包括北庄子、崇岗村、金达煤业、崇岗镇民生服务中心和崇岗村砖厂；当洪水演进12,h，结合图11可知，溃口分流洪水淹没大量农田耕地，影响范围扩展至农三队．

图16 研究区地形图

按不同淹没水深等级(见表1)统计分析溃漫堤洪水淹没区域面积与受影响道路长度，如表1所示．洪水淹没总面积为6.70,km2，其中溃堤洪水淹没面积为5.84,km2，漫堤洪水淹没面积为0.86,km2，淹没水深小于0.5,m的风险区域面积为6.58,km2，占淹没总面积的98.21%,；受影响道路总长度为4.00,km，其中积水深小于0.5,m的道路长度为3.16,km，占受影响道路总长度的79.00%,．

图17 洪水演进不同时刻淹没水深分布

表1 不同水深等级下淹没面积与受影响道路统计表

Tab.1 Statistics table of submergedarea under different depth levels and roads affected

4 结 语

本文基于圣维南方程与VOF法的标准双方程，建立了溃漫堤洪水多维耦联计算模型．以贺兰山汝箕沟为研究对象，分析模拟了其遭遇百年一遇超标准洪水情况下，沟道一维洪水演进、溃口三维分流以及平面区域二维动态淹没的多维洪水耦合联动过程．结果表明，溃漫堤多维洪水演进过程模拟效果较好，流态较为真实，计算结果基本合理可靠．所建模型在水面追踪与多维洪水动态耦联模拟等方面具有明显优势，横纵向水力要素模拟与多空间尺度处理能力较高，具有很好的推广应用前景．

［1］ Blade，Gmez-Valentn M，Dolz J，et al. Integration of 1D and 2D finite volume schemes for computations of water flow in natural channels[J].，2012，42：17-29.

［2］ 张靖雨，徐 佳，袁先江，等. 外洪溃堤一、二维耦合模型与内涝模型叠加在防洪保护区内的应用探讨[J]. 水利水电技术，2017，48(5)：87-94.

Zhang Jingyu，Xu Jia，Yuan Xianjiang，et al. Discussion on application of superposition of 1-D and 2-D coupling model for dike failure due to outside flood and water-logging model to flood control protection area[J].，2017，48(5)：87-94(in Chinese).

［3］ 苑希民，田福昌，王丽娜. 漫溃堤联算全二维水动力模型及应用[J]. 水科学进展，2015，26(1)：83-90.

Yuan Ximin，Tian Fuchang，Wang Lina. Comprehensive two-dimensional associate hydrodynamic models for overflow and levee-breach flood and its application[J].，2015，26(1)：83-90(in Chinese).

［4］ 施 露，董增川，付晓花，等. Mike Flood在中小河流洪涝风险分析中的应用[J]. 河海大学学报：自然科学版，2017，45(4)：351-358.

Shi Lu，Dong Zengchuan，Fu Xiaohua，et al. Application of Mike flood and waterlogging risks of medium and small rivers[J].：，2017，45(4)：351-358(in Chinese).

［5］ 练继建. 溃坝水流在复杂河道中传播的三维数值模拟[J]. 水利学报，2007，38(10)：1151-1157.

Lian Jijian. Three-dimensional numerical simulation of dam-break flow in complicated river sections[J].，2007，38(10)：1151-1157(in Chinese).

［6］ 王桂芬. 二、三维潮流数学模型嵌套连接技术[J]. 交通与计算机，1988(3)：39-40.

Wang Guifen. Second-dimensional and three-dimensional trend mathematical model nested connection technology [J].，1988(3)：39-40(in Chinese).

［7］ Zounemat K M，Sabbagh-Yazdi S R. Coupling of two-and three-dimensional hydrodynamic numerical models for simulating wind-induced currents in deep basins[J].，2010，39：994-1011.

［8］ 苑希民，薛文宇，冯国娜，等. 溃堤洪水分析的一、二维水动力耦合模型及应用[J]. 水利水电科技进展，2016，36(4)：53-58.

Yuan Ximin，Xue Wenyu，Feng Guona，et al. A couple one-and two-dimensional hydrodynamic model for analysis of levee-breach flood and its application[J].2016，36(4)：53-58(in Chinese).

［9］ 李 超. 长距离调水工程高填方渠道溃堤洪水演进情景仿真[D]. 天津：天津大学建筑工程学院，2014.

Li Chao. Scenario Simulation on Dike Breach Flood of High Embankment Channel in Long Distance Water Diversion Project[D]. Tianjin：School of Civil Engineering，Tianjin University，2014(in Chinese).

［10］ 黄玉新，张宁川. 二、三维耦合水动力模型研究I：模型的建立[J]. 水道港口，2013，34(4)：304-310.

Huang Yuxin，Zhang Ningchuan. Research on a couple 2D-3D hydrodynamic model I：Model establishment [J].，2013，34(4)：304-310(in Chinese).

［11］ Soares Frazao S，Zech Y. Dam break in channels with 90° bend[J].，2002，128(11)：956-968.

［12］ 贺新娟，徐国宾，苑希民. 溃口宽度和倒灌、退水对防洪保护区洪水演进的影响[J]. 水资源与水工程学报，2016，27(2)：142-147.

He Xinjuan，Xu Guobin，Yuan Ximin. Effect of levee-breach width and water intrusion and recession on flood routine in flood protected zone[J].，2016，27(2)：142-147(in Chinese).

(责任编辑：王新英)

Application of Multidimensional Coupling Numerical Model for Collapse Flood

Yuan Ximin，Wang Yadong，Tian Fuchang

(State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300350，China)

In view of the problems of low flood control standard of ravine，high risk of catastrophic dike and difficulty in accurate pridiction，the numerical model of flood multidimensional coupling was studied．Based on the principle of Saint-Vernon equation and VOF method，the Abbott six-point implicit scheme，the element center explicit finite volume method of unstructured grid Roe scheme and the advantages of structured grid PISO algorithm were fully integrated．A multi-dimensional computational model was built systematically．Taking the Rujigou valley of Helan Mountain as an example，the evolution process and influence of multidimensional water flow in the 1D river，2D model of calculated area and 3D segment breach flood inundated areas were analyzed when it was struck by 100-year flood．Results show better simulation for the multidimensional motion process and flow of collapse flood，and higher calculation accuracy for the numerical model with multidimensional flood coupling，which is of great significance for the accurate calculation of dynamic coupling of multi-stream flooding and the risk assessment of flooding risk.

collapse flood；multidimensional coupling；dynamic coupling；numerical model；Rujigou

10.11784/tdxbz201707017

TV131.2

A

0493-2137(2018)07-0675-09

2017-07-08；

2017-08-19.

苑希民（1968—），男，博士，教授，yxm@tju.edu.cn.

田福昌，tianfuchang@tju.edu.cn.

国家重点研发计划资助项目(2017YFC0405601)；高等学校学科创新引智计划资助项目(B14012)；国家自然科学基金委创新团队资助项目(51621092).

the National Key Research and Development Program of China(No. 2017YFC0405601)，the Program of Introducing Talents of Discipline to Universities(No. B14012)and the Science Fund for Creative Research Groups of the National Natural Science Foundation of China(No. 51621092).