莫里森Morison公式和结构动力学的再学习
一、背景
许多同学进行海洋工程结构物分析，都是从应用“莫里森公式”开始的。相比于《流体力学》中各种复杂的数学公式和推导，莫里森公式绝对是“独一挡”的存在，公式简明、工程适用性广。从导管架圆管构件，到自升式平台的圆柱腿，甚至再大一点到浮式风机三立柱，按照书上的说法，只要特征长度（如直径）不大于1/5的波长，基本都适用。笔者结合这些年来的一些工作经验，分享、讨论一些对莫里森公式及结构动力学的一点看法和理解。
二、问题
在应用莫里森公式这种半经验方法时，比较容易“一知半解”或者“稀里糊涂”，因为“规范”上有了比较明确的“建议”。翻看在课本教材上，莫里森公式都出现在“小尺度结构物”的章节，与“大尺度结构物”分的很开（或者说毫不相关）。这就好像如今的专业工作划分，做固定式结构的和做浮体结构的是两拨“最熟悉的陌生人”，“融会贯通”四个字在这里存在困难。
另一方面，莫里森公式虽然很便捷，但其中的“平方”运算（非线性项）却“惹”出了很多麻烦，“线性化”的概念理解起来总有点“玄乎”的感觉。
再有，对莫里森公式的理解，很多同学还止于“准静态”的概念。我们知道波浪载荷具有周期性特征，正确理解莫里森公式的“周期性”和结构的动力响应特征是真正用好莫里森公式的关键。
针对以上三点，笔者将展开一定的讨论，换个与书本不同的思路浅显地“再学习”一下，不推导公式、没有严格的数学过程和理论分析，面向的是结构方向的同学，难免有不严谨的地方还请包涵，参考就好。
三、莫里森公式中的水动力系数的取值
一般按书本的顺序，先给出莫里森公式，定义拖曳系数Cd、质量系数Cm两个系数，再讨论Cd、Cm与哪些因素有关，如粗糙度、 雷诺数Re、KC数等，最后讨论一下流动分离和横向力，这样一个章节就结束了。实际工程中，因为现在的规范对Cd，Cm的取值已经有明确的推荐，给工程师可选择的余地是不大的，但笔者认为，同学们仍需要明确这样几点，才能正确理解这个方法。
1） 莫里森公式的应用对象主要是波浪
莫里森公式的应用主要是针对周期性变化的“波浪”（规则波和不规则波），即流场中有“速度”和“加速度”分布，从而载荷中有相应的“拖曳”项和“惯性”项。在参数讨论中应用的Re和KC都是水质点运动参数的函数，即Re = umD/ν，KC = umT/D （um为水质点的速度幅值，T为波浪周期，ν为流体的运动粘度）。
书本中为了说明Cd（拖曳力项）与粗糙度和Re的关系，引用的是“定常流”下的Cds图。注意因为是定常流，流场水质点速度 = 流速 = 常数，不考虑杆件对流场的影响，这时流速的大小决定了Re的量级。因为ν是10-6量级（单位m2/s），假设海工中常见的流速量级u=1m/s，则无量纲的Re在量级106以上（属于高雷诺数），对于常规粗糙度k/D<0.01，再定常流下Cds稳定在1.0附近（Cds = Cd under Steady flow）。
2） 莫里森公式中的Cd，Cm是“平均值”的概念，
莫里森公式中Cd，Cm可以理解为一个波浪周期内的”平均”，而不是严格的“常数”（工程上做了简化处理，认为Cd，Cm在一个波浪下近似常数）。
3） Cd，Cm是Re，KC，粗糙度k/D的函数
理论上，莫里森公式中Cd，Cm与波浪参数相关，是Re，KC和粗糙度的函数。KC数反应了波浪的“震荡”特性。规范DNV-RP-205中，Cd由定常流下的Cds乘以一个震荡修正系数Ψ而得：
Cd = Cds • Ψ(KC)
设Cds = 1, KC = umT/D = πH/D，其中H为波高。对于小尺度结构物，典型的H/D~10代入得KC~30。可以看到典型的Ψ的取值范围为1.1~1.2。当KC比较大时，Ψ->1.0。故Cd 的取值也在1.1附近。
对于圆柱，Cm的取值主要也与KC数有关，如下图。在KC≈30，Cm≈1.2~1.6。一般保守地取Cm=1.8。
下图中给出了前辈们得到的光滑圆柱杆的Cd，Cm与Re，KC的关系图，可知在海洋工程关心的Re范围内（红色框范围），圆柱的Cd和Cm都是比较“稳定”的。Cd=0.6~0.7，Cm=1.6~1.8。与规范的建议值接近。
四、莫里森公式与达朗贝尔佯谬
相信所有学习《流体力学》的同学都通过达朗贝尔佯谬知道了势流理论的“局限性”。笔者个人认为在理解这部分内容时，应明确两点：
1） 达贝尔佯谬，即在流场中应用势流理论，固定圆柱受力为0的结论是对于“定常流”，“定常流”，“定常流”而言的（重要的事情说三遍），且与圆柱直径D的大小无关（已经考虑了绕射）。
2） 在规则波下的流场中，应用势流理论，考虑绕射，固定圆柱受力是不为0的。可将这个载荷写成与莫里森载荷惯性项相同的量纲形式，得到：
Cm = 4 • A(ka) / π • (ka)2
其中，a=D/2为圆柱的半径，k为波数（k=2π/L，L为波长），A为ka的函数，由两个特殊函数的导数确定（第一类Bessel函数和第一类Hankel函数）。可以固定k，让圆柱的半径a变化（D=2a），可得上式中Cm随D/L变化的关系图如下图所示。
可见，当D/L<0.2时，即认为绕射对流场影响可忽略，莫里森公式适用时，Cm的值的取值范围为1.8~2.1，变化不大。而D/L趋于0时Cm趋于2.0，这与莫里森方程中（光滑圆柱）的惯性项相同！
通过上面的讨论，可以把莫里森公式和势流理论联合起来考虑，而不是“非此即彼”，即：
- 莫里森公式中的惯性项是“等价于”势流理论的，尤其在D/L<0.2即所谓的“小尺度杆件”适用范围内是对势流理论的很好的“简化”；
- 莫里森公式中的拖曳项正是从工程的角度（而不是CFD的角度）解决势流理论中无法计算的，诸如粘性、流动分离等问题。Morrison公式应配合势流理论一起运用；
- 书本和规范中，给出了各种截面形状在大Re数下的Cds的推荐。在实际工程中，可以根据目标的KC值来确定修正系数得到Cd；
- Cm=1+Ca，本质上是附加质量系数Ca的求解问题。而附加质量是所谓“大尺度物体”波浪力计算的内容，是波浪周期的函数。但在初步评估那些“较大的结构物”时，仍可以根据D/L比值，“合理”地选择Cm（常数）应用莫里森公式进行计算（可先不考虑绕射），这样得到的总载荷误差总体上是可控的。
五、莫里森公式载荷作用下的结构动力响应
莫里森公式中的两项：拖曳力项和惯性力项都是随水质点速度周期变化的，由基础流体力学知识可知，当流体质点水平速度达到最大时（拖曳力最大），其水平加速度为零（惯性力为零），反之当流体质点水平加速度达到最大时（惯性力力最大），其水平速度为零（拖曳力为零）。两者不会同时达到最大，但可以先直观比较下两者在水线附近的最大值，参考DNV-RP-C205的推荐公式：
Fd,max / Fi,max = Cd/Cm • KC/π2
对于“小尺度杆件”，取KC=30，代入Cd=0.65, Cm=2.0，可得Fd,max / Fi,max ≈1.0，即拖曳力与惯性力是相当的，不可忽略；
对于“大尺度物体”，D/L>0.2，取平均波陡1/15，则H=1/15*L，即H/D<1/3，取H/D=0.3，KC=0.3π≈1，代入Cd=0.65，Cm=2.0，Fd,max / Fi,max ≈0.03 < 5%，惯性力占绝对主导。
对单个杆件来说，总的莫里森载荷达到最大，有三种可能（参考书本上具体的数学推导）：
- 惯性力达到最大时；
- 拖曳力达到最大时；
- 惯性力和拖曳力之和达到最大。
但是，因为一般的海工结构物有多个杆件，如自升式平台至少有3个桩腿，它们相互之间相互作用产生所谓“群柱”效应。下图反应了平台在一个全波浪相位下的受波浪力的历程。总的波浪力（莫里森载荷）是前排腿和后排腿的叠加，相位差=2π•Δ/L（式中Δ为前后排腿间距，L为波长）。
也就是说，对平台整体而言，总的莫里森载荷是所有桩腿受力的叠加。拖曳力项和惯性力项的相互关系也比单个杆件要复杂。主要考虑：
- 桩腿的间距，注意不同浪向下桩腿之间间距是不同的
- 波浪的波长（周期）和间距的关系（相对相位）
可以看到，由于桩腿相对间距的存在，波浪在各桩腿上的受力可能相互“叠加”，也可能相互“抵消”，即所谓的reinforcement和cancellation效应。在做结构强度分析时，应注意避免cancellation效应，否则将大大地低估所受载荷。限于文章的篇幅，这里不做过多展开，具体可查看ISO-19905有关的部分。
回到本节的重点，关于莫里森载荷下的结构动力响应。莫里森载荷的周期性来自水质点运动的周期性，而平台的“群柱”效应会对这种周期性做一定的“修正”。下图是一个典型的自升式平台的莫里森载荷拆解，总载荷Fh = Fd（拖曳力）+Fi（惯性力）。可以看到波动部分的主要来源是惯性力成分，拖曳力成分虽然大小和惯性力相当，但波动幅值不如惯性力。故莫里森载荷整体可以近似认为是以波浪周期Twave变化的“正弦激励”（但均值≠0）。
讨论所谓的结构动力响应，首先就是“共振”判别，也就是结构固有周期Tn和激励周期Twave的关系，将平台整体看作是SDOF系统，相信下图中左边所示的“幅值放大系数”即DAF图大家都很熟悉了，右边这个相位（差）图往往容易被忽视。事实上，这两幅图可以说是同等重要的，左边表示了动力响应幅值的大小，右边表示了动力响应的相位（即动力响应达到最大时和输入是最大时的时间差or相位差），这里做一梳理。
- 一般认为Tn是系统的固有属性，当Tn/Twave较小时（<0.5，即激励Twave较长），系统受力是接近“准静态”模式，这时响应和输入几乎是“同时”的，输入最大时响应几乎同时达到最大；
- 当Twave减小，即Tn/Twave从0.5增大到接近0.9时，响应达到最大时和输入最大时开始有明显的“时间差”（响应最大值时刻滞后于激励最大值时刻），当共振时（Tn/Twave=1），这个时间差达到1/4激励周期Twave（即360度/4=90度的相位差）。注意到，对于一个正弦激励，从最大激励到激励减少到0，正好相差90度相位。
- 若Twave继续减少到小于Tn，最后响应和输入将相位相反（相差半个激励周期Twave）。对于一个正弦输入，从最大输入到最小输入正好相差180度相位。
在ISO-19905中很好地阐述了这个现象，下面给出了一个4腿自升式平台在不同的Tn/Twave下，结构总动力响应和波浪载荷的幅值和相位关系的例子，与规范符合程度较高。
六、莫里森公式的“线性化”
莫里森公式中的拖曳力项与流体质点的速度的平方项成正比，也就是莫里森载荷与波浪的波高不是线性关系。在海工结构分析中，常用谱分析这种成熟的线性化方法来求解疲劳问题。这时，就需要对莫里森公式基于一定的“等效”原则，进行“线性化”处理。
在最新2023年版的DNV-RP-C104对其进行了讨论和说明（指向了ISO-19902中的一些推荐做法）。诚然，对莫里森公式进行线性化看似有些“刻意”（为了要用谱分析这个工具，强行把问题对象套进“线性化”系统，有些不严谨），但是前辈们努力寻找的“等效原则”而提出的工程上可操作的方法，还是值得大家学习和借鉴，甚至提出适合自己项目的方法。这里对莫里森公式拖曳项的线性化处理做一个简单的、大白话式的梳理：
- 方法1用一个常数C（Borgaman, 1958）
用一个经验的常数C来取代莫里森公式中的平方项v•v，变成C•v。
等效性原则：设速度v(t)为一个平稳随机过程，符合正太分布p(v)，v(t)方差为σv2。让(|v|•v-C•v)2• p(v)这个量的积分最小，可计算得到C=(8/π)0.5•σv，而σv可由速度谱Sv通过积分求得，而Sv可由波浪谱Sω求得（详见书本教材）。当然，对于大型结构物或者拖曳力项占比很小时，也可以简单取一个较小的C，如C=1m/s来进行线性化而不至于引起较大误差。
- 方法2 用一个波陡S（ISO-19902）
区别于浮式结构用单位波幅波浪计算各波浪周期下的响应，方法2另辟蹊径（或者说向后转移矛盾），用一个经验的（或者是平均的）波陡S=H/L（ISO-19902推荐1/15~1/25），来确定不同波浪周期下的波高，进而计算出响应，再把得到响应除以对应波幅来得到传递函数RAO。这种方法避免了对莫里森公式的拖曳力项直接进行线性化，而是用了“归一”法得到RAO。有了RAO，谱分析自然就可以进行下去了。但是这里的等效性原则不强，仅仅是这个S是“平均的”或者有“代表性”的，主观较强且缺乏验证。
- 方法3 选一个波陡S（ISO-19902）
方法2中，RAO的结果（进而谱分析的结果）取决于S的选取。矛盾虽然转移了，但问题没有解决，如何选一个合理的S？一个保守的方法就是选取一系列S，取一个保守的结果。但“保守”≠符合等效性原则。
ISO-19902针对疲劳分析，建议先从波浪散步图中，选择一个对疲劳最大贡献的seastate，认为其Damage 正比于 N•Hs2m/Tz，其中N为seastate在波浪散步图中出现的次数，Hs和Tz分别式seastate的有义波高和平均过零周期，m为疲劳SN曲线的斜率（一般取m=3）。其等效性原则是：选取的S（从而得到的RAO）应能使分析对象在该seastate下由谱分析方法得到的最大响应与采用设计波方法（不采用线性化处理）得到的一致。
- 方法4 选一个更好的波陡S（Simon）
我们从方法3上看到其对方法2的改进，及其更强的等效性原则。通过一些实际的计算分析总结，笔者个人给出进一步改进的方法4。改进点：在方法3选择seastate时，仅仅考虑了波浪自身的因素，而没有考虑结构特性（固有周期Tn）的因素。因此，在选择seastate时，在其Damage 中引入DAF的因素，即：
Damage正比于 N•(DAF•Hs)2m/Tz，其中DAF根据Tn和谱峰周期Tp的关系确定。
七、小结
本文从对莫里森公式的水动力系数的讨论，到莫里森载荷作用下的结构动力响应，再到莫里森公式的线性化，对莫里森公式进行了再学习和探讨，希望对广大同仁有所借鉴。
更多讨论欢迎关注笔者在仿真秀的课程：
“海洋平台强度分析 - 工程师进阶课程”，同时也欢迎大家和笔者讨论规范和工程问题。
