潜艇破冰上浮近场动力学模型

叶礼裕，王超，郭春雨，常欣

哈尔滨工程大学船舶工程学院，黑龙江哈尔滨150001

0 引 言

随着全球气候变暖，北极地区潜在的政治、经济、科技和军事价值日益突显，资源开采和航道开通的需求也在不断提高。目前，极地地区已成为国际政治经济竞争与合作的热点地区，吸引着许多国家的关注。随着人们对北极地理和气候环境研究的深入，越来越深刻认识到潜艇在北极地区所能发挥的政治和军事价值。北极海冰的平均厚度可达到3 m，其特殊的环境为潜艇的活动提供了特殊的隐蔽条件。基于北极地区特殊的地理位置，掌握潜艇在冰层下自由航行和活动的技能，将形成特殊的地缘战略控制能力，从而在大国之间的战略竞争中发挥特殊作用。要想胜任在北极地区的作战任务，潜艇必须具有破冰上浮的能力。由文献［1］可知，为了满足北极的发展战略，美、俄两国现阶段设计和建造的大部分核潜艇都需要满足破冰要求。但是，潜艇在破冰上浮的过程中有可能会陷入危险状态，从而使潜艇的安全受到威胁。因此，开展潜艇破冰上浮过程的科学研究至关重要。

由于潜艇冰区航行的研究涉及到一个国家的关键技术，公开资料非常少，可查阅的资料均为国家层面的军事竞争。自二战结束后，美、苏两国潜艇在冰面下的争夺愈演愈烈。1946～1952年，美国派出了多艘常规动力潜艇前往北极侦察和收集资料，均未能完成使命。常规动力潜艇因续航力有限，需要频繁浮出水面充电，故不适合在冰下长期作战，而核潜艇则具有无限续航能力，可长时间潜航，因而成为角逐北冰洋的新战略力量［2］。自核潜艇诞生以来，军事家和战略家们就开始探讨将其部署到北冰洋以发动战略攻击。1958年，美国第1艘核动力潜艇“鹦鹉螺”号成功到达北极点。自此以后，美国海军便开始强化其在北极的水下力量。2009年，美国海军在北极地区进行了特殊的军事演习，训练潜艇在北极条件下的行动能力。截至目前，美国海军的多艘核潜艇已在北极军事演练中从水下破冰而出。同样，俄罗斯海军在北极冰下也具有较强的实力，北极一直是其潜艇活动的战略性海域。1995年8月，俄罗斯海军利用“台风”级核潜艇在北极冰下进行了SS-N-20导弹发射试验。2006年，“台风”级核潜艇浮出冰面，沿着冰层裂缝破冰航行，并进行了导弹试射。2009年，俄罗斯海军“德尔塔-Ⅳ”型战略核潜艇成功躲避美国的侦察，发射了2枚弹道导弹。

一直以来，人们都比较重视极地航行船舶的科学研究，已建立了可靠的实验和数值研究手段，但潜艇在极地地区运行的安全性研究却并未受到重视，到目前为止还未建立理论计算、数值模拟和实验测量等研究手段。对潜艇破冰上浮过程模拟的关键是对海冰破碎过程的模拟。目前，在海冰破碎过程模拟方法中，应用最广的是有限元法和离散元法。但由于海冰在不同应变率下会呈现出2种不同的力学性质，故采用有限元法难以准确模拟冰的真实本构以及处理冰的极端大变形问题。离散元法适用于非连续介质力学问题的计算方法［3-6］，但在连续体阶段的结果偏差制约了其在动态破坏问题上的应用。一些学者还将光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH）方法应用于海冰动力学问题的求解。尽管如此，采用数值方法预报和分析冰—艇体的相互作用问题依然是一个比较有效的方法，必须采用更为先进的数值计算方法来解决此类问题。近场动力学方法是一种无网格方法，在计算断裂和大尺度变形问题方面具有独特的优势［7-8］，非常适合海冰破碎问题［9］。近场动力学方法与SPH方法均为粒子法，但两者的求解思路不同：近场动力学方法是借助物质点法和分子动力学思想，对连续介质力学进行分析，并以积分形式构建物体的运动方程；而SPH方法则是通过在计算域内布置粒子，将这些粒子的变量通过核函数进行插值近似来实现积分计算。

