水下滑翔机水动力系数数值模拟研究

Numerical Simulation of Hydrodynamic Co-efficient of Underwater Gliders

摘要: 采用三维Navier-Stokes方程开展了三维回转椭球体的垂荡及纵摇水动力系数试验数值模拟研究，在此基础上进一步开展了SUBOFF模型潜艇横荡及艏摇试验的数值仿真。分析了水下滑翔机垂荡、纵摇、横荡及艏摇实验的运动形式并推导了对应的数学方程，同时开展了水下滑翔机水动力系数数值预报仿真。结果表明，本文发展的水下滑翔机水动力系数数值模拟方法计算结果与实验结果对比一致，相关推导过程及计算方法对水动力系数数值预报方法有一定的参考借鉴意义。

Abstract: The numerical simulation of the vertical and longitudinal hydrodynamic coefficient test of the three-dimensional rotary ellipsoid was carried out by using the three-dimensional Navier-Stokes equation, and the numerical simulation of the transverse and bow swing test of the SUBOFF model submarine was further carried out on this basis. The motion forms of the underwater glider hanging, longitudinal, transverse and bow rocking experiments are analyzed, and the corresponding mathematical equations are derived, and the numerical prediction simulation of the hydrodynamic coefficient of the underwater glider is carried out. The results show that the calculation results of the numerical simulation method of hydraulic coefficient of underwater gliders developed in this paper are consistent with the experimental results, and the relevant derivation process and calculation method have certain reference significance for the numerical prediction method of hydrodynamic coefficient.

1. 引言

水下滑翔机的垂荡试验是一种典型的垂直面平面运动机构试验，该机构由振荡机构、驱动电控系统和测力数据处理系统等三部分组成。在拖车匀速直行的同时，调节两支振荡杆做正弦振荡，保持同位相、同振幅和同频率，这样滑翔机在垂直面内作正弦运动，其纵倾角始终为零，综合两种运动就实现了垂荡运动。水下滑翔机纵摇运动也是垂直面较为典型的运动之一，纵摇运动是在垂荡运动的基础上增加了绕轴转动，增加了角速度q，但是模型在随体坐标系中垂向速度与加速度为0。作为水下滑翔机在水平面的典型运动，横荡及艏摇运动试验及其数值模拟是研究操纵性能水动力的重要部分 [1] 。

本文采用三维Navier-Stokes方程开展了三维回转椭球体的垂荡及纵摇水动力系数试验数值模拟研究，在此基础上进一步开展了SUBOFF模型潜艇横荡及艏摇试验的数值仿真。分析了水下滑翔机垂荡、纵摇、横荡及艏摇实验的运动形式并推导了对应的数学方程，同时开展了水下滑翔机水动力系数数值预报仿真。

2. 垂荡试验数值模拟验证

2.1. 垂荡运动及数学方程

对于垂荡运动，有：

${\xi }_1={\xi }_2=a\sin \left(\omega t\right)$ (1)

式中：

${\xi }_1,{\xi }_2$ ——两支杆端部的垂向位移，也即滑翔机的垂向位移。

$\omega$ ——偏心轮的角速度，即支杆的振荡频率，也即滑翔机垂荡运动的圆频率。

a——滑翔机垂荡运动的振幅。

此时滑翔机的运动参数为：

$\{\begin{array}{l}\theta =\stackrel{\dot{}}{\theta }=0\\ w=\stackrel{\dot{}}{\xi }=a\omega \cos \left(\omega t\right)\\ \stackrel{\dot{}}{w}=-a{\omega }^2\sin \left(\omega t\right)\end{array}$ (2)

式中：

$\theta ,\stackrel{\dot{}}{\theta }$ ——滑翔机绕y轴的倾斜角度和角速度。

$w,\stackrel{\dot{}}{w}$ ——滑翔机的垂向速度和加速度。

设在运动过程中滑翔机受到的垂向力为Z，绕Y轴的力矩为纵倾力矩M，滑翔机在Z、M作用下，做垂荡运动。垂荡运动是模拟滑翔机做弱机动运动，即运动参数w，a等均为小量，同时由于仅模拟纯升沉运动，所以角速度q等参数为0。针对水动力模型并结合水动力系数的意义 [2] ，得到简化的线性运动方程为：

$\{\begin{array}{l}Z=Z_{\stackrel{\dot{}}{w}}\stackrel{\dot{}}{w}+Z_{w}w+Z_0=-a{\omega }^2{Z}_{\stackrel{\dot{}}{w}}\sin \left(\omega t\right)+a\omega Z_w\cos \left(\omega t\right)+Z_0\\ M=M_{\stackrel{\dot{}}{w}}\stackrel{\dot{}}{w}+M_{w}w+M_0=-a{\omega }^2{M}_{\stackrel{\dot{}}{w}}\sin \left(\omega t\right)+a\omega M_w\cos \left(\omega t\right)+M_0\end{array}$ (3)

将(3)式进行无量纲化，垂向力方程遍除 $0.5\rho V^2{L}^2$ ，纵摇力矩方程遍除 $0.5\rho V^2{L}^3$ ，得到下式：

$\{\begin{array}{l}{Z}^{\prime }=-\frac{a{\omega }^{2}L}{V^2}{Z^{\prime }}_{\stackrel{\dot{}}{w}}\sin \left(\omega t\right)+\frac{a\omega }{V}{Z^{\prime }}_w\cos \left(\omega t\right)+{Z^{\prime }}_0\\ M^{\prime }=-\frac{a{\omega }^{2}L}{V^2}{M^{\prime }}_{\stackrel{\dot{}}{w}}\sin \left(\omega t\right)+\frac{a\omega }{V}{M^{\prime }}_w\cos \left(\omega t\right)+{M^{\prime }}_0\end{array}$ (4)

式中：L为参考长度，V为来流速度。

采用三维回转椭球体开展垂荡实验的数值模拟验证。根据试验情况，数值计算时入口处的来流为V = 0.8 m/s，给定正弦振荡振幅α = 2.5 m，模拟工况振荡频率f为0.01 Hz。椭球体在均匀流场中保持正弦振荡运动，这样合成的结果即为椭球体的升沉运动，求解参数设置及边界条件与斜航实验保持一致。

2.2. 水动力系数求解分析

针对三维回转椭球体的垂荡实验，采用积分法求解椭球的水动力系数。令：

$\{\begin{array}{l}{Z}_a=-\frac{a{\omega }^{2}L}{V^2}{Z^{\prime }}_{\stackrel{\dot{}}{w}}{Z}_b=\frac{a\omega }{V}{Z^{\prime }}_w\\ M_a=-\frac{a{\omega }^{2}L}{V^2}{M^{\prime }}_{\stackrel{\dot{}}{w}}{M}_b=\frac{a\omega }{V}{M^{\prime }}_w\end{array}$ (5)

将(5)代入(4)得到：

$\{\begin{array}{l}{Z}^{\prime }=Z_a\sin \left(\omega t\right)+Z_b\cos \left(\omega t\right)+{Z^{\prime }}_0\\ M^{\prime }=M_a\sin \left(\omega t\right)+M_b\cos \left(\omega t\right)+{M^{\prime }}_0\end{array}$ (6)

采用积分法求解系数 $Z_a,Z_b,M_a,M_b$ 。

记 $Z^{\prime }$ 随时间变化的函数为 $f\left(t\right)$ ，并记垂荡运动的周期为2l，由傅立叶级数定理展开周期为2l的周期函数，得到：

$f\left(t\right)=\frac{a_0}{2}+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left(a_n\cos \frac{n\pi t}{l}+b_n\sin \frac{n\pi t}{l}\right)}$ (7)

其中系数 $a_n,b_n$ 为：

$\{\begin{array}{l}{a}_n=\frac{1}{l}{\displaystyle {\int }_{-1}^{1}f\left(t\right)\cos \frac{n\pi t}{l}\text{d}t\left(n=0,1,2,\cdots \right)}\\ b_n=\frac{1}{l}{\displaystyle {\int }_{-1}^{1}f\left(t\right)\sin \frac{n\pi t}{l}\text{d}t\left(n=1,2,3,\cdots \right)}\end{array}$ (8)

根据水动力系数的含义， $Z_a,Z_b$ 的表达式如下：

$\{\begin{array}{l}{Z}_a=b_1=\frac{1}{l}{\displaystyle {\int }_{-l}^{l}f\left(t\right)}\sin \left(\omega t\right)\text{d}t\\ Z_b=a_1=\frac{1}{l}{\displaystyle {\int }_{-l}^{l}f\left(t\right)}\cos \left(\omega t\right)\text{d}t\end{array}$ (9)

采用同样的方法可以得到 $M_a,M_b$ 的值。

计算得到的水动力系数计算值与理论值对比如表1所示。垂荡试验计算值、斜航试验计算值以及理论值差别不大。较好的验证了垂荡试验数值模拟方法及水动力系数公式推导的正确性。

Table 1. Comparison of the calculated hydrodynamic coefficient of the vertical test with the theoretical value

表1. 垂荡试验水动力系数计算值与理论值对比

3. 纵摇试验数值模拟验证

3.1. 纵摇运动及数学方程

纵摇运动能够获得关于角速度q与角加速度 $\stackrel{\dot{}}{q}$ 的水动力系数。

纵摇运动时，滑翔机运动速度与重心轨迹曲线相切，滑翔机冲角、滑翔机动坐标系的Z向速度与加速度均为零，即 $\alpha =w=\stackrel{\dot{}}{w}=0$ 。当保持前后支杆的振幅相等，调节其相位差，使后杆对于前杆有一定的滞后角 $\epsilon$ ：

$\epsilon =2\arctan \frac{l_0\omega }{V}$ (10)

此时两支杆的位移各为：

$\{\begin{array}{l}{\xi }_1=\alpha \cos \left(\omega t+\frac{\epsilon }{2}\right)\\ {\xi }_2=\alpha \cos \left(\omega t-\frac{\epsilon }{2}\right)\end{array}$ (11)

此时滑翔机的运动规律如下：

$\{\begin{array}{l}\theta ={\theta }_0\sin \left(\omega t\right)\\ q=\stackrel{\dot{}}{\theta }={\theta }_0\omega \cos \left(\omega t\right)\\ \stackrel{\dot{}}{q}=-{\theta }_0{\omega }^2\sin \left(\omega t\right)\\ w=\stackrel{\dot{}}{w}=0\end{array}$ (12)

式中：

$\theta$ ——滑翔机运动过程中的纵倾角。

${\theta }_0$ ——滑翔机的最大纵倾角。

$\omega$ ——偏心轮的角速度，即支杆的振荡频率，也即滑翔机纵摇运动的圆频率。

$q,\stackrel{\dot{}}{q}$ ——滑翔机的纵摇角速度和纵摇角加速度。

设运动过程中滑翔机受到的垂向力为Z，绕Y轴的力矩为纵倾力矩M。根据水动力模型，并考虑到纵摇运动时唯一的机动参数为角速度q [3] ，可得到纵摇运动线性运动方程：

$\{\begin{array}{l}Z=Z_{\stackrel{\dot{}}{q}}\stackrel{\dot{}}{q}+Z_{q}q+Z_0=-{\theta }_0{\omega }^2{Z}_{\stackrel{\dot{}}{q}}\sin \left(\omega t\right)+{\theta }_0\omega Z_q\cos \left(\omega t\right)+Z_0\\ M=M_{\stackrel{\dot{}}{q}}\stackrel{\dot{}}{q}+M_{q}q+M_0=-{\theta }_0{\omega }^2{M}_{\stackrel{\dot{}}{q}}\sin \left(\omega t\right)+{\theta }_0\omega M_q\cos \left(\omega t\right)+M_0\end{array}$ (13)

将(13)式进行无量纲化，垂向力方程遍除 $0.5\rho V^2{L}^2$ ，纵摇力矩方程遍除 $0.5\rho V^2{L}^3$ ，得到下式：

$\{\begin{array}{l}{Z}^{\prime }={Z^{\prime }}_{\stackrel{\dot{}}{q}}\stackrel{\dot{}}{q}+{Z^{\prime }}_{q}q+{Z^{\prime }}_0\\ M^{\prime }={M^{\prime }}_{\stackrel{\dot{}}{q}}\stackrel{\dot{}}{q}+{M^{\prime }}_{q}q+{M^{\prime }}_0\end{array}$ (14)

采用三维回转椭球体开展纵摇实验的数值模拟验证。根据试验情况，数值计算时入口处的来流为V = 0.8 m/s，给定正弦振荡振幅 ${\theta }_0=\text{1}\text{rad}$ ，模拟工况振荡频率f为0.1 Hz。椭球体在均匀流场中保持正弦振荡运动，这样合成的结果即为椭球体的纵摇运动，求解参数设置及边界条件与斜航实验保持一致。

3.2. 水动力系数求解分析

采用频域变换法求解水动力系数 ${Z^{\prime }}_q,{Z^{\prime }}_{\stackrel{\dot{}}{q}},{M^{\prime }}_q,{M^{\prime }}_{\stackrel{\dot{}}{q}}$ 。频域分析有许多优势并在不同的领域中得到应用。通过作用在稳定谐波的一个周期上的频域变换可以获得水动力的频谱特性及其同相/异相分量。一个连续的时间标量函数 $x\left(t\right)$ 在一个有限的时间间隔内 $t\in \left[0,T\right]$ 的有限傅里叶积分定义为：

$\mathbb{F}\left[x\left(t\right)\right]\equiv \bar{x}\left(i\omega \right)={\displaystyle {\int }_0^{T}x\left(t\right){\text{e}}^{-i2\pi ft}\text{d}t}$ (15)

其中i是虚数单位，f是空间频率。

有限傅里叶积分可以在离散的频率f_k上得到计算。f_k等间隔的从零频率到奈奎斯特频率f_N，其频率分辨率为时间步长的倒数。给定波段上的任意频率求解可以采用线性调频Z变换得到。系统响应中与输入信号密切相关的振幅比和相位角可以用输入和输出的传递函数来确定。对垂向力Z，其传递函数为：

$G\left(i\omega \right)=\frac{\mathbb{F}\left[Z^{\prime }\left(t\right)\right]}{\mathbb{F}\left[q^{\prime }\left(t\right)\right]}=R\left(\omega \right){\text{e}}^{i\varphi \left(\omega \right)}$ (16)

其中 $R\left(\omega \right)$ 和 $\varphi \left(\omega \right)$ 分别表示振幅比和相位角，其定义为：

$R\left(\omega \right)=\frac{\Vert {\bar{Z}}^{\prime }\left(i\omega \right)\Vert }{\Vert \overline{q^{\prime }}\left(i\omega \right)\Vert }$ (17)

$\varphi \left(\omega \right)=\angle {\bar{Z}}^{\prime }\left(i\omega \right)-\angle \overline{q^{\prime }}\left(i\omega \right)$ (18)

从而求得水动力系数：

${Z^{\prime }}_q=R\left(\omega \right)\cos \varphi \left(\omega \right)$ (19)

${Z^{\prime }}_{\stackrel{\dot{}}{q}}=R\left(\omega \right)\sin \varphi \left(\omega \right)/\omega$ (20)

采用同样方法可以求得 ${M^{\prime }}_q,{M^{\prime }}_{\stackrel{\dot{}}{q}}$ 。

计算得到的水动力系数计算值与理论值对比如表2所示。表2中还给出了其他文献的计算结果。本文计算结果优于其他文献的计算结果，部分文献的计算甚至出现了垂向力系数导数反号的现象。本文纵摇试验计算值与理论值差别不大。较好的验证了纵摇试验数值模拟方法及水动力系数公式推导的正确性。

Table 2. Comparison of the calculated value of the hydrodynamic coefficient of the longitudinal shaking test with the theoretical value

表2. 纵摇试验水动力系数计算值与理论值对比

4. 横荡试验数值模拟验证

4.1. 横荡运动及数学方程

横荡运动的数学形式如下：

$\{\begin{array}{l}\xi =a\sin \left(\omega t\right)\\ \psi =\stackrel{\dot{}}{\psi }=0\\ v=\stackrel{\dot{}}{\xi }=a\omega \cos \left(\omega t\right)\\ \stackrel{\dot{}}{v}=-a{\omega }^2\sin \left(\omega t\right)\end{array}$ (21)

其中， $\xi$ 为滑翔机的横向位移，a为滑翔机横荡运动的振幅， $\omega$ 为滑翔机横荡运动的圆频率， $\psi ,\stackrel{\dot{}}{\psi }$ 为滑翔机绕z轴的倾斜角度和角速度， $v,\stackrel{\dot{}}{v}$ 为滑翔机横荡横向速度和加速度。

设在运动过程中滑翔机受到的横向力为Y，绕Z轴的力矩为艏摇力矩N，滑翔机在Y、N作用下，做横荡运动，其简化的线性运动方程为：

$\{\begin{array}{l}Y=Y_{\stackrel{\dot{}}{v}}\stackrel{\dot{}}{v}+Y_{v}v+Y_0=-a{\omega }^2{Y}_{\stackrel{\dot{}}{v}}\sin \left(\omega t\right)+a\omega Y_v\cos \left(\omega t\right)+Y_0\\ N=N_{\stackrel{\dot{}}{v}}\stackrel{\dot{}}{v}+N_{v}v+N_0=-a{\omega }^2{N}_{\stackrel{\dot{}}{v}}\sin \left(\omega t\right)+a\omega N_v\cos \left(\omega t\right)+N_0\end{array}$ (22)

将(22)式进行无量纲化，横向力方程遍除 $0.5\rho V^2{L}^2$ ，艏摇力矩方程遍除 $0.5\rho V^2{L}^3$ ，得到下式：

$\{\begin{array}{l}{Y}^{\prime }=-\frac{a{\omega }^{2}L}{V^2}{Y^{\prime }}_{\stackrel{\dot{}}{v}}\sin \left(\omega t\right)+\frac{a\omega }{V}{Y^{\prime }}_v\cos \left(\omega t\right)+{Y^{\prime }}_0\\ N^{\prime }=-\frac{a{\omega }^{2}L}{V^2}{N^{\prime }}_{\stackrel{\dot{}}{v}}\sin \left(\omega t\right)+\frac{a\omega }{V}{N^{\prime }}_v\cos \left(\omega t\right)+{N^{\prime }}_0\end{array}$ (23)

式中：L为参考长度，V为来流速度。

采用SUBOFF潜艇模型开展垂荡实验的数值模拟验证。根据试验情况，数值计算时入口处的来流为V = 4.5 kn，给定正弦振荡振幅α = 1.0 m，模拟工况振荡频率f为0.01 Hz。SUBOFF潜艇模型在均匀流场中保持正弦振荡运动，这样合成的结果即为SUBOFF潜艇模型的升沉运动。

4.2. 水动力系数求解分析

采用最小二乘法求解水动力系数。令：

$y=y\left(t\right)=Y^{\prime }-{Y^{\prime }}_0$ (24)

两个基函数分别为：

${\phi }_1\left(t\right)=v=v_{\max }\cos \left(\omega t\right)$ (25)

${\phi }_2\left(t\right)=\stackrel{\dot{}}{v}=-\omega v_{\max }\sin \left(\omega t\right)$ (26)

稳定性参数分别记为：

${\alpha }_1={Y^{\prime }}_v$ (27)

${\alpha }_2={Y^{\prime }}_{\stackrel{\dot{}}{v}}$ (28)

设 ${\phi }^{*}\left(t\right)={\alpha }_1{\phi }_1\left(t\right)+{\alpha }_2{\phi }_2\left(t\right)$ 为离散数据的拟合函数，为使该函数为最小二乘解，即满足：

${\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left[{\phi }^{*}\left(t_i\right)-y_i\right]}^2}=\min {\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left[\phi \left(t_i\right)-y_i\right]}^2}$ (29)

若令：

$S={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left[\left({\displaystyle \underset{k=1}{\overset{2}{\sum }}{\alpha }_k{\phi }_k\left(t_i\right)}\right)-y_i\right]}^2}$ (30)

则S应该满足：

$\frac{\partial S}{\partial {\alpha }_k}=0,k=1,2$ (31)

即：

$\{\begin{array}{l}{\alpha }_1{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_1\left(t_i\right){\phi }_1\left(t_i\right)}+{\alpha }_2{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_1\left(t_i\right){\phi }_2\left(t_i\right)}={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_1\left(t_i\right)y\left(t_i\right)}\\ {\alpha }_1{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_2\left(t_i\right){\phi }_1\left(t_i\right)}+{\alpha }_2{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_2\left(t_i\right){\phi }_2\left(t_i\right)}={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_2\left(t_i\right)y\left(t_i\right)}\end{array}$ (32)

作如下变量代换：

$A={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_1\left(t_i\right){\phi }_1\left(t_i\right)}$ (33)

$B={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_1\left(t_i\right){\phi }_2\left(t_i\right)}$ (34)

$C={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_2\left(t_i\right){\phi }_2\left(t_i\right)}$ (35)

$X={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_1\left(t_i\right)y\left(t_i\right)}$ (36)

$Y={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\phi }_2\left(t_i\right)y\left(t_i\right)}$ (37)

则方程组简化为：

$\{\begin{array}{l}{\alpha }_{1}A+{\alpha }_{2}B=X\\ {\alpha }_{1}B+{\alpha }_{2}C=Y\end{array}$ (38)

解之可得：

${Y^{\prime }}_v={\alpha }_1=\frac{XC-YB}{AC-B^2}$ (39)

${Y^{\prime }}_{\stackrel{\dot{}}{v}}={\alpha }_2=\frac{YA-XB}{AC-B^2}$ (40)

同理可以求得 ${N^{\prime }}_v,{N^{\prime }}_{\stackrel{\dot{}}{v}}$ 。

计算得到的水动力系数计算值与理论值对比如表3所示。横荡试验计算值与实验值差别不大。较好的验证了横荡试验数值模拟方法及水动力系数公式推导的正确性。

Table 3. Comparison of the calculated hydrodynamic coefficient of the transverse test with the test value

表3. 横荡试验水动力系数计算值与试验值对比

5. 艏摇试验数值模拟验证

5.1. 艏摇运动及数学方程

滑翔机做水平面内漂角 $\beta =0$ 的运动，称为艏摇运动。其运动方程如下：

$\{\begin{array}{l}\psi ={\psi }_0\sin \left(\omega t\right)\\ r=\stackrel{\dot{}}{\psi }={\psi }_0\omega \cos \left(\omega t\right)\\ \stackrel{\dot{}}{r}=-{\psi }_0{\omega }^2\sin \left(\omega t\right)\\ v=\stackrel{\dot{}}{v}=0\end{array}$ (41)

式中：

$\psi$ ——滑翔机运动过程中的艏向角。

${\psi }_0$ ——滑翔机的最大艏向角。

$\omega$ ——滑翔机艏摇运动的圆频率。

$r,\stackrel{\dot{}}{r}$ ——滑翔机的艏摇角速度和艏摇角加速度。

设运动过程中滑翔机受到的横向力为Y，绕z轴的力矩为纵倾力矩N。根据水动力模型，并考虑到艏摇运动时唯一的机动参数为角速度r，可得到艏摇运动线性运动方程：

$\{\begin{array}{l}Y=Y_{\stackrel{\dot{}}{r}}\stackrel{\dot{}}{r}+Y_{r}r+Y_0=-{\theta }_0{\omega }^2{Y}_{\stackrel{\dot{}}{r}}\sin \left(\omega t\right)+{\theta }_0\omega Y_r\cos \left(\omega t\right)+Y_0\\ N=N_{\stackrel{\dot{}}{r}}\stackrel{\dot{}}{r}+N_{r}r+N_0=-{\theta }_0{\omega }^2{N}_{\stackrel{\dot{}}{r}}\sin \left(\omega t\right)+{\theta }_0\omega N_r\cos \left(\omega t\right)+N_0\end{array}$ (42)

将(42)式进行无量纲化，垂向力方程遍除 $0.5\rho V^2{L}^2$ ，纵摇力矩方程遍除 $0.5\rho V^2{L}^3$ ，得到下式：

$\{\begin{array}{l}{Y}^{\prime }={Y^{\prime }}_{\stackrel{\dot{}}{r}}\stackrel{\dot{}}{r}+{Y^{\prime }}_{r}r+{Y^{\prime }}_0\\ N^{\prime }={N^{\prime }}_{\stackrel{\dot{}}{r}}\stackrel{\dot{}}{r}+{N^{\prime }}_{r}r+{N^{\prime }}_0\end{array}$ (43)

采用SUBOFF潜艇模型开展艏摇实验的数值模拟验证。根据试验情况，数值计算时入口处的来流为V = 4.5 kn，给定正弦振荡振幅 ${\theta }_0=\text{1}\text{rad}$ ，模拟工况振荡频率f为0.1 Hz。SUBOFF潜艇模型在均匀流场中保持正弦振荡运动，这样合成的结果即为SUBOFF潜艇模型的艏摇运动，求解参数设置及边界条件与SUBOFF潜艇模型直航实验保持一致。

5.2. 水动力系数求解分析

强迫滑翔机做小振幅低频率简谐振动时，将非定常水动力/力矩进行傅里叶展开，以横向力为Y为例：

$Y\left(t\right)=Y_0+\bar{Y}\sin \left(\omega t+\tau \right)+u\left(t\right)$ (44)

式中 $Y_0$ 为基准状态的横向力， $\bar{Y}$ 为基频谐波分量的水动力幅值， $\tau$ 为滑翔机强迫振动时位移和水动力之间的相位差， $u\left(t\right)$ 为高次谐波分量。

当非定常计算时间足够长时，n达到一个足够大的值，就可以忽略初始效应的影响，水动力/力矩达到一个周期性稳态值。

令 $\omega t=2n\pi +\frac{1}{2}\pi \left(n=1,2,3,\cdots \right)$ ，则有：

${\left(\frac{\partial Y}{\partial r}\right)}_0=\frac{\bar{Y}\sin \left(\tau +\frac{1}{2}\pi \right)}{r_{\max }}=\frac{Y\left[\left(4n+1\right)\pi /2\omega \right]-Y_0}{r_{max}}$ (45)

令 $\omega t=2n\pi \left(n=1,2,3,\cdots \right)$ ，则有：

${\left(\frac{\partial Y}{\partial \stackrel{\dot{}}{r}}\right)}_0=\frac{\bar{Y}\sin \tau }{-\omega r_{\max }}=\frac{Y\left(2n\pi /\omega \right)-Y_0}{-\omega r_{\max }}$ (46)

其中， $\bar{Y}$ 、 $\tau$ 、 $Y\left[\left(4n+1\right)\pi /2\omega \right]$ 、 $Y\left(2n\pi /\omega \right)$ 可以通过对非定常流场进行数值模拟得到的非定常水动力/力矩随时间的变化曲线得到。采用同样方法可以计算关于纵倾力矩N的水动力系数。

计算得到的水动力系数计算值与理论值对比如表4所示。艏摇试验计算值与实验值差别不大。较好的验证了艏摇试验数值模拟方法及水动力系数公式推导的正确性。

Table 4. Comparison of the calculated hydrodynamic coefficient of the bow shake test with the test value

表4. 艏摇试验水动力系数计算值与试验值对比