通常我们希望在人工地震区域上释放强度均匀的地震波, 因为变化剧烈的地震波会影响判断的效果.因此,我们希望各电缆落到海底的 位置最好呈均匀分布状态(如 图 1(a) 所示). 然而,由于海水的冲击作用, 电缆在海底的落点往往会分布不均匀(如 图 1(b)所示).因此,研究海底电缆铺设 过程中电缆线的运动轨迹,从而为探测船设计合理的施工路线和航行速度,以使得电缆 在海底的落点位置接近均匀分布状态,提高探测的准确性有着很重要的实际意义.
自Zajac[9]在这个领域的开创性工作以来,许多学者作出了相应的研究. Lenord和 Karnoski[10]提出了探测船在匀速运动施工状态下电缆运动的计算算法. Huang 和 Vassalos[11]提出了基于弹簧-质点方法的数值算法来预测施工过程中电缆的断点位置. Pate和 Vaz[3, 6]提出了一个计算电缆短时间行为的数值模型,后来他们又将此模型推 广到研究施工过程中电缆的弹性延展行为.Prpić-Oršić 和 Nabergoj[5]研究了电缆与船的交互作用. Xu[12]也提出了一个海底电缆运动数值模型.
上述大多数文献都集中在分析电缆施工过程中,缆绳张力场(实际上是弹性力学上的一 种反抗形变的应力)的分布状态、探测船的运动 状态以及海水深度对缆绳姿态的影响.除此之外,这些工作都忽略了流体升力的作用, 而且均是在移动坐标系的基础上研究电缆的局部($a(t)\leq s \leq b(t)$)行为. 为了全局$(0\leq s \leq b(t))$地分析电缆的行为,我们曾建立了一个基于固定坐标系 (GPS大地坐标系)的电缆运动模型[1],在此模型中,我们忽略了流体升力的作用 (受力分析如图 2(a)所示).现在为了更准确地分析电缆的形态,本文将额外考虑流体升力$L$ 的作用(受力分析 如图 2(b)所示).
本文的主要目的是对文献[1]中模型进行有效的改进,进而对改进后的模型方程的全局和局部 光滑的静态解的存在条件加以证明并通过数值模拟的方法来研究电缆的运动姿态. 对电缆铺设过程中设计合理的探测船航行速度和路线以保证电缆沉放到设定位置的研究提 供理论上的指导.
文章内容安排如下: 首先,介绍课题的研究背景和意义.在第二节中,我们将在固定坐标系 (大地坐标系)下建立描述电缆运动过程的微分方程模型. 在第三节中,我们讨论全局和局部光滑静态解的存在条件.在第四节中,我们提出求解微分 方程模型的数值计算算法,并通过数值试验的方法来证明数值算法的高精度性. 最后将数值算法推广到应用层面来动态模拟二维电缆的短时间行为.
在文献[1]的模型中,对于电缆上任意微元$\Delta s$,我们考虑了重力$\rho g \Delta s$, 浮力$-\rho_0 g \Delta s$,流体阻力$D$以及微元两端点的张力$T$ (如图 2(a)所示).现在 为了更加准确地描述电缆的运动轨迹,我们额外地考虑了流体升力$L$ (如图 2(b)所示)的作用.
为了简化问题,我们假设海底是平滑的,海水处于静止状态,且电缆线是由连续、弹性、非延展的材料构成.基于以上假设,我们来建立二维电缆的运动模型.首先我们建立坐标系 $xOy$ (如图 2 所示),其中$Ox$轴表示海底,直线$y=h$ (图中虚线)表示海面. 令$t$时刻描述电缆质点位置和张力的参数方程如下 $ x = x(s,t),y = y(s,t),T = T(s,t), $ (2.1) 其中,$s$ 是从$O$点算起的电缆总长度,$t$ 表示时间.
基于微元法思想,取悬垂电缆上任意一个微元并分析其受力情况(如图 2(b)所示). 现在我们需要将各个力分解到$x$,$y$方向并在这两个方向上建立模型. 因此,张力$T$与$x$轴之间的夹角 $\alpha(s,t) $ 需要表示成$s$和$t$的函数.
令 $$ \Delta x =x(s+ \Delta s ,t)-x(s,t),\Delta y =y(s+ \Delta s ,t)-y(s,t),$$ 则 ${\Delta s }= \sqrt{(\Delta x)^{2}+(\Delta y)^{2}} $. 根据几何关系(如图 3(a) 所示),可得 $ \cos\alpha(s,t) =\frac{x(s+\Delta s,t)-x(s,t)}{\Delta s}=\frac{\Delta x}{\Delta s}, \sin\alpha(s,t) =\frac{y(s+\Delta s,t)-y(s,t)}{\Delta s}=\frac{\Delta y}{\Delta s}. $ (2.2)
再令 $ \Delta s \rightarrow 0$,可得 $ \cos\alpha(s,t)=\frac{\partial x(s,t)}{\partial s},\ \sin\alpha(s,t) =\frac{\partial y(s,t)}{\partial s}. $ (2.3)
注2.1根据物理学原理以及图 3(a)可知, $x(s,t)<x(s+\Delta s,t),y(s,t)<y(s+\Delta s,t)$.即 $$0\leq \cos\alpha(s,t)=\frac{\partial x(s,t)}{\partial s}\leq 1,0\leq \sin\alpha(s,t) =\frac{\partial y(s,t)}{\partial s}\leq 1. $$
为了书写方便,我们记 $$\rho=\pi r^{2}\rho_{c},~~\rho_{0}=\pi r^{2}\rho_{w}, $$ 其中, $\rho_{c},\rho_{w}$ 表示电缆线和水的密度,$\rho,\rho_{0}$ 表示它们线密度(单位长度的质量),$ r=\frac {d}{2}$ 表示电缆的半径.
设$\vec{\mathbf{i}}$,$\vec{\mathbf{j}}$分别表示$x,y$轴正向的单位向量. 取坐标轴正向为正方向,则微元$\Delta s$ 的重力与浮力的合力为$-\Delta s (\rho-\rho_{0})g\vec{\mathbf{j}}$.
根据物理学原理[4]:绕流流体阻力$D$的大小与流体速度$U$的平方以及物体的有 效面积$A_p$成正比,方向与流体速度方向相同.同时,流体升力$L$的大小与流体速度 $U$的平方以及物体的有效面积$A$成正比,方向与流体速度方向垂直.即 $ D=C_{D}\frac {\rho_{f}}{2} U^2 A_p,L=C_{L}\frac {\rho_{f}}{2} U^2 A, $ (2.4) 其中,$ C_{D} $ 是总阻力系数,$ C_{L} $ 是总升力系数,$\rho_{f}$ 是流体密度, $U$ 是流体速度,$A_p$是物体沿速度方向的投影面积,$A$ 是物体在垂直于速度方向的投影面积.
根据(2.4)式以及图 3(a)--图 3(c),可得 $D=D_x+D_y$,升力 $L=L_x+L_y$. 其中, $ \ D_x=f_x\vec{\mathbf{i}}=-\frac{1}{2}C_D\rho_w d \bigg(\frac {\partial x}{\partial t}\bigg)^2 \Delta y \vec{\mathbf{i}},~~ L_x=g_x\vec{\mathbf{j}}=\frac{1}{2}C_L\rho_w d \bigg(\frac {\partial x}{\partial t}\bigg)^2 \Delta x \vec{\mathbf{j}}, $ (2.5) $ D_y=f_y\vec{\mathbf{j}}=\ \frac{1}{2}C_D\rho_w d \bigg(\frac {\partial y}{\partial t}\bigg)^2 \Delta x \vec{\mathbf{j}},~~ L_y=g_y\vec{\mathbf{i}}=-\frac{1}{2}C_L\rho_w d \bigg(\frac {\partial y}{\partial t}\bigg)^2 \Delta y \vec{\mathbf{i}}, $ (2.6) $D_x,L_x$分别为考虑微元仅在$x$方向上运动时所受的流体阻力和升力; $D_y,L_y$分别为考虑微元仅在$y$方向上运动时所受的流体阻力和升力.
类似的原理,在\textit{y}方向上有 $ T(s+\Delta s,t)\sin\alpha(s+\Delta s,t)-T(s,t)\sin\alpha(s,t)+f_y+g_x-\Delta s (\rho-\rho_{0})g =\Delta s\rho\frac{\partial^{2} y}{\partial t^{2}}. $ (2.10)
由(2.3)式以及(2.5)--(2.6)式,可得 $ T(s+\Delta s,t)\frac{\partial y}{\partial s}(s+\Delta s,t)-T(s,t) \frac{\partial y}{\partial s}(s,t)+f_y+g_x-\Delta s(\rho-\rho_{0})g =\Delta s\rho\frac{\partial^{2} y}{\partial t^{2}}. $ (2.11)
(2.11)式两端除以 $\Delta s$ 同时令 $\Delta s \rightarrow 0$,有 $ \frac{\partial }{\partial s}\bigg(T\frac{\partial y}{\partial s}\bigg) +\frac{1}{2}C_D\rho_w d \bigg(\frac {\partial y}{\partial t}\bigg)^2\frac{\partial x}{\partial s} +\frac{1}{2}C_L\rho_w d \bigg(\frac {\partial x}{\partial t}\bigg)^2\frac{\partial x}{\partial s} -(\rho-\rho_0)g=\rho\frac{\partial^{2} y}{\partial t^{2}}. $ (2.12)
除此以外,由(2.3)式可得电缆线的几何约束方程 $ \bigg(\frac{\partial x}{\partial s}\bigg)^{2}+\bigg(\frac{\partial y}{\partial s}\bigg)^{2}=1. $ (2.13)
由(2.9),(2.12)以及(2.13)式可得,二维海底电缆铺设模型可以写为 $ \left\{\begin{array}{lllll} \frac{\partial }{\partial s}\bigg(T\frac{\partial x}{\partial s}\bigg) -K_1\bigg(\frac {\partial x}{\partial t}\bigg)^2\frac{\partial y}{\partial s} -K_2\bigg(\frac {\partial y}{\partial t}\bigg)^2\frac{\partial y}{\partial s} =\rho\frac{\partial^{2} x}{\partial t^{2}} ,\\ \frac{\partial }{\partial s}\bigg(T\frac{\partial y}{\partial s}\bigg) +K_1\bigg(\frac {\partial y}{\partial t}\bigg)^2\frac{\partial x}{\partial s} +K_2\bigg(\frac {\partial x}{\partial t}\bigg)^2\frac{\partial x}{\partial s} -(\rho-\rho_0)g =\rho\frac{\partial^{2} y}{\partial t^{2}} ,\\ \bigg (\frac{\partial x}{\partial s}\bigg)^{2}+ \bigg(\frac{\partial y}{\partial s}\bigg)^{2}=1, \end{array}\right. $ (2.14) 其中,$K_1=\frac{1}{2}C_D\rho_w d,K_2=\frac{1}{2}C_L\rho_w d$.
除此以外,从实际背景出发,我们可以得到相应的初始条件和边界条件.
(1)初始条件
假设 $a(t)$ 表示$t$时刻悬垂电缆与海底接触点(TDP)的$x$ 坐标值, $b(t)$表示$t$时刻的电缆总长度,则有 $ x(s,t)\big|_{t=0}=0,\ \ \ \ \ y(s,t)\big|_{t=0}=s, $ (2.15) $ \frac{\partial x(s,t)}{\partial t}\Big|_{t=0}=0,\ \ \ \ \ \frac{\partial y(s,t)}{\partial t}\Big|_{t=0}=0, $ (2.16) $ T(s,t)\big|_{t=0}=(\rho-\rho_{0})gs,\ \ \ \ \ a(t)\big|_{t=0}=0, $ (2.17) $ b(t)=\left\{\begin{array}{ll} h+0.5a_0t^{2} , & 0\leq t<\frac {\sigma}{a_0},\\ h+ \frac {{\sigma}^2}{2a_0}+\sigma (t-\frac {\sigma}{a_0}),~~ & t \geq \frac {\sigma}{a_0}, \end{array}\right. $ (2.18) 其中,$\sigma$ 表示探测船最终的匀速速度,$a_0$ 表示探测船加速过程中的加速度, $h$ 表示海底的深度.
(2)边界条件 $ x(s,t)\big|_{s\leq a(t)}=s,\ \ \ y(s,t)\big|_{s\leq a(t)}=0,\ \ \ T(s,t)\big|_{s\leq a(t)}=f_{0}, $ (2.19) $ x(s,t)|_{s=b(t)} =\left\{\begin{array}{ll} 0.5a_0t^{2} ,& 0\leq t<\frac {\sigma}{a_0} ,\\ \frac {{\sigma}^2}{2a_0}+\sigma (t-\frac {\sigma}{a_0}) ,~~ &t \geq \frac {\sigma}{a_0}, \end{array}\right. $ (2.20) $ y(s,t)\big|_{s=b(t)}=h,\ \ \ \ \ \ \ T(s,t)\big|_{s=b(t)}=f(t). $ (2.21)
注2.2对(2.14)式第一式关于$s$在$(a(t),b(t))$上积分得 \begin{eqnarray*} &&T(b(t),t)\frac{\partial x(b(t),t)}{\partial s} \\ &=&T(a(t),t)\frac{\partial x(a(t),t)}{\partial s}+\int_{a(t)}^{b(t)}K_1 \bigg(\frac {\partial x}{\partial t}\bigg)^2\frac{\partial y}{\partial s}+K_2 \bigg(\frac {\partial y}{\partial t}\bigg)^2\frac{\partial y}{\partial s}+ \rho\frac{\partial^{2} x}{\partial t^{2}}{\rm d}s\\ & =&A_1 . \end{eqnarray*}
对(2.14)式第二式关于$s$在$(a(t),b(t))$上积分得 \begin{eqnarray*} &&T(b(t),t)\frac{\partial y(b(t),t)}{\partial s} \\ &=&T(a(t),t)\frac{\partial y(a(t),t)}{\partial s}-\int_{a(t)}^{b(t)}K_1 \bigg(\frac {\partial y} {\partial t}\bigg)^2\frac{\partial x}{\partial s}+K_2\bigg(\frac {\partial x} {\partial t}\bigg)^2\frac{\partial x}{\partial s} -(\rho-\rho_0)g+\rho\frac{\partial^{2} y}{\partial t^{2}}{\rm d}s \\ &=& A_2 . \end{eqnarray*}
再由(2.14)式第三式得 $ T(b(t),t)=\sqrt{A_1^2+A_2^2}. $ (2.22)
在给定的运动状态下,$T(b(t),t),T(a(t),t)$必定满足(2.22)式,即 模型中只需给定$T(b(t),t)$,$T(a(t),t)$ 其中一个即可.因此(2.19)与(2.21)式 中的最后两个边界条件只需给定其中一个即可.鉴于(2.19)式中最后一个条件是基于物 理意义的假设(事实上,本文是在此假设的基础上做的研究),而 (2.21)式中的最后一个 边界条件目前还是很多学者感兴趣的研究课题,因此本文将这两个条件同时列出.
在文献[7]中,Faltinsen从物理分析的角度给出了静态情况下悬垂电缆线上各点的纵坐标$x$与 $y$, $T$的解析关系.在文献[7]中假设条件下,本节的定理3.1利用数学的方法严格推导出了 与文献[7]等价的$x$,$y$,$T$与弧长参数$s$的解析关系.并将定理3.1的结果延伸到更一般的情形 (定理 3.2),提出了相应的解的存在条件.同时,通过找出静态 解的解析式,我们可以利用(4.2)--(4.8)式计算相同条件下的数值解,通过数值解与解析式的比较, 一定程度上可以体现数值格式(4.2)--(4.8)式的高精度性.
现在我们来考虑(2.14)式的一种特殊情况下的解的存在性. \begin{equation}\label{eqn 3.1} \left\{\begin{array}{llll} \frac{\rm d}{{\rm d}s}\bigg(T\frac{{\rm d}x}{{\rm d}s}\bigg)=0,\\ \frac{\rm d}{{\rm d}s}\bigg(T\frac{{\rm d}y}{{\rm d}s}\bigg)-(\rho-\rho_{0})g =0 ,\\ \bigg(\frac{{\rm d}x}{{\rm d}s}\bigg)^{2}+\bigg(\frac{{\rm d}y}{{\rm d}s}\bigg)^{2}=1 . \end{array}\right. \end{equation} (3.1)
注 3.1 很明显(3.1)式是(2.14)式在 $$ \frac{\partial x}{\partial t}\equiv 0,~~~\frac{\partial y}{\partial t}\equiv 0 , $$ $\forall s \in [x_*,s_*],t\in [0,+\infty)$下的简化模型.
首先,考虑电缆线全局$(0\leq s \leq b(t))$光滑的情况,此时相应的边界条件如下 \begin{equation}\label{eqn 3.2} T\frac{{\rm d}x}{{\rm d}s}\Big|_{s=x_*}=f_0,\ \ \ \ T\frac{{\rm d}y}{{\rm d}s}\Big|_{s=x_*}=0,\ \ \ \ x(x_*)=x_*,\ \ \ \ y(x_*)=0, \end{equation} (3.2) 其中,$x_*$ 是悬垂电缆与海底接触点(TDP)的 $x$ 坐标.
我们得到如下结论.
定理3.1在边界条件(3.2)下,方程组 (3.1 )存在唯一全局光滑解 \begin{equation}\label{eqn 3.3} \left\{\begin{array}{llll} x(s)=x_*+\frac{1}{K}\ln[K(s-x_*)+\sqrt{K^2(s-x_*)^2+1}],\\ y(s)=\frac{1}{K}[\sqrt{K^2(s-x_*)^2+1}-1] ,\\ T(s)=f_0\sqrt{K^2(s-x_*)^2+1} , \end{array}\right. \end{equation} (3.3) 其中,$K=\frac {(\rho-\rho_0)g}{f_0}$.
证由(3.1)式的第一个方程关于$s$从$x_*$到$s$积分以及 (3.2)式,可得 \begin{equation}\label{eqn 3.4} T\frac{{\rm d}x}{{\rm d}s}=T(x_*)\frac{{\rm d} x(x_*)}{{\rm d}s}=f_0. \end{equation} (3.4)
由(3.1)式的第二个方程关于$s$从$x_*$到$s$积分以及(3.2)式,可得 \begin{equation}\label{eqn 3.5} T\frac{{\rm d}y}{{\rm d}s}=T(x_*)\frac{{\rm d}y(x_*)}{{\rm d}s}+(\rho-\rho_0)g(s-x_*)=(\rho-\rho_0)g(s-x_*). \end{equation} (3.5)
由(3.1)式的第三个方程以及(3.4)--(3.5)式,可得 \begin{equation}\label{eqn 3.6} T(s)=f_0\sqrt{1+K^2(s-x_*)^2}, \end{equation} 其中,$K=\frac{(\rho-\rho_0)g}{f_0}$. (3.6)
将(3.6)式代入(3.4)式,可得 \begin{equation}\label{eqn 3.7} \frac{{\rm d}x}{{\rm d}s}=\frac{1}{\sqrt{1+K^2(s-x_*)^2}}, \end{equation} (3.7) 即 \begin{eqnarray}\label{eqn 3.8} x(s)&=&x_*+\int_{x_*}^{s}\frac{1}{\sqrt{1+K^2(s-x_*)^2}}{\rm d}s \nonumber\\ &=&x_*+\frac{1}{K}\ln[K(s-x_*)+\sqrt{K^2(s-x_*)^2+1}]. \end{eqnarray} (3.8)
将(3.6)式代入(3.5)式,可得 \begin{equation}\label{eqn 3.9} \frac{{\rm d}y}{{\rm d}s}=\frac{K(s-x_*)}{\sqrt{1+K^2(s-x_*)^2}}, \end{equation} 即 \begin{equation}\label{eqn 3.10} y(s)=\int_{x_*}^{s}\frac{K(s-x_*)}{\sqrt{1+K^2(s-x_*)^2}}{\rm d}s =\frac{1}{K}[\sqrt{K^2(s-x_*)^2+1}-1]. \end{equation} (3.9)
在施工过程中,最理想的状态是让电缆全局光滑,因为只有全局光滑的电缆线才能使海底接触点 $(a(t),0)$处的张力为$0$.然而,由(3.2)式可以看出,电缆线全局光滑的 条件是很苛刻的,而且电缆线施工过程中在船速太大的情况下, 电缆绳大多数不会保持全局光滑姿态,而是在$(a(t),0)$处产生尖点. 因此,研究电缆线局部光滑情况下 解的存在条件有着很重要的意义.
接下来,我们考虑电缆线局部光滑$(a(t)\leq s \leq b(t))$的情况.此时相应的边界条件如下 \begin{equation}\label{eqn 3.11} x(x_*)=x_*,\ \ \ \ \ \ \ y(x_*)=0,\ \ \ \ \ \ \ T(x_*)=f_0,\ \ \ \ \ \ \ y(s_*)=h, \end{equation} (3.10) 其中,$s_*$满足 $x(s_*)=x_0$.
引理3.1假设$K(x_0-x_*)\geq\ln 2$,以及 $\frac{1}{2}{\rm e}^{K(x_0-x_*)}-\frac{1}{2}{\rm e}^{-K(x_0-x_*)}<\sqrt{K^2h^2+2Kh}$,则函数 $$ \varphi(u)=\frac{1}{2}(\sqrt{1-u^2}+1){\rm e}^{\frac{K(x_0-x_*)}{u}} -\frac{u^2}{2(\sqrt{1-u^2}+1){\rm e}^{\frac{K(x_0-x_*)}{u}}}-\sqrt{(Kh+1)^2-u^2} $$ 在$(0,1)$上存在唯一的零点.
证 令$m=K(x_0-x_*)\geq \ln 2$,当$u\in (0,1)$时,通过简单计算可得 \begin{eqnarray*}\label{E:1.1} \varphi'(u)&=&-\frac{1}{2}\bigg[{\rm e}^{\frac{m}{u}}+\frac{u^2}{(\sqrt{1-u^2}+1)^2 {\rm e}^{\frac{m}{u}}}\bigg]\bigg[\frac{u}{\sqrt{1-u^2}}+\frac{(\sqrt{1-u^2}+1)m}{u^2}\bigg] \\ &&-\frac{u}{(\sqrt{1-u^2}+1){\rm e}^{\frac{m}{u}}} +\frac{u}{\sqrt{(Kh+1)^2-u^2}} \nonumber \\ &\leq & -\frac{1}{2}{\rm e}^{\frac{m}{u}}\frac{u}{\sqrt{1-u^2}}+\frac{u}{\sqrt{(Kh+1)^2-u^2}} \nonumber \\ &\leq &-\frac{1}{2}({\rm e}^{\frac{m}{u}}-2)\frac{u}{\sqrt{1-u^2}} \\ & <& -\frac{1}{2}({\rm e}^{m}-2)\frac{u}{\sqrt{1-u^2}} \\& <& -\frac{1}{2}({\rm e}^{ln2}-2)\frac{u}{\sqrt{1-u^2}} \\ & =& 0 . \end{eqnarray*}
因此,$\varphi(u)$在$(0,1)$上单调递减.
此外还有 $$ \lim \limits_{u\rightarrow 0}\varphi(u)=+\infty, \lim\limits_{u\rightarrow 1}\varphi(u)=\frac{1}{2}{\rm e}^{K(x_0-x_*)} -\frac{1}{2}{\rm e}^{-K(x_0-x_*)}-\sqrt{K^2h^2+2Kh}< 0. $$
由介值定理,可得$\varphi(u)$在$(0,1)$上存在唯一的零点.
定理 在引理3.1的条件以及边界条件(3.11)下,方程组(3.1)存在唯一局部光滑解 \begin{equation}\label{eqn 3.12} \left\{\begin{array}{llll} x(s)=x_*+\frac{1}{K}\frac{{\rm d}x(x_*)}{{\rm d}s}\ln\frac{\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*) +\sqrt{[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)]^2+(\frac{{\rm d}x(x_*)}{{\rm d}s})^2}}{\frac{{\rm d}y(x_*)}{{\rm d}s}+1}, \\ y(s)=\frac{1}{K} \Bigg[\sqrt{\bigg[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)\bigg]^2+\bigg(\frac{{\rm d}x(x_*)}{{\rm d}s}\bigg)^2}-1\Bigg] , \\ T(s)=f_0\sqrt{\bigg(\frac{{\rm d}x(x_*)}{{\rm d}s}\bigg)^2+ \bigg[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)\bigg]^2}, \end{array}\right. \end{equation} (3.12) 其中,$K=\frac {(\rho-\rho_0)g}{f_0},\frac{{\rm d}x(x_*)}{{\rm d}s}$ 是引理3.1中$\varphi(u)$的唯一零点,$\frac{{\rm d}y(x_*)}{{\rm d}s}=\sqrt{1-(\frac{{\rm d}x(x_*)}{{\rm d}s})^2}$.
证由(3.1)式第一个方程以关于$s$从$x_*$到$s$及(3.11)式,可得 \begin{equation}\label{eqn 3.13} T\frac{{\rm d}x}{{\rm d}s}=T(x_*)\frac{{\rm d}x(x_*)}{{\rm d}s}=f_0\frac{{\rm d}x(x_*)}{{\rm d}s}. \end{equation} (3.13)
由(3.1)式第二个方程关于$s$从$x_*$到$s$以及(3.11)式,可得 \begin{equation} \label{eqn 3.14} T\frac{{\rm d}y}{{\rm d}s}=T(x_*)\frac{{\rm d}y(x_*)}{{\rm d}s}+(\rho-\rho_0)g(s-x_*)=f_0\frac{{\rm d}y(x_*)}{{\rm d}s}+(\rho-\rho_0)g(s-x_*). \end{equation} (3.14)
由(3.1)式第三个方程以及(3.11)式,可得 \begin{equation}\label{eqn 3.15} T(s)=f_0\sqrt{ \bigg[\frac{{\rm d}x(x_*)}{{\rm d}s}\bigg]^2+\bigg[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)\bigg]^2}, \end{equation} (3.15) 其中,$K=\frac{(\rho-\rho_0)g}{f_0}$.
将(3.15)式代入(3.13)式,可得 \begin{equation}\label{eqn 3.16} \frac{{\rm d}x}{{\rm d}s}=\frac{\frac{{\rm d}x(x_*)}{{\rm d}s}}{\sqrt{[\frac{{\rm d}x(x_*)}{{\rm d}s}]^2+[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)]^2}}, \end{equation} (3.16) 即 \begin{eqnarray}\label{eqn 3.17} x(s)&=&x_*+\int_{x_*}^{s}\frac{\frac{{\rm d}x(x_*)}{{\rm d}s}}{\sqrt{[\frac{{\rm d}x(x_*)}{{\rm d}s}]^2+[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)]^2}}{\rm d}s \nonumber\\ &=&x_*+\frac{\frac{{\rm d}x(x_*)}{{\rm d}s}}{K}\ln\frac{\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)+\sqrt{[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)]^2+[\frac{{\rm d}x(x_*)}{{\rm d}s}]^2}}{\frac{{\rm d}y(x_*)}{{\rm d}s}+1}. \end{eqnarray} (3.17)
将(3.15)式代入(3.14)式,可得 \begin{equation}\label{eqn 3.18} \frac{{\rm d}y}{{\rm d}s}=\frac{\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)}{\sqrt{[\frac{{\rm d}x(x_*)}{{\rm d}s}]^2+[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)]^2}}, \end{equation} (3.18) 即 \begin{eqnarray}\label{eqn 3.19} y(s)&=&\int_{x_*}^{s}\frac{\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)}{\sqrt{[\frac{{\rm d}x(x_*)}{{\rm d}s}]^2+[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)]^2}}{\rm d}s \nonumber \\ &=&\frac{1}{K} \Bigg\{\sqrt{\bigg[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s-x_*)\bigg]^2+ \bigg[\frac{{\rm d}x(x_*)}{{\rm d}s}\bigg]^2}-1\Bigg\}. \end{eqnarray} (3.19)
由(3.17)式可得 $$ x_0=x(s_*)=x_*+\frac{\frac{{\rm d}x(x_*)}{{\rm d}s}}{K}\ln\frac{\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s_*-x_*)+\sqrt{[\frac{{\rm d}y(x_*)}{{\rm d}s}+K(s_*-x_*)]^2+[\frac{{\rm d}x(x_*)}{{\rm d}s}]^2}}{\frac{{\rm d}y(x_*)}{{\rm d}s}+1}, $$ 则有 \begin{equation}\label{eqn 3.20} s_*=x_*+\frac{1}{2K}\bigg(\frac{{\rm d}y(x_*)}{{\rm d}s}+1\bigg){\rm e}^{\frac{K(x_0-x_*)}{\frac{{\rm d}x(x_*)}{{\rm d}s}}}-\frac{(\frac{{\rm d}x(x_*)}{{\rm d}s})^2}{2K(\frac{{\rm d}y(x_*)}{{\rm d}s}+1){\rm e}^{\frac{K(x_0-x_*)}{\frac{{\rm d}x(x_*)}{{\rm d}s}}}}-\frac{1}{K}\frac{{\rm d}y(x_*)}{{\rm d}s}. \end{equation} (3.20)
由(3.19)式,令 $y(s_{*})=h$,则有 \begin{equation}\label{eqn 3.21} s_{*}=x_*+\frac{1}{K}\sqrt{(Kh+1)^2-\bigg(\frac{{\rm d}x(x_*)}{{\rm d}s}\bigg)^2}-\frac{1}{K}\frac{{\rm d}y(x_*)}{{\rm d}s}. \end{equation} (3.21)
令 $u=\frac{{\rm d}x(x_*)}{{\rm d}s},$则根据(3.1)式的第一个方程可得$\sqrt{1-u^2}=\frac{{\rm d}y(x_*)}{{\rm d}s}$. 再由(3.20)--(3.21)式可得 $$ \varphi(u)=\frac{1}{2}(\sqrt{1-u^2}+1){\rm e}^{\frac{K(x_0-x_*)}{u}}-\frac{u^2}{2(\sqrt{1-u^2}+1){\rm e}^{\frac{K(x_0-x_*)}{u}}}-\sqrt{(Kh+1)^2-u^2}=0. $$
由注2.1以及引理3.1可知$\frac{{\rm d}x(x_*)}{{\rm d}s} $是$\varphi(u)$的唯一零点, $\frac{{\rm d}y(x_*)}{{\rm d}s}=\sqrt{1-(\frac{{\rm d}x(x_*)}{{\rm d}s} )^2} $.
自由边界问题(2.14)--(2.21)相当复杂,以至于目前还无法直接表达它的解. 然而,在电缆通信以及石油探测方面的应用要求我们至少能找到它的近似解. 因此,本节采用数值计算的方法[13]来近似求解自由边界问题(2.14)--(2.21).
首先,我们将弧长参数$s$ 以及时间参数$t$ 离散化如下 \begin{equation}\label{eqn 4.1} 0=s_0<s_1<s_2<\cdots<s_{I_0+N},0=t_0<t_1<\cdots <t_n<\cdots <t_N; \end{equation} (4.1)
在节点 $(s_i,t_{n+1})(n=1,2,\cdots,N-1;i=1,2,\cdots,I_0+n)$处,将方程组(2.14)离散化[8]如下 \begin{eqnarray}\label{eqn 4.2} F_i(\mathbf{U})& =&T_{i+\frac{1}{2}}^{n+1}\alpha_{i+1}(x_{i+1}^{n+1}-x_{i}^{n+1})-T_{i-\frac{1}{2}}^{n+1}\beta_{i+1}(x_{i}^{n+1}-x_{i-1}^{n+1}) \nonumber\\ & & -K_1\bigg(\frac{x_{i}^{n+1}-x_{i}^{n}}{\Delta t_{n+1}}\bigg)^2 \frac{y_{i+1}^{n+1}-y_{i}^{n+1}}{\Delta s_{i+1}}-K_2 \bigg(\frac{y_{i}^{n+1}-y_{i}^{n}}{\Delta t_{n+1}}\bigg)^2\frac{y_{i+1}^{n+1}-y_{i}^{n+1}}{\Delta s_{i+1}} \nonumber\\ && -\rho[\gamma_{n+1}(x_{i}^{n+1}-x_{i}^{n})-\theta_{n+1}(x_{i}^{n}-x_{i}^{n-1})] \nonumber\\ &=&0, \end{eqnarray} (4.2) \begin{eqnarray}\label{eqn 4.3} G_i(\mathbf{U})& =&T_{i+\frac{1}{2}}^{n+1}\alpha_{i+1}(y_{i+1}^{n+1}-y_{i}^{n+1})-T_{i-\frac{1}{2}}^{n+1}\beta_{i+1}(y_{i}^{n+1}-y_{i-1}^{n+1}) \nonumber\\ && +K_1\bigg(\frac{y_{i}^{n+1}-y_{i}^{n}}{\Delta t_{n+1}}\bigg)^2 \frac{x_{i+1}^{n+1}-x_{i}^{n+1}}{\Delta s_{i+1}}+K_2 \bigg(\frac{x_{i}^{n+1}-x_{i}^{n}}{\Delta t_{n+1}}\bigg)^2\frac{x_{i+1}^{n+1}-x_{i}^{n+1}}{\Delta s_{i+1}} \nonumber\\ & & -(\rho-\rho_0)g-\rho[\gamma_{n+1}(y_{i}^{n+1}-y_{i}^{n})-\theta_{n+1}(y_{i}^{n}-y_{i}^{n-1})] \nonumber\\ &=&0, \end{eqnarray} (4.3) \begin{equation}\label{eqn 4.4} \Psi_i({\mathbf U})=(x_{i+1}^{n+1}-x_{i}^{n+1})^2+(y_{i+1}^{n+1}-y_{i}^{n+1})^2-(\Delta s_{i+1})^2=0 , \end{equation} (4.4) 其中, $$ \alpha_{i}=\frac {2}{\Delta s_{i+1}(\Delta s_{i+1}+\Delta s_{i})},~~\beta_{i}=\frac {2}{\Delta s_{i}(\Delta s_{i+1}+\Delta s_{i})}, $$ $$ \gamma_{n+1}=\frac {2}{\Delta t_{n+1}(\Delta t_{n+1}+\Delta t_{n})},~~\theta_{n+1}=\frac {2}{\Delta t_{n}(\Delta t_{n+1}+\Delta t_{n})}, $$ $$ {\mathbf U}=(x_1^{n+1},x_2^{n+1},\cdots,x_{I_0+n}^{n+1},y_1^{n+1},y_2^{n+1},\cdots,y_{I_0+n}^{n+1},T_{\frac{3}{2}}^{n+1},T_{\frac{5}{2}}^{n+1},\cdots,T_{I_0+n+{\frac{1}{2}}}^{n+1})^{T} . $$
接下来,我们在初始节点$(s_i,t_0)$处,将(2.15)--(2.18)式离散化如下 \begin{equation}\label{eqn 4.5} x_{i}^{0}=0 ,y_{i}^{0}=s_{i},x_{i}^{1}-x_{i}^{0}=0,y_{i}^{1}-y_{i}^{0}=0, T_{i+\frac{1}{2}}^{0}=(\rho-\rho_0)gs_{i+\frac{1}{2}} , \ i=0,1,\cdots,I_{0}. \end{equation} (4.5)
设 $r$表示已经沉到海底的节点总数,则在边界节点处将(2.19)--(2.21)式离散化如下 \begin{equation}\label{eqn 4.6} y_{r}^{n+1}=h ,T_{r+\frac{1}{2}}^{n+1}=f_{0} ,x_{I_0+n}^{n+1}=xup(t_{n+1}) , y_{I_0+n}^{n+1}=0 ,T_{I_0+n+\frac{1}{2}}^{n+1}=f(t_{n+1}) , \ n=0,1,\cdots,N-1. \end{equation} (4.6)
注 4.1 对于任意时刻 $t$,怎样获取$r$的值?首先, 当$t=0$时,令 $r=0$,并设定任意小的常量$\varepsilon >0$.假设在$t=t_{n}$时刻有 $r=r_{*}$.则在$t=t_{n+1}$时刻,我们比较 $y_i^{n+1}(i=r+1,r+2,\cdots,I_0+n)$ 与 $\varepsilon$ 的大小,并令 $l=\overline{\overline{\{i_{*}|y_{i_*}^{n+1}<\varepsilon, i_{*}=r+1,r+2,\cdots,I_0+n\}}}. $ 在$t=t_{n+1}$时刻 $r$的值为$r=r_{*}+l$. 其中,对任意集合$A$,$\overline{\overline{A}}$ 表示$A$中元素的个数.
由(4.2)--(4.4)式可知,在$t_{n+1}$时刻,当 $i=I_0+n$时, $x_{I_0+n}^{n-1},y_{I_0+n}^{n-1},z_{I_0+n}^{n-1}$ 是未知的(因为$n-1$时刻,第$I_0+n$个节点还未沉放入海底). 因此,我们需要给这几个变量取定一个合适的近似值,根据几何关系有 \begin{equation}\label{eqn 4.7} x_{I_0+n}^{n-1}=x_{I_0+n-1}^{n-1}+\frac {\Delta s_{I_0+n}}{\Delta s_{I_0+n-1}}(x_{I_0+n-1}^{n-1}-x_{I_0+n-2}^{n-1}), \end{equation} (4.7) \begin{equation}\label{eqn 4.8} y_{I_0+n}^{n-1}=y_{I_0+n-1}^{n-1}+\frac {\Delta s_{I_0+n}}{\Delta s_{I_0+n-1}}(y_{I_0+n-1}^{n-1}-y_{I_0+n-2}^{n-1}). \end{equation} (4.8)
这种近似处理以后,$x_{I_0+n}^{n-1},y_{I_0+n}^{n-1}$就可以参与计算了.
由(2.17)和(2.18)式可知,在每个时刻 $t$,方程组(2.14)式的弧长参数取值范围为 $a(t) \leq s \leq b(t)$,其中$a(t)$ 随着沉入海底的部分越来越多而逐渐增大, $ b(t)$随着探测船行走过程中释放缆绳也逐渐增大,这就涉及到一个自由边界计算的问题[2]. 自由边界是个很复杂的问题.在这篇论文中,上边界$b(t)$由(2.18)式定义,下 边界$a(t)$是未知的.怎样取$a(t)$的值呢?本文采用如下方法: 为了让已经坠入海底的节点退出计算 ,我们在(4.2)--(4.4)式中将 $x_{i}^{n+1},y_{i}^{n+1},T_{i+\frac{1}{2}}^{n+1}~(0 \leq i\leq r)$作 为边界值,并令 $x_{i}^{n+1}=x_{i}^{n},y_{i}^{n+1}=y_{i}^{n}=0,$ $T_{i+\frac{1}{2}}^{n+1}=f_0~(0 \leq i\leq r). $ 同时,为了让刚刚从船上 放入海水中的节点参与计算,我们每一个时刻在缆绳顶端增加一个节点,即 增加$x_{I_0+n}^{n+1},y_{I_0+n}^{n+1},T_{I_0+n+\frac{1}{2}}^{n+1}$ 三个变量.除此以外,我们还用前面注释4.1的方法来取$r$的值,故有,$a(t_{n+1})=x_{r}^{n+1}$.
差分格式的离散化已经完成,下面我们采用迭代算法来求解非线性方程组(4.2)--(4.4)式,具体过程如下.
在$t=t_{n+1}$时刻,取如下初始迭代值: $ x_{i}^{(0)}=x_{i}^{n},$ $y_{i}^{(0)}=y_{i}^{n},$ $T_{i+\frac{1}{2}}^{(0)} =T_{i+\frac{1}{2}}^{n}$,$ i=r+1,\cdots,$ $I_0+n. $
(1) $x$迭代计算 \begin{eqnarray}\label{eqn 4.9} &&T_{i+\frac{1}{2}}^{(s)}\alpha_{i+1}x_{i+1}^{(s+1)}-(T_{i+\frac{1}{2}}^{(s)}\alpha_{i+1}+T_{i-\frac{1}{2}}^{(s)}\beta_{i+1}+\rho\gamma_{n+1})x_{i}^{(s+1)}+T_{i-\frac{1}{2}}^{(s)}\beta_{i+1}x_{i-1}^{(s+1)} \nonumber\\ &=&K_1\bigg(\frac{x_{i}^{(s)}-x_{i}^{n}}{\Delta t_{n+1}}\bigg)^2\frac{y_{i+1}^{(s)} -y_{i}^{(s)}}{\Delta s_{i+1}}+K_2 \bigg(\frac{y_{i}^{(s)}-y_{i}^{n}}{\Delta t_{n+1}}\bigg)^2\frac{y_{i+1}^{(s)}-y_{i}^{(s)}}{\Delta s_{i+1}} \nonumber\\ &&-\rho\gamma_{n+1}x_{i}^{n}-\rho\theta_{n+1}(x_{i}^{n}-x_{i}^{n-1}). \end{eqnarray} (4.9)
$x_{i}^{(s)},y_{i}^{(s)},T_{i+\frac{1}{2}}^{(s)}(i=r+1,\cdots,I_0+n)$ 的值已经在前面的迭代中算得,那么(4.9)式就是关于 $x_{i}^{(s+1)}$ $(i=r+1,\cdots,I_0+n)$ 的线性方程组.求解这样一个代数方程组, 可以得到$x_{i}^{(s+1)}$ $(i=r+1,\cdots,I_0+n) $的值.
(2) $y$迭代计算 \begin{equation}\label{eqn 4.10} y_{i+1}^{(s+1)}-y_{i}^{(s+1)}=\sqrt{|(\Delta s_{i+1})^2-(x_{i+1}^{(s+1)}-x_{i}^{(s+1)})^2|}. \end{equation} (4.10)
$x_{i}^{(s+1)}$的值已由 (4.9)式获得,则(4.10) 式就成为关于$y_{i}^{(s+1)}$ $(i=r+1,\cdots,I_0+n)$的线性方程组.这样就可以计算出$y_{i}^{(s+1)}$ $ (i=r+1,\cdots,I_0+n)$ .
(3) $T$迭代计算 \begin{eqnarray}\label{eqn 4.11} &&\alpha_{i+1}(y_{i+1}^{(s+1)}-y_{i}^{(s+1)})T_{i+\frac{1}{2}}^{(s+1)}-\beta_{i+1}(y_{i}^{(s+1)}-y_{i-1}^{(s+1)})T_{i-\frac{1}{2}}^{(s+1)} \nonumber\\ &=&-K_1\bigg(\frac{y_{i}^{(s+1)}-y_{i}^{n}}{\Delta t_{n+1}}\bigg)^2\frac{x_{i+1}^{(s+1)} -x_{i}^{(s+1)}}{\Delta s_{i+1}}-K_2 \bigg(\frac{x_{i}^{(s+1)}-x_{i}^{n}}{\Delta t_{n+1}}\bigg)^2\frac{x_{i+1}^{(s+1)}-x_{i}^{(s+1)}}{\Delta s_{i+1}} \nonumber\\ &&+(\rho-\rho_0)g+\rho[\gamma_{n+1}(y_{i}^{(s+1)}-y_{i}^{n})-\theta_{n+1}(y_{i}^{n}-y_{i}^{n-1})]. \end{eqnarray} (4.11)
$x_{i}^{(s+1)},y_{i}^{(s+1)}$ 已由(4.9)--(4.10)式获得,则 (4.11)式是关于 $T_{i+\frac{1}{2}}^{(s+1)}$ $(i=r+1,\cdots,I_0+n)$ 的线性方程组. 因此,可以计算出 $T_{i+\frac{1}{2}}^{(s+1)}$,$i=r+1,\cdots,I_0+n$.
现在 $x_{i}^{(s+1)},y_{i}^{(s+1)},T_{i+\frac{1}{2}}^{(s+1)}$的值已获得, 重复 (4.9)--(4.11)式 依此取$s=0,1,2,\cdots ,M$, 经过有限次迭代,得到$t=t_{n+1}$时刻的数值解如下 $$ x_{i}^{n+1}=x_{i}^{(M)},y_{i}^{n+1}=y_{i}^{(M)}, T_{i+\frac{1}{2}}^{n+1}=T_{i+\frac{1}{2}}^{(M)},~ i=r+1,r+2,\cdots,I_0+n. $$
由定理3.1和定理3.2,我们知道,在解的存在性条件下,(3.3)和 (3.12)式是 (3.1)式 加上某些特定边界条件的解.现在为了检验我们提出的算法的有效性, 我们用此格式来数值模拟静态方程 (3.1)的解,并与其解析解(3.3),(3.12)进行对比, 通过下面的两组数值实验可以看出格式(4.2)--(4.8) 有很高的精度.
试验1 考虑方程组(3.1)在边界条件(3.2)下的解.各物理参数取值如表 1所示.模拟结果见图 4.
试验2 考虑方程组(3.1)在边界条件(3.11)下的解. 各物理参数取值如表 1所示.模拟结果见图 5.
二维电缆短时间实际应用模拟. 由图 4--图 5可以看出,我们提出的数值算法有很高的精度. 下面我们利用该数值算法来动态模拟二维电缆铺设过程中电缆线的运动轨迹 (定解问题(2.14)--(2.21)). 各参数取值如表 1 所示.模拟结果如图 6. 其中,彩色线表示电缆线.
2015, Vol. 35 Issue (4): 780-793
PDF (623 KB)
2015, Vol. 35 Issue (4): 780-793
PDF (623 KB)