本文将采用近场动力学方法解决冰—艇体的接触问题，并模拟海冰破碎过程。首先，将艇体表面看成是刚性固体且不发生变形，采用接触检测算法识别冰—艇体的接触位置。然后，基于近场动力学方法和接触识别算法，开发潜艇破冰上浮求解程序，通过跟踪每一个时刻的破冰状态来计算海冰对潜艇壳体的动态冰载荷力。最后，以美国DARPA潜艇模型SUBOFF为研究对象，开展潜艇破冰过程的数值预报和特性研究，并对计算结果进行分析。

1 近场动力学理论

1.1 基本方程

近场动力学是一种非局部连续方法，其提供了一种更加实用的粒子方法，能从宏观角度模拟大尺度破坏和变形问题。近场动力学方法将连续介质离散为均匀的物质点。在参考坐标系下，物质点x的位置与时间t的物质点运动方程为［8］

式中：Hx为由临近物质点x的所有其他物质点所构成的域；Vx′为物质点x′的体积；u为物质点x的位移；ρ为材料密度；f为物质点x′与物质点x之间相互作用的力密度；b为物质点受到的外力。由式（1）可知，近场动力学运动方程采用的是空间积分形式，故可在任何位置处进行计算。

为表述方便，将2个物质点之间的相对位置定义为ξ=x′-x，2个物质点之间的相对位移定义为η=u(x′，t)-u(x，t)。于是，2个物质点在 t时刻的相对位置可以表示为ξ+η。力密度可以表示为f(η，ξ)，其大小依赖于ξ和η。2个物质点之间的相互作用可以称为“键”，类似于一对相互作用的弹簧力，如图1所示。近场动力学定义了一个近场范围尺寸δ，物质点x只与距离δ半径的球形范围内的物质点x′相互作用，即

对于一个连续体，成对的力可以视为对物质点之间的相互作用，并且满足线动量和角动量的关系：

1.2 微观弹脆性（PMB）模型

对于由微弹性材料构成的物体，力密度函数可以表示为

式（1）和式（5）共同构成了键型近场动力学材料的本构模型。为了描述物质点之间键的伸长程度，近场动力学方法使用了非线性弹性材料的键伸长率：

当键处于拉伸状态时，s为正值。由于s不依赖于ξ的方向，因此这种材料是各向同性的。为了在本构模型中引入材料破坏的概念，使用了一个比较简单的判定条件，即当某一时刻的键伸长率超过一定值时，认为物质点之间的键永久发生断裂。从此刻起，此对物质点之间的力密度为0。该条件可称为近场动力学的破坏准则。因此，对于典型的微观弹脆性（Prototype Micro-elastic Brittle，PMB）材料，力密度函数可以简化为

式中，g为线性标量函数，可定义为

式中，c为均匀各向同性材料的标量微模量常数，可定义为

式中，κ为体积模量。

式（7）中的μ为历史变形判断标准，当键伸长率达到破坏准则时，键发生断裂，μ=0，反之，μ=1，可以定义为

式中，s0为键破坏时键伸长率的极限值，可称为极限伸长率。对于典型的PMB材料，s0可由下式计算得到：

式中，G0为裂缝扩展时的能量释放率。尽管弹脆性材料在初始条件下是各向同性的，但是有些键断裂之后会引起材料的各向异性。为此，引入键的破坏水平参数

式中，Vξ为物质点ξ的体积。

1.3 方程离散化

为了便于开展数值计算，需要将连续体材料离散成物质点，然后对每个物质点进行积分运算［9］。离散后，每个物质点x的运动方程可以表示为

式中：m为时间步；下标i和j为物质点的编号，其中i为要计算的物质点，j为临近x的物质点；Vj为物质点j的体积。

通过式（13），可以求得每个物质点的加速度。为了得到每个物质点的位移，可以通过中央差分法计算得到：

式中，Δt为时间步长。

2 潜艇破冰上浮计算方法

2.1 问题描述

北极大部分地区被海冰覆盖，极大地限制了水面舰船的航行，但潜艇却可以潜入冰层之下自由航行，不仅能躲避反潜巡逻机和反潜水面舰船的侦察与攻击，而且极地地区冰层的膨胀、破裂和因摩擦产生的环境噪声还能掩盖潜艇本身的噪声。考虑到北极地区作战任务的需要，要求潜艇不仅要能在冰下长久航行，关键时刻还需要能够破冰上浮，发射弹道导弹或者对外通信联络或者接受救援补给。

图2所示为潜艇破冰上浮情况。在潜艇破冰上浮的过程中，主要受力部位是指挥室围壳和艇的上层建筑，通常需要对这些部位进行结构优化设计，以满足强度要求。极地地区的冰面通常无限宽广，潜艇在上浮过程中只能使壳体周围的海冰破碎。

参考潜艇破冰上浮的实际情况，建立如图3所示的数值模型。该模型假定潜艇以一定的速度上浮，并与上方静止的冰面接触。为便于数值计算，将平整冰简化为长方体，而且长度和宽度远大于厚度。由于实际的海冰尺寸相对于艇体要大得多，因此在进行模拟时冰块的长度和宽度要足够大，两边加固定边界条件以代替无限大域的边界问题。另外，考虑到潜艇在破冰上浮过程中流体对冰载荷的影响较小，忽略了流体的干扰作用。本文的主要目的是计算冰对艇体的接触力并模拟海冰的破碎过程，故将潜艇壳体看成是刚性固体，暂不考虑壳体变形的影响。

为了考虑海冰粒子重力和浮力的影响，根据浮力和重力计算公式，当物质点位于冰排上方时,式（1）中外力的取值为b(x，t)=ρVg；当物质点位于冰排下时，外力的取值为b(x，t)=0.1ρVg，其中V为物质点的体积。

2.2 模型数值离散化

本文的海冰模型是采用近场动力学方法进行模拟，在数值计算过程中，需要将冰块连续体离散为物质点形式。但是，需要特别注意物质点间距Δx和邻域半径δ的选择，它们是影响计算精度的重要参数，需要综合考虑计算精度和速度，以选取最优值。当邻域半径过小时，裂缝扩展过程对网格有依赖性，难以模拟材料的大尺度变形，因此建议采用δ=3Δx作为邻域半径［10］。

对于艇体模型，本文将艇体简化为刚性表面，并用四边形面元来划分。本文以美国的SUBOFF潜艇模型为研究对象，该模型的主要参数见第3.1节。将艇体三维模型离散为一系列四边形面元形式，面元越小，越能逼近艇体的真实形状。图4所示为艇体表面网格划分结果。

2.3 海冰物质点与潜艇壳体的接触识别

当海冰与潜艇壳体表面离散化后，可将海冰与壳体表面的接触判断转换为物质点和四边形面元的相对位置判断，从而将复杂的接触过程简化为一系列点和面的数学问题，有效降低接触检测计算的复杂度。

为减少计算量和运算的复杂度，同时保证接触检测的精度和效率，基于接触检测中包围盒的思想［11］，本文给出了一种简单高效的方法来剔除不可能与艇体表面发生接触的物质点。建立一个长方体将整个艇体包围在内，该长方体的长度等于艇长Ls，宽度等于艇宽Bs，高度等于艇体最高点到最低点z方向（垂向）的距离。海冰的物质点位于长方体内部被认为可能与艇体接触的物质点。

对于每一个可能与艇体接触的物质点，若与艇体发生接触，只能有一个四边形面元与物质点接触。为了提高搜索效率，通过物质点和四边形4个点的坐标位置关系来确定可能与该物质点接触的面元。假设粒子的坐标为(x0，y0，z0)，对于艇体表面的所有四边形面元，可以找到其4个点在x轴方向上最小值xmin和最大值xmax，在y轴方向上的最小值ymin和最大值ymax，以及在z轴方向上的最小值zmin和最大值zmax。若面元与点的关系 为 ：xmin＜x0＜xmax且zmin＜z0＜zmax；xmin＜x0＜xmax且ymin＜y0＜ymax，则认为这些面元可能与该物质点发生碰撞。

对于每一个可能与艇体接触的物质点，可计算得到可能与其发生碰撞面元的控制点之间的距离，从而找到距离最短的面元。通过点与面的数学公式，判断物质点与其距离最短面元的位置关系，如式（15）所示。

式中，A1，B1，C1，D1为四边形面元所在空间平面方程的系数。

于是，与艇体表面发生接触的物质点和与其发生接触的面元即可确定。

2.4 接触作用力计算

潜艇在上浮的过程中，海冰物质点将与潜艇壳体发生接触。接触以后，物质点就会渗入到艇体内，如图5（b）所示。为了反映真实的物理情况，渗入到艇体内的物质点需要进行位置再分配，将该物质点分配到临近处的艇体表面，如图5（c）所示。通过下式，可以计算出渗入到艇体内的物质点和新分配物质点的距离为

新分配物质点的坐标可以通过下式计算：

式中：下标k为物质点标号；v0为物质点k的速度；n为接触点处的法向量。

物质点x(k)在t+Δt时刻的速度可用式（18）计算：

3 潜艇破冰上浮算例分析

3.1 计算模型及参数设置

海冰破碎方式的多样性以及不同海域海冰的物理和力学性质有很大的不同，导致海冰的物理和力学性质不太容易确定，较好的方式是通过采样相关海域的海冰来试验和测量其物理性质，例如压缩强度、弹性模量和密度等。由于本文关注的重点不在于海冰的物理性质，而是对海冰破坏现象的模拟，因此将海冰简化为匀质材料，而海冰的物理和力学参数与文献［12］相同，即弹性模量E=1.8 GPa，泊松比υ=0.25，临界伸长率s0=0.02，密度ρ=900 kg/m3。

本文的研究对象选用的是美国DARPA SUBOFF实尺度潜艇模型，如图6所示。模型的型值取自文献［13］，全附体潜艇长104.544 m，艇身最大直径12.192 m，指挥室围壳长8.832 m，平行中体为53.496 m，尾翼的翼型后缘位于96.168 m处，4个尾翼剖面均为NACA 0020翼型，并呈对称布置［14-15］。

3.2 计算模型的收敛性分析

在开展潜艇破冰上浮特性研究之前，首先针对物质点间隔Δx和时间步长Δt对计算结果的影响进行分析。将冰层尺度定为长度L=150 m、宽度B=36 m、厚度T=0.9 m的长方体，分析潜艇以0.5 m/s速度破冰上浮过程中冰载荷的动态响应特性。由于潜艇在破冰上浮过程中z轴方向（垂向）的力远大于x轴方向（纵向）和y轴方向（横向）的力，故本文重点分析z轴方向的冰载荷。

为了分析物质点间距对计算结果的影响，在其他工况相同的情况下，将冰层离散成间距分别为Δx=L/250，L/500和L/750的物质点，时间步长控制为Δt=0.006 25 s。图7给出了冰载荷Fz随时间变化的曲线。随着物质点间距的减小，计算曲线逐渐收敛。物质点间距Δx=L/250的冰载荷大小要比间距Δx=L/500和L/750的稍小，曲线振荡也比较剧烈。由于冰和结构物相互作用有很大的随机性，冰载荷将发生剧烈的动态变化，因此相同工况下的冰载荷时程曲线无法完全重合。由图可知，物质点间距Δx=L/500和L/750的计算曲线基本一致。因此可以认为，当物质点间距小于Δx=L/500时，计算结果已能够达到收敛。

同理，控制物质点间距为Δx=L/500不变，时间步长分别取为 Δt=0.000 948，0.000 625和0.000 312 s。图8所示为不同时间步长时冰载荷Fz随时间变化的曲线。由图8可知，尽管3条曲线在某些时间点上有一些振荡，但总体来看基本还是一致的。因此可以认为当时间步长小于0.000 948 s时，计算结果可达到收敛。

综上所述，本文计算模型的收敛性可以得到保证。在后续的计算中，兼顾计算速度和精度，物质点间距取为 Δx=L/500，时间步长取为 Δt=0.000 625 s。

3.3 潜艇破冰上浮动态特性分析与验证

本节以潜艇以0.5 m/s速度连续破冰上浮撞破0.9 m厚的冰层为例，来分析潜艇破冰上浮过程中海冰的破碎特点和冰载荷的动态响应特性。选择的冰层是长度L=150 m、宽度B=36 m的长方体。将离散后的海冰物质点之间的间隔取为Δx=L/500，将时间步长设定为Δt=0.000 625 s。

图9给出了不同时刻潜艇破冰上浮时海冰的挤压和弯曲破坏过程。图9中云图的数值表示冰物质点的破碎水平，颜色的深度表示冰粒子的损伤度，其中蓝色表示冰粒子处于无损状态，红色代表冰粒子处于完全损伤状态。由图9可知，潜艇指挥室围壳上部首先撞击到冰面，由于它垂直于冰面，围绕着指挥室围壳的海冰破碎，指挥室围壳首先破冰而出，并且其上部还出现了一块独立于冰面的海冰，伴随着指挥室围壳一起向上运动。随后，潜艇尾翼与海冰接触，冲破冰层露出冰面。随着破冰过程的进行，潜艇上层建筑也与海冰接触，冰面出现了大面积的挤压和弯曲破坏。

下面，将通过与网上可找到的美、俄核潜艇在北极冰下破冰而出的画面进行对比来验证本文计算模型的有效性。图10所示为美国核潜艇在北极地区破冰而出的最终画面，将其与图9（f）进行比较后发现，潜艇上层建筑的海冰分布基本一致。图9（f）给出了海冰物质点损伤度的分布，可以看出，潜艇上层建筑的海冰受到潜艇的冲击作用发生了破碎。潜艇破冰上浮后，只有指挥室围壳和尾翼露出冰面，而潜艇上层建筑则被掩盖在海冰之下。

图11所示为本文计算的潜艇指挥室围壳和美国核潜艇指挥室围壳破冰而出的局部放大图。通过比较可知，两者的海冰破坏形式基本一致，指挥室围壳上部均出现了一块独立于冰面的海冰。由图11（a）可知，在当前计算工况下，指挥室围壳冲击海冰的过程中，其前端的海冰发生了纵向剪切破坏，而后端的海冰则以弯曲破坏为主。可以推测，在指挥室围壳冲击海冰的过程中，海冰的破坏模式与指挥室围壳形状有关，当然也可能与潜艇的上浮速度和海冰的厚度等有关。

图12所示为本文计算的潜艇尾翼和美国核潜艇尾翼破冰而出的局部放大图。通过比较可知，在当前计算工况下，两者尾翼在冲击海冰的过程中均以弯曲破坏为主，破坏形式基本一致，这种破坏模式不同于指挥室围壳冲击海冰的过程。这主要是由于与指挥室围壳相比，尾翼更加窄小，降低了与海冰接触的面积，因而受力面积更小，使得海冰更加易于破碎。

图13所示为当前计算工况下潜艇上浮过程中冰载荷Fz的时程曲线。在刚开始阶段，由于指挥室围壳对海冰的冲击作用，潜艇受到的冰载荷出现小幅度的波动。3.0 s以后，潜艇受到的冰载荷进入一个比较平稳的过程，大小约为168.3 kN。在该过程中，海冰受到指挥室围壳的冲击作用结束，指挥室围壳破冰而出，此时潜艇主要受指挥室围壳上部一块独立于冰面的海冰的作用，因此该过程中潜艇受到的冰载荷约等于这块海冰的重力。通过测量，得到指挥室围壳顶部面积S约为28.1 m2，而海冰的厚度T=0.9 m，从而可以大概估计出这块海冰的重力为228.2 kN，稍大于该过程所受到的冰载荷，从而验证了本文接触力计算方法的有效性。在8.5～9.0 s之间，冰载荷在波动中上升，该时间段海冰受到尾翼和潜艇上壳体的共同作用，因此冰载荷要大于指挥室围壳与海冰作用的峰值。在9.0～10.5 s之间，冰载荷迅速增大，该时刻曲线比较平滑，可以推测在该过程中潜艇与海冰主要发生了挤压作用。在10.5 s以后，冰载荷开始出现剧烈的振荡，这可能是由于在该过程中海冰不断发生挤压和弯曲破坏所引起，这将会导致艇体的激烈振动，有可能会对艇体结构安全性构成威胁。当潜艇上壳体上升到一定程度以后，冰载荷开始迅速减小，表明破冰过程将要结束。

为了更好地分析潜艇指挥室围壳与海冰相互作用的冰载荷特性，单独绘制了该过程冰载荷Fz的时程曲线，如图14所示。指挥室围壳上部与海冰发生接触后，海冰在未发生破碎之前主要以挤压破坏为主，冰载荷迅速增大，此时曲线比较平滑。随后，海冰发生弯曲或者剪切破坏，冰载荷曲线发生剧烈振荡。当潜艇指挥室围壳破冰而出时，此时海冰受到指挥室围壳的冲击作用结束，冰载荷迅速减小。

4 结 论

本文基于近场动力学方法，建立了潜艇破冰上浮数值计算模型，考虑了海冰浮力和重力的影响，并将其加入到了近场动力学基本方程中。通过运用计算程序对潜艇的破冰上浮过程和冰载荷动态响应进行计算，分析得到了以下结论：

1）开展了潜艇破冰上浮过程中物质点间隔和时间步长的收敛性分析，计算结果表明，物质点间隔小于Δx=L/500、时间步长小于0.000 948 s时，计算结果可达到收敛。

2）基于本文建立的数值模型对潜艇破冰上浮过程中海冰的破坏过程进行模拟，很好地再现了海冰的挤压、纵向剪切破坏或弯曲破坏过程，证明了本文计算方法的有效性以及接触识别算法的可靠性。

3）对潜艇破冰上浮过程冰载荷时程曲线的分析表明，在海冰挤压破坏过程中，冰载荷迅速增大，曲线比较平滑，但海冰弯曲破坏过程一旦发生，冰载荷将出现剧烈的振荡，对艇体结构的安全性不利。

本文建立的潜艇破冰上浮过程计算模型，和对潜艇破冰上浮的动态特性进行的计算分析，可以用于指导潜艇壳体结构的优化设计。本文并未考虑流体与艇体和海冰之间相互作用的影响，计算结果可能存在误差。下一步，将更全面地考虑各种因素的影响，建立更加可靠的数值预报方法。

参考文献：

［1］刘征鲁，李春迪.核潜艇破冰要有哪些“硬功夫”？［N］.中国国防报，2016-05-06（014）.

［2］陆俊元.北极地缘政治与中国应对［M］.北京：时事出版社，2010.LU J Y.Geopolitics in the Arctic and China's response［M］.Beijing：Current Affairs Press，2010（in Chi⁃nese）.

［3］WILLIAMS J R，HOCKING G，MUSTOE G G W.The theoretical basis of the discrete element method［C］//Proceeding of the International Conference on Numeri⁃cal Methods in Engineering：Theory and Applications.Rotterdam，Netherlands，1985.

［4］LAU M，LAWRENCE K P，ROTHENBURG L.Dis⁃crete element analysis of ice loads on ships and struc⁃tures［J］.Ships and Offshore Structures，2011，6（3）：211-221.

［5］HAN K，PERIC D，CROOK A J L，et al.A combined fi⁃nite/discrete element simulation of shot peening pro⁃cesses-Part I：studies on 2D interaction laws［J］.Engi⁃neering Computations，2000，17（5）：593-620.

［6］HASHASH Y M A.Large-scale numerical simulations via parallel computing：an application using discrete element modeling of granular material［C］//NCSA Stra⁃tegic Applications Program Status Report 03/28/2006.Champaign，IL：Center for Super Computing Applica⁃tions，University of Illinois，2006.

［7］SILLING S A.Reformulation of elasticity theory for dis⁃continuities and long-range forces［J］.Journal of the Mechanics and Physics of Solids，2000，48（1）：175-209.

［8］SILLING S A，ASKARI E.A meshfree method based on the peridynamic model of solid mechanics［J］.Com ⁃puters&Structures，2005，83（17/18）：1526-1535.

［9］YE L Y，WANG C，CHANG X，et al.Propeller-ice con⁃tact modeling with peridynamics［J］.Ocean Engineer⁃ing，2017，139：54-64.

［10］MADENCI E，OTERKUS E.Peridynamic theory and its applications［M］.New York：Springer，2014.

［11］NEZAMI E G，HASHASH Y M A，ZHAO D W，et al.A fast contact detection algorithm for 3-D discrete el⁃ement method［J］. Computers and Geotechnics，2004，31（7）：575-587.

［12］ZHAO G L，XUE Y Z，LIU R W，et al.Numerical sim⁃ulation of ice load for icebreaker based on peridynam⁃ic［C］//Proceedings of the ASME 35th International Conference on Ocean，Offshore and Arctic Engineer⁃ing.Busan，South Korea，2006.

［13］GROVES N C，HUANG T T，CHANG M S.Geometric characteristics of DARPA SUBOFF models（DTRC ModelNos.5470 and 5471）：reportDTRC/SHD-1298-01［R］.Bethesda，USA：David Taylor Research Center，1989.

［14］操盛文，吴方良.尺度效应对全附体潜艇阻力数值计算结果的影响［J］. 中国舰船研究，2009，4（1）：33-37，42.CAO S W，WU F L.Investigation of scaling effects on numerical computation of submarine resistance［J］.Chinese Journal of Ship Research，2009，4（1）：33-37，42（in Chinese）.

［15］吴方良，吴晓光，马运义，等.潜艇实艇阻力预报方法研究［J］. 中国舰船研究，2009，4（3）：28-32.WU F L，WU X G，MA Y Y，et al.Method of predict⁃ing the total submerged resistance of submarines［J］.Chinese Journal of Ship Research，2009，4（3）：28-32（in Chinese）.