-
随着国家政策对新能源行业的倾斜和风能行业的蓬勃发展，海上风电由于风资源充足，地形影响小，距离电网消纳端近等特点，受到日益广泛的关注，加之技术发展降低了海上风电的度电成本，海上风电场建设正在进入快车道。海上风电快速发展的同时，仍然面临一些技术性问题有待解决，其中，环境载荷对风机基础结构的冲击作用，及其引发的风机结构安全性问题需要给予一定的重视。根据欧洲海上风电场的基础冲刷实测数据，包括Barrow、Scroby Sands和Otzumer Balje在内的多个风电场场址内，部分机位的基础冲刷深度达到或超过1.2倍桩径[1]。对于我国南方部分省份的沿海风电场，由于位于太平洋热带气旋和南海热带气旋登陆中国的主要路径上，极端天气引起的恶劣海况对基础附近海床的冲刷作用十分明显，随之形成的冲刷坑对当地风电场的稳定运行造成了严峻挑战。张冬冬[2]、杨博[3]等人通过不同范围和深度的冲刷坑对单桩风机结构进行了结构模态和响应模拟，发现冲刷坑的存在和发展会降低风机结构的自振频率，增加塔顶结构的水平位移，从而额外增加基础结构的负荷，降低结构整体的安全性。目前的风机基础设计中，对于由海流、波浪等环境载荷导致的基础冲刷问题，一般通过解析公式、缩比试验和数值模拟等方法对冲刷深度和冲刷坑形态进行预估。其中，解析公式能够实现对冲刷深度的快速预测，但计算结果难以通过精度分析的方式进行衡量，且计算结果无法给出冲刷坑的具体形态，因而也难以预测冲刷坑的存在对基础结构的动力学影响。缩比试验能够在给定的环境条件下，得到基础冲刷的真实数据，同时能够结合拟采用的防护措施进行防冲刷试验，可信度较高。但缩比试验所采用的环境条件较为单一，难以复现实际工程中的复杂流动环境，且由于试验条件的限制，实验结果往往受到边界条件的影响，因此缩比试验结果与工程实测数据仍然存在一定的差异。数值模拟能够克服缩比试验中的场地条件限制，通过数学模型实现复杂环境条件下的冲刷模拟，计算数据的可视化程度高，且计算结果具有较高的精度，因此能够弥补解析公式和缩比试验的不足。数值方法目前被广泛用于对冲刷问题开展细致地研究，以了解基础周围的流动环境、冲刷产生的机理和冲刷坑的发展规律，进而提出更合适的冲刷预测方法和防冲刷方案。
国内外对于桩基冲刷的数值研究开展较晚，主要集中于圆柱形单桩基础的冲刷模拟。Brilay和McDonald[4]首先通过三维不可压缩N-S方程研究了硬化床面上桩基础周围的流动结构。Olsen和Melaaen等人[5]使用稳态流动方程研究了无粘沙粒床面上，圆柱单桩的清水冲刷问题，床面切应力使用k－ε湍流模型求解，冲刷计算则采用泥沙输运方程。Roulund等人[6]实施了稳流冲刷水槽实验，在稳流条件下开展了硬化海床上圆柱形桩的定常和非定常动床冲刷数值模拟，对比了两种方法的计算结果；并将沙粒床面上的单桩基础冲刷数值模拟结果与实验数据进行了对比，在冲刷坑上游得到了较好的一致性，但数值模拟时间较长，需要在计算过程中逐步调整迭代时间步。Baykal[7]等人研究了波浪运动对单桩基础冲刷和海床回填作用的影响，在计算方法上，针对Roulund方法进行了改进，加入了悬移质冲刷算法，改善了桩后冲刷深度计算结果与实验值的一致性，但同样存在数值模拟时间较长的问题。
国内研究学者中，陈兵[8]等人使用切割单元法对海底管道冲刷进行了二维数值模拟，该方法在计算过程中不需要进行网格运动和重构操作，只需要施加合适的边界近似。数值计算结果与实验数据的对比显示，切割单元法的模型准确度高，数值结果符合较好。李绍武[9]等人使用Flow-3D软件对不同直径的圆桩冲刷进行了数值模拟，总结了冲刷深度与桩径的关系并研究了基础近壁面网格宽度对计算结果的影响。齐梅兰[10]等人通过三维数值模拟，研究了不同流速下，推移质和悬移质泥沙输运对冲刷坑发展的影响，提高了桩下游冲刷深度的计算精度。
本文基于海上风机单桩基础冲刷问题，在有限体积法求解框架下构建了流动-冲刷耦合计算模型。模型使用StarCCM+进行流动模拟，将床面切应力计算结果作为冲刷计算的输入，在自编程序中进行床面冲刷计算和计算域内网格的运动控制，随后将更新后的计算域返回到流动计算中进行下一时间步的计算，从而通过流动和冲刷的交叉迭代实现耦合问题的实时模拟。为解决冲刷计算周期过长的问题，在冲刷计算中引入了加速算法，在保证求解过程鲁棒性的基础上大幅缩短了计算周期。本文创新点在于构建了便于工程应用的桩基础冲刷数值模拟流程，在计算精度和资源消耗上实现了一定平衡，并通过算例验证了流程的可靠性。
-
流动计算采用三维非定常不可压缩N-S方程，湍流模型为SST k-ω模型，控制方程表达式如下：
连续性方程：
$$ \frac{{\partial {u_{i}}}}{{\partial {x_{i}}}} = 0 $$ (1) 动量方程：
$$ \frac{{\partial \rho {u_{i}}}}{{\partial t}} + \frac{{\partial \rho {u_{i}}{u_{\text{j}}}}}{{\partial {x_{j}}}} = \frac{\partial }{{\partial {x_{j}}}}\left( {\mu \frac{{\partial {u_{i}}}}{{\partial {x_{j}}}}} \right) - \frac{{\partial {p^ + }}}{{\partial {x_{i}}}} + {S_{i}} $$ (2) 湍动能方程：
$$ \begin{gathered} \frac{{\partial pk}}{{\partial t}} + \frac{{\partial pk{u_{j}}}}{{\partial {x_{j}}}} - \frac{\partial }{{\partial {x_{j}}}}\left[ {\left( {\mu + {\sigma _{k}}{\mu _{{\rm{T}}}}} \right)\frac{{\partial k}}{{\partial {x_{j}}}}} \right] = \\ {P_{k}} - \rho {\beta ^*}{f_{{{\beta }^*}}}\left( {\omega k - {\omega _0}{k_0}} \right) \\ \end{gathered} $$ (3) 耗散率方程：
$$ \begin{gathered} \frac{{\partial \rho \omega }}{{\partial t}} + \frac{{\partial \rho \omega {u_{j}}}}{{\partial {x_{j}}}} - \frac{\partial }{{\partial {x_{j}}}}\left[ {\left( {\mu + {\sigma _{\text{ω }}}{\mu _{\text{T}}}} \right)\frac{{\partial \omega }}{{\partial {x_{j}}}}} \right] = \\ {P_{\omega}} - \rho \beta {f_{\beta}}\left( {{\omega ^2} - \omega _0^2} \right) \\ \end{gathered} $$ (4) 式中：
$\rho $ ——清水密度（kg/m3）；$u$ ——流速矢量（m/s）；$\mu $ ——水的动力粘度（kg/m·s）；p+——包括静压在内的总压（Pa），
${{p}^{+}}=p+\rho gh $ ；${P_{k}} = {G_{k}} + {G_{{\rm{b}}}} + {G_{{{\rm{nl}}}}}$ 与${P_{\text{ω }}} = {G_{\text{ω }}} + {D_{\text{ω }}}$ 为结果项。湍流方程中各项的表达式和常数取值详见StarCCM+软件用户说明书。 -
冲刷运动主要受到床面切应力和泥沙粒径的影响。冲刷过程的每个时间步内，可以假设单元边界的床面切应力与单元内泥沙输运量存在平衡关系，泥沙输运处于守恒状态。因此海床高程的变化可以通过泥沙输运方程表示为：
$$ \frac{{\partial h}}{{\partial t}} = - \frac{1}{{1 - n}}\frac{{\partial {q_{{\text{b}i}}}}}{{\partial {x_{i}}}} $$ (5) 式中：
$n$ ——泥沙孔隙率，取0.4；$ i $ ——坐标系内的方向分量；${q_{{\text{b}i}}}$ ——推移质沿坐标系各方向的单宽输沙率矢量，由Engelund和Fredsøe公式得到：$$ {q_{\text{b}}} = \frac{1}{6}\text{π} {d^3}\frac{{{p_{{\text{EF}}}}}}{{{d^2}}}{u_{\text{b}}} $$ (6) 式中：
$ {\mu _{\text{d}}} $ ——动摩擦系数，取0.51。不难看出，单宽输沙率
${q_{{\text{bi}}}}$ 主要受到沙粒起动几率$ {p_{{\text{EF}}}} $ 和平均输沙速度矢量$ {u_b} $ 的影响，其中${p_{{\rm{EF}}}}$ 可表达为：$$ {p_{{\text{EF}}}} = {\left[ {1 + {{\left( {\frac{{\text{π} {\mu _{\text{d}}}}}{{6\left( {\theta - {\theta _{{\text{cr}}}}} \right)}}} \right)}^4}} \right]^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}}} $$ (7) 式中：
$\theta = {{{u_{\text{f}}}^2} \mathord{\left/ {\vphantom {{{u_{\text{f}}}^2} {\left( {s - 1} \right)}}} \right. } {\left( {s - 1} \right)}}{\rm{g}}d$ ——希尔兹数，表征静止海床上推移质抵抗剪切起动的能力；$s = {{{\rho _{\text{s}}}} \mathord{\left/ {\vphantom {{{\rho _{\text{s}}}} \rho }} \right. } \rho }$ ——推移质与清水的密度比；$ {u_{\text{f}}} = \sqrt {{\tau \mathord{\left/ {\vphantom {\tau \rho }} \right. } \rho }} $ ——基于壁面切应力$ \tau $ 的海床摩擦速度（cm/s）。任意床面坡度下的临界希尔兹数可通过下式计算：
$$ {\theta _{{\text{cr}}}} = {\theta _{{\text{c}}0}}\left( {\sqrt {1 - \frac{{{{\sin }^2}\alpha {{\tan }^2}\beta }}{{{\mu _{\text{s}}}^2}}} - \frac{{\cos \alpha \sin \beta }}{{{\mu _{\text{s}}}}}} \right) $$ (8) 式中：
$ {\mu _{\text{s}}} $ ——静摩擦系数，取0.63；$ {\theta _{{\text{c}}0}} $ ——水平海床条件下的临界希尔兹数，取0.05；α ——坡向和摩擦速度矢量的夹角（°），初始状态下为0；
β ——床面坡度（°）。
泥沙起动时，
$ {u_{\text{b}}} $ 的大小可由下式得到：$$ \frac{{{u_{\text{b}}}}}{{{u_{\text{f}}}}} = \alpha \left[ {1 - 0.7\sqrt {\frac{{{\theta _{{\text{cr}}}}}}{\theta }} } \right] $$ (9) 泥沙起动后，α的大小发生变化，无法通过解析方式求解，此时需要通过牛顿迭代法求解床面泥沙平衡方程组求解相关参数，具体的计算过程可参考文献[6]。
-
冲刷过程在每个时间段上，可以看作是海床的垂向运动，因此网格变化集中于垂直方向，在来流平面上无变化。本文根据单元边界处的输沙率矢量计算单元格心处的高度变化，格点处的高度变化则通过周围单元格心值的反距离加权得到：
$$ {\rm{d}}{H_n} = {{\sum\limits_{m = 1}^N {{w_m}{\rm{d}}{H_m}} } \mathord{\left/ {\vphantom {{\sum\limits_{m = 1}^N {{w_m}d{H_m}} } {\sum\limits_{m = 1}^N {{w_m}} }}} \right. } {\sum\limits_{m = 1}^N {{w_m}} }} $$ (10) 式中：
${\rm{d}}{H_n}$ ——格点高度变化（m）；${\rm{d}}{H_m}$ ——格心高度变化（m）；$ {w_m} = {1 \mathord{\left/ {\vphantom {1 {\left\| {{r_{mn}}} \right\|}}} \right. } {\left\| {{r_{mn}}} \right\|}} $ ——格心到格点距离的倒数。得到床面边界格点高度变化后，更新计算域底部边界，通过底部边界修改计算域内的格点位置，输入到流动计算模块中开始下一个迭代步的计算。
-
基础周围的床面遭受冲刷后形成冲刷坑，上游坑道的坡度往往大于下游坑道。对于无黏沙粒构成的海床，上游坑道坡度范围一般为28°～32°，一般情况下当地坡度不会超过该范围，因此也称为临界坡度。文献[6]在实验过程中发现，上游冲刷坑的部分区域内，坑道坡度存在短时间内超过临界值的情况，随后出现局部滑坡，且坑道坡度在重新稳定后回归正常值范围。数值模拟中同样需要考虑这一过程，因此采用韦雁机[11]等提出的方法，在完成全局冲刷计算后，搜索床面边界内存在局部坡度超过34°的单元，将单元节点沿水平轴翻转，直至当地坡度小于30°，随后重新形成完整的海床地貌。如图1所示，计算过程中，单元节点的坐标变化遵循下式：
Figure 1. Sand slide model
$$ \overrightarrow{{N_i}{N'_i}} = \frac{{\overrightarrow {OA} \times \overrightarrow {O{N_i}} }}{{\overrightarrow {OA} }} \cdot \frac{{\sin \left( {\beta - \phi } \right)}}{{\cos \phi }} $$ (11) 式中：
$ \phi $ ——滑坡后的稳定坡度（°），取30°；$ O $ ——单元形心；${N_{{i}}}$ ——滑坡前的单元格点；${N'_{{i}}}$ ——滑坡后的单元格点；A ——形心与格点连线向水平轴的垂线交点。
-
冲刷问题是流动和海床形变互相作用的耦合问题，本文采用弱耦合分步交叉进行求解，在数值模拟过程中需要考虑流动迭代时间步和冲刷时间步长的适配性。海床形变相较于水流流动是个缓慢发展的过程，如果冲刷时间步长与流动步长相同，会导致迭代次数大幅增加，计算周期十分漫长，从而失去工程应用的价值；若冲刷步长过大，不仅导致冲刷和流动各自的计算结果丧失相关性，还会因为冲刷深度变化过大而造成计算的发散。文献[6]在冲刷开始时限定冲刷计算步长，并随着冲刷的进行逐步增大步长。Liang[12]等人通过数值模拟发现，要兼顾耦合计算的稳定性和计算成本的经济性，需要每次迭代计算的深度变化不超过0.001D。本文采用文献[12]所述方法，首先给定冲刷计算的初始时间步长，在计算得到全场深度变化后，将全场最大坑深变化其放大至限定值并反求对应的时间步长，以此缩短耦合计算的次数，节约计算资源。
-
本文使用文献[6]中的实验作为算例，验证计算模型的真实性与可靠性。实验采用无黏沙粒，沙粒的50%中位数粒径d50=2.6 mm，计算时不考虑粒径分级，采用该数值作为代表。海床边界层厚度δ=0.2 m，圆桩直径D=0.1 m，近底平均流速Vref=0.46 m/s，边界层内速度廓线如图2所示，表示为：
Figure 2. Velocity profile
$$ \frac{u}{{{u_{\text{f}}}}} = 2.5{\rm{ln}}\frac{{30z}}{{{k_{\text{s}}}}} $$ (12) 式中：
z ——水线高度（m）；
ks——尼古拉斯平衡沙粒粗糙度（cm），取
$ 15 $ ；uf——海床摩擦速度（cm/s），取4.8。
为保证算例中计算域入口处的湍流充分发展，需要先在相同粗糙度的无扰动水平床面上，以廓线指定的入流速度分布进行试算，得到出口处流速垂向分布及湍流耗散情况，并将出口数据作为冲刷算例的入口条件加入求解器。
流动计算求解前，需要在StarCCM+中设置计算参数和计算条件。计算模型采用隐式非定常，柱面边界粗糙度设为2.5d50，壁面函数采用高Y+壁面处理，入口处的流速分布，湍动能和耗散率均从入口场函数中提取。迭代时间设为0.02 s，最大内迭代步设为10。
-
本文采用多块分区六面体网格的划分方式，计算域为沿流向18D，垂直流向15D，高度2D，域内包含12个网格分区，每个分区内的网格数量为323。桩和海床边界附近网格加密，边界附近的网格进行稀疏处理，网格具体形式如图3所示。对于边界条件的设置，入口条件设置为速度进口，出口条件指定为压力出口；考虑边界层顶部垂向运动对算例的影响较小，因此顶部及两侧边界条件设为对称边界条件，柱面及底面边界指定为无滑移壁面边界条件。
Figure 3. Detail of mesh for scour calculation Meshing
-
初始时刻的床面切应力计算结果如图4所示。通过与文献[6]中计算结果进行对比可以看出，最大区别在于本文流动计算中，海床切应力核心区沿桩体后移，靠近桩体沿流向的最大厚度位置，同时核心区的后移导致桩后涡街的运动范围受到挤压，因此桩侧剪切主导了桩周围的切应力分布。而文献[6]中的切应力核心区则较为靠前，接近桩前45°方向，且影响区域较为集中，与桩后涡街运动区域互不影响。这种区别的主要原因来自于软件和文献关于床面切应力计算方法和湍流相关计算项在公式选用上的差别。不同的海床切应力分布数据会导致冲刷计算结果产生一定差异，影响冲刷演化的进程，但由于冲刷演化是极为缓慢的过程，中间结果的微小扰动并不会对冲刷坑平衡形态产生太大的影响。
Figure 4. Seabed shear stress diagram
完成模型设置和网格划分后，对模型进行了用时两周的数值模拟，模拟结果如图5所示。模型计算时间达到16 min时停止计算，此时的冲刷深度已达到0.81D。
Figure 5. Scour pit development
结合图5所示的冲刷过程和图6所示的桩前坑深时变曲线，不难看出对于本算例中的桩基础，30 s时上游冲刷坑已经展现基本形态，且由于泥沙向下游堆积，桩两侧出现沙丘，沙丘的位置随着冲刷的进行也在不断后移，而沙丘背后伴随产生的凹坑则是由于沙丘对水流的遮挡导致水流在沙丘上游积蓄能量，到达下游后加剧了对局部床面的冲击。冲刷的前2 min内，冲刷速度最快，但随着坑深的增加和坑道范围的扩展，冲刷速度从5 min开始逐渐变慢，这是因为桩前冲刷坑在冲刷初始阶段尚未形成，此时桩前下洗涡和侧向马蹄涡的强度很大，对当地床面的冲刷效应十分显著。桩前冲刷坑逐渐形成后，坑深的增加和坑道坡度共同削弱了下洗涡和马蹄涡对床面的掏空作用，使得冲刷速度放缓。但坑道整体由于冲刷的积累仍然处于扩展状态，冲刷深度在16 min的计算时间内达到0.81D，说明桩土结构遭受了十分显著的破坏。
Figure 6. Time-varying curve of pit depth in front of pile
如图7所示的冲刷坑纵向对称轴剖面可以看出，坑道下游在桩后脱落涡的影响下逐渐成形，但由于上游泥沙的持续补给，下游坑深的变化始终小于上游坑深，坑道向下游发展，坑深缓慢增加，但变化速度远小于上游坑道。根据实验结果[6]，冲刷平衡状态下的桩下游坑深达到0.9D，稍低于上游1.2D的坑深。本文数值计算虽未达到平衡状态，但根据下游坑深发展趋势判断，桩下游的冲刷计算速度偏小。这一问题也出现在文献[13]中，原因是桩后冲刷由悬移质冲刷主导，泥沙在涡流起动作用下离开床面，随水流运动至下游。而对于推移质冲刷，由于桩体阻挡了水流的冲击，对桩后泥沙产生了遮挡效应。两方面原因共同造成了推移质模型下，桩后冲刷速度的计算结果偏小。
Figure 7. Longitudinal symmetry axis profile
-
本文通过结合CFD方法和泥沙输运模型，以不可压缩雷诺时均应力模型作为流动模块，以自编程序为海床冲刷及网格运动模块，建立了一套圆柱形单桩的实时冲刷动网格数值模拟流程。
1）通过本文模型与参考实验进行对比，验证了计算方法的可行性。
2）针对冲刷计算和流动计算在时间尺度上的不同，使用数值优化算法加快了迭代速度，并保证了算法的鲁棒性。
3）模拟过程中发现，桩基冲刷的最初阶段是冲刷发展最快，预防效果最好的时段，实际施工过程中应当考虑改进基础施工流程，在桩基入土后就开始冲刷的预防，从而更好地保护桩土结构的完整性。
4）后续工作将进行算法的改进，优化流动计算和网格运动算法，改进冲刷计算模块的程序结构，从而提高计算模型的整体精度，缩短预测时间。此外还将尝试扩展程序功能，使之具备多桩基础在浪、流荷载复杂流动条件下的数值模拟能力。
Numerical Simulation of Scour Process Around Offshore Wind Turbine Monopile Foundation in Steady Flow
doi: 10.16516/j.gedi.issn2095-8676.2023.01.010
- Received Date: 2021-10-27
- Rev Recd Date: 2021-11-18
- Available Online: 2022-11-29
- Publish Date: 2023-01-11
-
Key words:
- offshore wind farm /
- foundation scour /
- numerical simulation /
- nonlinear coupling /
- scour depth
Abstract:
| Citation: | XIA Bing, MIAO Desheng, ZHANG Min, LIU Huaixi, GE Wenpeng. Numerical Simulation of Scour Process Around Offshore Wind Turbine Monopile Foundation in Steady Flow[J]. SOUTHERN ENERGY CONSTRUCTION, 2023, 10(1): 81-87. doi: 10.16516/j.gedi.issn2095-8676.2023.01.010
|
DownLoad:
