CN102445318A - 一种顶张式立管顺流向振动分析方法 - Google Patents

Abstract

本发明涉及海洋深水立管的研究方法，具体涉及一种顶张式立管顺流向振动分析方法。该方法在立管的顺流向涡激振动分析中，考虑恒定拖曳力引起的立管动平衡位置偏移对立管顺流向振动的影响，从而使分析结果与立管的实际服役状态更吻合，并且考虑脉动拖曳力在锁定区和非锁定区的差异，改进现有顺流向振动分析模型，改善顺流向振动分析精度。

Description

一种顶张式立管顺流向振动分析方法

技术领域

本发明涉及海洋深水立管的研究方法，具体涉及一种顶张式立管顺流向振动分析方法。

背景技术

涡激振动是深水立管设计中的主要问题之一，由于我国的深水开发刚刚起步，尚没有深水立管的设计分析软件，目前均采用国外的同类型软件进行设计分析，如Orcaflex和Shear7等。

涡激振动是圆柱体在定常流作用下的一种奇特振动现象，是圆柱体尾流处的交替涡旋泄放引起的振动。对于刚性圆柱体和柔性不太大的圆柱体(一般的工程结构)，其垂直于流速方向(横流向)的振幅一般大于平行于流速方向(顺流向)的振幅，而振动频率为流速方向振动频率的二分之一。因此，现有的涡激振动研究一般忽略顺流向振动，而只研究横流向振动。现有的涡激振动分析方法只进行横流向的分析，如Iwan尾流振子模型：

$\frac{d^2{\stackrel{\‾}{y}}_n}{d{t}^2}+2{\mathrm{\ζ}}_n^T{\mathrm{\ω}}_n\frac{d{\stackrel{\‾}{y}}_n}{\mathrm{dt}}+{\mathrm{\ω}}_n^2{\stackrel{\‾}{y}}_n=\left(\frac{{\mathrm{\α}}_4}{v_n}\right)\frac{d{\stackrel{\‾}{z}}_n}{\mathrm{dt}}$ 式中：为结构横流向位移模态参数；

为尾流流体的模型参数。

和涡激升力模型：

式中：C_L--升力系数；

      ρ--流体密度；

      D--圆柱体直径；

U--流速；

ω_s--Strouhal频率，

其中：St为斯特罗哈数；

t--时间。

Orcaflex是目前被海洋工程界普遍接受的深水立管分析软件，它的涡激振动分析模块采用了尾流振子模型，计算结果与工程实测结果误差较大，所以，工程设计分析仅采用它的波浪分析功能，而涡激振动分析通常采用Shear7。

Shear7是深水立管涡激振动分析的专业软件，它采用工程实测数据修正计算参数，因此，计算结果与工程实测结果比较吻合，普遍为工程界所接受。Shear7采用的是频域分析方法。即：将立管的时域运动方程：

$m\frac{d^{2}y(z,t)}{d{t}^2}+c\frac{\mathrm{dy}(z,t)}{\mathrm{dt}}+\mathrm{EI}\frac{d^{4}y(z,t)}{d{z}^4}-T\frac{d^{2}y(z,t)}{d{z}^2}=\frac{1}{2}{\mathrm{cC}}_L\mathrm{\ρD}{U}^2\cos {\mathrm{\ω}}_{s}t$

式中：m--立管单位长度的质量；

      y(z，t)--立管的横流向弯曲位移，其中，z为立管轴向坐标；

      c--阻尼系数，包括水动力阻尼；

      EI--立管截面弯曲刚度；

      T--立管有效张力。

利用坐标变换：

$y(z,t)=\underset{r}{\mathrm{\Σ}}{Y}_r\left(z\right)q_r\left(t\right)$

式中：Y_r(z)--立管第r阶振型；

      q_r(t)--立管第r阶振型的坐标。

转化为模态(振型)坐标方程：

$M_r{\stackrel{..}{q}}_r\left(t\right)+C_r{\stackrel{.}{q}}_r\left(t\right)+K_r{q}_r\left(t\right)=P_r\left(t\right)(r=\mathrm{1,2},\·\·\·,n)$

式中：M_r--立管第r阶模态质量， $m_r={\∫}_0^L{Y}_r^2\left(z\right)\mathrm{mdz};$

      C_r--立管第r阶模态阻尼， $C_r={\∫}_0^L{Y}_r^2\left(z\right)(c+c_a)\mathrm{dz};$

      K_r--立管第r阶模态刚度，

$K_r={\∫}_0^L\mathrm{EI}{\left(\frac{d^2{Y}_r\left(z\right)}{{\mathrm{dz}}^2}\right)}^2\mathrm{dz}-{\∫}_0^{L}T{\left(\frac{d{Y}_r\left(z\right)}{\mathrm{dz}}\right)}^2\mathrm{dz}$

P_s(t)--第r阶模态力，

$P_r\left(t\right)={\∫}_0^L\frac{1}{2}{C}_L\mathrm{\ρD}{U}^2{Y}_r\left(z\right)\mathrm{dz}\sin {\mathrm{\ω}}_{s}t$

然后采用随机振动的谱分析方法(公知技术)求解上述模态(振型)坐标方程，即采用傅利叶变换将时域的力函数转换为频域的谱函数，从而求得立管的响应谱。因此，Shear7的计算结果是以响应的平均值和均方差表示的，不能给出响应时程，因而不能给出立管响应的最大值。

上述两个软件采用的方法是目前比较有代表性的两种涡激振动分析方法，它们都没有考虑立管的顺流向振动问题，完全基于传统的圆柱体涡激振动理论和分析方法。而深水立管由于柔性较大，其涡激振动为高阶(＞10阶)多模态(5～20个模态)振动(一般圆柱体的涡激振动为1～3阶模态振动)，研究表明，深水立管的顺流向涡激振动幅值与横流向振动幅值之比大于50％。国外也有学者提出深水立管涡激振动分析应考虑顺流向振动问题。

近年来。关于立管顺流向振动的研究越来越多，但尚没有考虑定常流引起的恒定拖曳力对振动的影响。恒定拖曳力将引起立管振动平衡位置的偏移，从而造成动应力的均值不为零。而动应力的均值对结构的疲劳寿命有较大的影响，疲劳理论中有三种应力均值对疲劳损伤的修正方法：

1、Goodman直线模型

${\mathrm{\σ}}_a={\mathrm{\σ}}_{-1}[1-\left(\frac{{\mathrm{\σ}}_m}{{\mathrm{\σ}}_b}\right)]$

2、Gerber抛物线模型

${\mathrm{\σ}}_a={\mathrm{\σ}}_{-1}[1-{\left(\frac{{\mathrm{\σ}}_m}{{\mathrm{\σ}}_b}\right)}^2]$

3、Soderberg直线模型

${\mathrm{\σ}}_a={\mathrm{\σ}}_{-1}[1-\left(\frac{{\mathrm{\σ}}_m}{{\mathrm{\σ}}_s}\right)]$

式中：σ_a--循环应力幅；

      σ_-1--对称循环应力幅；

      σ_m--平均应力。

与本发明最相近的解决方案是申请人的前一项发明专利“一种深水顶张式立管涡激振动与疲劳分析的方法”ZL200910136583.X，该专利提出了一种涡激振动分析方法包括顺流向振动分析，如下式，但是并没有考虑恒定拖曳力的影响，也没有考虑脉动拖曳力在锁定区和非锁定区的区别。

$(m+m_a)\frac{{\∂}^{2}x(z,t)}{\∂t^2}+(c+c_a)\frac{\∂x(z,t)}{\∂t}$

$+\mathrm{EI}\frac{{\∂}^{4}x(z,t)}{\∂z^4}-T\frac{{\∂}^{2}x(z,t)}{\∂z^2}=\frac{1}{2}{C}_D\mathrm{\ρD}{(U-\stackrel{.}{x})}^2\cos 2{\mathrm{\ω}}_s^{\′}t$

式中：x(z，t)--顺流向振动响应函数；

      m_a--附加质量系数；

      c_a--附加阻尼系数；

      C_D--拖曳力系数；

      ω′_s--考虑流固耦合的涡旋泄放频率，

总的来说，现有技术存在如下缺陷：

1、没有考虑恒定拖曳力引起的动平衡位置偏移，因此，位移响应和应力响应均与立管的实际值有一定的误差，计算结果偏于不安全；

2、没有考虑脉动拖曳力在锁定区和非锁定区的差别，而全部采用锁定区的计算方法，与立管的实际情况有较大区别。因为，立管的设计会避开锁定区，即在服役海域的主流速条件下，立管的涡激振动将处于非锁定区。

发明内容

本发明的目的在于针对现有技术的缺陷，提供一种顶张式立管顺流向振动分析方法，在立管的顺流向涡激振动分析中，考虑恒定拖曳力引起的立管动平衡位置偏移对立管顺流向振动的影响，从而使分析结果与立管的实际服役状态更吻合，并且考虑脉动拖曳力在锁定区和非锁定区的差异，改进现有顺流向振动分析模型，改善顺流向振动分析精度。

本发明的技术方案如下：一种顶张式立管顺流向振动分析方法，采用了同时考虑脉动拖曳力在锁定区和非锁定区的差异和恒定拖曳力的顺流向振动分析模型：

a)非锁定区

$(m+m_a)\stackrel{..}{u}+c\stackrel{.}{u}+{\mathrm{EIu}}^{\mathrm{IV}}-{\mathrm{Tu}}^{\′\′}=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t]$

b)锁定区

$(m+m_a)\stackrel{..}{u}+c\stackrel{.}{u}+{\mathrm{EIu}}^{\mathrm{IV}}-{\mathrm{Tu}}^{\′\′}=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t]$

式中，u--顺流向振动位移响应；

--顺流向振动速度响应；

      ü--顺流向振动加速度响应；

      u^IV--u对x的4次导数，x为立管的轴向坐标；

      u＂--u对x的2次导数；

      m--立管单位长度的质量；

      m_a--附加质量系数，顶张式立管表面流体随立管振动的惯性力，

$m_a=\mathrm{\ρ}\frac{\mathrm{\π}{D}^2}{4};$

      ρ--流体密度；

      D--顶张式立管直径；

      U--流速；

      EI--立管截面弯曲刚度；

      T--立管有效张力；

      C_D--恒定拖曳力系数，通常取1.2；

      C′_D--脉动拖曳力系数，通常取0.7；

     ω′_s--考虑流固耦合的涡旋泄放频率，

其中，St为斯托哈尔数，通常取0.2；

t--时间；

针对顶张式立管顺流向振动速度和加速度，采用迭代方法对以上公式进行计算，分别得到非锁定区和锁定区的顺流向响应。

进一步，如上所述的顶张式立管顺流向振动分析方法，其中，所述的采用迭代方法进行计算的具体步骤如下：

1)给定顶张式立管顺流向振动速度和加速度及计算时间的初值：

t_j＝t₀＝0， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_0^{\left(0\right)}=0$

式中：j--时间步数，计算开始时j＝0；

      i--迭代次数，每个时间步开始时i＝0；

2)计算给定流速下的约化速度：

$V_r=\frac{U}{f_{n}D}$

式中：V_r--约化速度；

      U--流体流速；

      f_n--顶张式立管的固有频率；

      D--顶张式立管直径；

3)当V_r＜5或V_r＞7时，将t_j，

代入下式计算第j时间步内第i+1次迭代的顺流向振动响应：

$(m+m_a){\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+\mathrm{EI}{\left(u^{\mathrm{IV}}\right)}_j^{(i+1)}-T{\left(u^{\′\′}\right)}_j^{(i+1)}$

$=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}t]$

当5≤V_r≤7时，将t_j，

代入下式计算第j时间步内第i+1次迭代的顺流向振动响应：

$(m+m_a){\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+\mathrm{EI}{\left(u^{\mathrm{IV}}\right)}_j^{(i+1)}-T{\left(u^{\′\′}\right)}_j^{(i+1)}$

$=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^{2}2\cos {\mathrm{\ω}}_s^{\′}t]$

4)如果

ε为预先设定的计算精度，则继续进行迭代计算，令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_j^{(i+1)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_j^{(i+1)}$

然后，重复步骤3)～4)的计算；

如果

则开始下一个时间步的计算，令：

t_j＝t_j+1＝t_j+Δt， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_{j+1}^{\left(0\right)}={\stackrel{.}{u}}_j^{\left(n\right)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_{j+1}^{\left(0\right)}={\stackrel{..}{u}}_j^{\left(n\right)}$

n为第j时间步的总迭代次数；

重复步骤3)～4)的计算，直至计算时长满足需要。

本发明的有益效果如下：本发明采用了考虑恒定拖曳力的立管顺流向振动分析模型，采用了迭代技术求解立管的顺流向振动方程。与现有顺流向振动分析方法相比具有下述优点：

(1)考虑了恒定拖曳力引起的立管顺流向动平衡位置变化，比现有方法更符合立管的顺流向振动和受力状态。

(2)考虑了脉动拖曳力在锁定区和非锁定区的差异，采用了立管在锁定区或非锁定区的顺流向振动方程，比现有技术仅采用锁定区的顺流向振动方程更符合立管顺流向涡激振动的本质特征。

附图说明

图1为本发明的方法流程图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细描述。

顶张式立管在定常流场作用下将发生涡激振动，其中的顺流向振动是由于涡旋脱落引起的拖曳力波动造成的。由于拖曳力的脉动量远远小于拖曳力的定常值，因此，立管顺流向振动的平衡位置不是立管的静平衡位置(流速为零时的位置)。现有技术的立管顺流向振动分析是以静平衡位置为振动的平衡位置，从而使得动应力的均值为零。而立管顺流向振动的动应力均值不为零，且应力均值对立管的疲劳寿命有较大的影响，因此，现有技术的计算结果将造成较大的疲劳计算误差，结果偏于不安全。本发明建立了一种立管顺流向振动的分析方法，使分析结果更符合实际。

本发明采用了同时考虑脉动拖曳力在锁定区和非锁定区的差异和恒定拖曳力的顺流向振动分析模型：

a)非锁定区

$(m+m_a)\stackrel{..}{u}+c\stackrel{.}{u}+{\mathrm{EIu}}^{\mathrm{IV}}-{\mathrm{Tu}}^{\′\′}=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t]---\left(1\right)$

b)锁定区

$(m+m_a)\stackrel{..}{u}+c\stackrel{.}{u}+{\mathrm{EIu}}^{\mathrm{IV}}-{\mathrm{Tu}}^{\′\′}=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t]---\left(2\right)$

式中，u--顺流向振动位移响应；

--顺流向振动速度响应；

      ü--顺流向振动加速度响应；

      u^IV--u对x的4次导数，x为立管的轴向坐标；

      u＂--u对x的2次导数；

      m--立管单位长度的质量；

      m_a--附加质量系数，顶张式立管表面流体随立管振动的惯性力，

$m_a=\mathrm{\ρ}\frac{\mathrm{\π}{D}^2}{4};$

      ρ--流体密度；

      D--顶张式立管直径；

      U--流速；

      EI--立管截面弯曲刚度；

      T--立管有效张力；

      C_D--恒定拖曳力系数，通常取1.2；

      C′_D--脉动拖曳力系数，通常取0.7；

      ω′_s--考虑流固耦合的涡旋泄放频率，

其中，St为斯托哈尔数，通常取0.2；

t--时间；公式(1)和(2)分别用于非锁定区和锁定区的顺流向振动分析。

由于公式(1)和公式(2)为隐式方程，因此，必须采用迭代方法(公知方法)计算，如图1所示，具体计算步骤如下：

1)给定顶张式立管顺流向振动速度和加速度及计算时间的初值(公知技术)：

t_j＝t₀＝0， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_0^{\left(0\right)}=0$

式中：j--时间步数，计算开始时j＝0；

i--迭代次数，每个时间步开始时i＝0；

2)计算给定流速(流速是设计或分析给定的条件，为已知值)下的约化速度(公知技术)：

$V_r=\frac{U}{f_{n}D}$

式中：V_r--约化速度；

      U--流体流速；

      f_n--顶张式立管的固有频率；

      D--顶张式立管直径；

3)当V_r＜5或V_r＞7时，将t_j，

代入下式计算第j时间步内第i+1次迭代的顺流向振动响应(包括位移、速度、加速度)：

$(m+m_a){\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+\mathrm{EI}{\left(u^{\mathrm{IV}}\right)}_j^{(i+1)}-T{\left(u^{\′\′}\right)}_j^{(i+1)}$

$=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}t]$

当5≤V_r≤7时，将t_j，

代入下式计算第j时间步内第i+1次迭代的顺流向振动响应(包括位移、速度、加速度)：

$(m+m_a){\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+\mathrm{EI}{\left(u^{\mathrm{IV}}\right)}_j^{(i+1)}-T{\left(u^{\′\′}\right)}_j^{(i+1)}$

$=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^{2}2\cos {\mathrm{\ω}}_s^{\′}t]$

上述两个公式分别对应非锁定区和锁定区的顺流向振动分析模型，各参数的含义与上面公式(1)、(2)中对应参数的含义相同。

4)如果

，ε为预先设定的计算精度(根据需要设定，如ε＝1×10^-5)，则继续进行迭代计算，令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_j^{(i+1)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_j^{(i+1)}$

然后，重复步骤3)～4)的计算；

如果

则开始下一个时间步的计算，令：

t_j＝t_j+1＝t_j+Δt， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_{j+1}^{\left(0\right)}={\stackrel{.}{u}}_j^{\left(n\right)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_{j+1}^{\left(0\right)}={\stackrel{..}{u}}_j^{\left(n\right)}$

n为第j时间步的总迭代次数；Δt一般取0.02秒；

重复步骤3)～4)的计算，直至计算时长满足需要(时长根据需要设定，如100秒)。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若对本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (2)

1.一种顶张式立管顺流向振动分析方法，其特征在于：采用了同时考虑脉动拖曳力在锁定区和非锁定区的差异和恒定拖曳力的顺流向振动分析模型：

a)非锁定区

$(m+m_a)\stackrel{..}{u}+c\stackrel{.}{u}+{\mathrm{EIu}}^{\mathrm{IV}}-{\mathrm{Tu}}^{\′\′}=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t]$

b)锁定区

$(m+m_a)\stackrel{..}{u}+c\stackrel{.}{u}+{\mathrm{EIu}}^{\mathrm{IV}}-{\mathrm{Tu}}^{\′\′}=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t]$

式中，u--顺流向振动位移响应；

--顺流向振动速度响应；

ü--顺流向振动加速度响应；

      u^IV--u对x的4次导数，x为立管的轴向坐标；

      u＂--u对x的2次导数；

      m--立管单位长度的质量；

      m_a--附加质量系数，顶张式立管表面流体随立管振动的惯性力，

$m_a=\mathrm{\ρ}\stackrel{\mathrm{\π}}{-}\frac{D^2}{4};$

      ρ--流体密度；

      D--顶张式立管直径；

      U--流速；

      EI--立管截面弯曲刚度；

      T--立管有效张力；

      C_D--恒定拖曳力系数；

      C′_D--脉动拖曳力系数；

      ω′_s--考虑流固耦合的涡旋泄放频率，

其中，St为斯托哈尔数；

t--时间；

针对顶张式立管顺流向振动位移、速度和加速度，采用迭代方法对以上公式进行计算，分别得到非锁定区或锁定区的顺流向响应。

2.如权利要求1所述的顶张式立管顺流向振动分析方法，其特征在于：所述的采用迭代方法进行计算的具体步骤如下：

1)给定顶张式立管顺流向振动速度和加速度及计算时间的初值：

t_j＝t₀＝0， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_0^{\left(0\right)}=0$

式中：j--时间步数，计算开始时j＝0；

i--迭代次数，每个时间步开始时i＝0；

2)计算给定流速下的约化速度：

$V_r=\frac{U}{f_{n}D}$

式中：V_r--约化速度；

      U--流体流速；

      f_n--顶张式立管的固有频率；

       D--顶张式立管直径；

3)当V_r＜5或V_r＞7时，将t_j，代入下式计算第j时间步内第i+1次迭代的顺流向振动响应：

$(m+m_a){\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+\mathrm{EI}{\left(u^{\mathrm{IV}}\right)}_j^{(i+1)}-T{\left(u^{\′\′}\right)}_j^{(i+1)}$

$=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}t]$

当5≤V_r≤7时，将t_j，代入下式计算第j时间步内第i+1次迭代的顺流向振动响应：

$(m+m_a){\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+\mathrm{EI}{\left(u^{\mathrm{IV}}\right)}_j^{(i+1)}-T{\left(u^{\′\′}\right)}_j^{(i+1)}$

$=\frac{1}{2}\mathrm{\ρD}[C_D{U}^2+C_D^{\′}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^{2}2\cos {\mathrm{\ω}}_s^{\′}t]$

4)如果

ε为预先设定的计算精度，则继续进行迭代计算，令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_j^{(i+1)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_j^{(i+1)}$

然后，重复步骤3)～4)的计算；

如果则开始下一个时间步的计算，令：

t_j＝t_j+1＝t_j+Δt， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_{j+1}^{\left(0\right)}={\stackrel{.}{u}}_j^{\left(n\right)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_{j+1}^{\left(0\right)}={\stackrel{..}{u}}_j^{\left(n\right)}$

n为第j时间步的总迭代次数；

重复步骤3)～4)的计算，直至计算时长满足需要。

Images