2013年全国水动力学研讨会：海面油污扩散运动的力学建模及相关数值研究

二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 海面油污扩散的力学建模及相关数值研究 张大鹏,史倩,陈瑜,谢锡麟 (复旦大学力学与工程科学系,上海,200433) 摘要:本文基于几何形态为曲面的连续介质有限变形理论,建立海面油污扩散过程的动力 学模型。模型中视海面上油污的几何形态为曲面,通过面密度刻画油污层的实际厚度;油污 面可相对海面发生相对运动,二者之间可以考虑摩擦作用;可考虑油污面内的表面张力作用 以及摩擦作用;以及考虑油污面的重力作用等。就此模型,本文通过 Lagrange观点进行求解 并获得了若干油污扩散的特征。 关鍵词:几何形态为曲面的连续介质有限变形理论;海面油污扩散;油污面密度 1引言 连续介质根据其几何结构的差异可以分为两种,用现代几何学的观点,第一种可以看作 是 Euclid流形而另外一种则可以看作是欧氏空间中的二维曲面及 Riemann流形。谢锡麟等 研究了几何形态对应为 Euclid流形和 Riemann流形的连续介质的有限变形理论,指出二者 力学及数学力学机制上的差异。几何形态为曲面( Riemann流形)的连续介质的有限变形理论 按一般有限变形理论3(实际对应 Euclid流形)平行发展,主要包括构型构造及曲线坐标 系,变形梯度及其基本性质,变形刻画,输运方程以及守恒律控制方程等 二维流动一般是指在三维欧氏空间中的平面或曲面上的运动。物理上来讲,如果连续介 质的流动限制在薄层中,那么都可以看作是一种二维流动。流体力学中肥皂膜的流动应该可 以看作是经典的二维流动。而平坦肥皂膜对于之正交的圆柱体的钝体扰流也已经有前人做过 详细的研究。当前对肥皂膜的流动显示一般采用光干涉技术来实现。而肥皂膜流动的涡结构 般都与膜的面密度,或者曲面的厚度相耦合。密度相对于其平均值变化的最大振幅可以达 到平均值的20%,但是密度变化的机制鲜有人涉足。厚度足够薄的膜自身的有限振幅震动也可 以看作是二维流动。已有前人在此领域做了很多研究,相关的理论将会对纳米技术中膜理论 的应用产生正面的作用。不论是限制在薄层内的流动,还是薄膜自身的有限振幅震动,他们 的几何性质都可以看作是三维欧氏空间中的二维曲面,也就是三维欧氏空间中的二维 Riemann 流形。而 Euclid流形与 Riemann流形几何结构上的差别也反映了二维流体所满足的有限变形 理论与普通三维连续介质有限变形理论的差别 石油及其炼制品在开采、炼制、贮运和使用过程中进入海洋环境而造成的污染是目前 种世界性的严重的海洋污染。石油入海后即发生一系列复杂变化,包括扩散,蒸发,溶解 乳化,光化学氧化,微生物氧化,沉降,形成沥青球,以及沿着食物链转移等过程。入海石 本研究受国家自然科学基金面上项目(11172069),上海市教委2011年上海高校本科重点教学改革项目的资 助。通讯作者:谢锡麟复旦大学力学与工程科学系:Emal:xiexin@fudanedu.cn

第二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 - 1 - 海面油污扩散的力学建模及相关数值研究* 张大鹏, 史倩，陈瑜，谢锡麟* (复旦大学 力学与工程科学系，上海，200433) 摘要：本文基于几何形态为曲面的连续介质有限变形理论，建立海面油污扩散过程的动力 学模型。模型中视海面上油污的几何形态为曲面，通过面密度刻画油污层的实际厚度；油污 面可相对海面发生相对运动，二者之间可以考虑摩擦作用；可考虑油污面内的表面张力作用 以及摩擦作用；以及考虑油污面的重力作用等。就此模型，本文通过 Lagrange 观点进行求解， 并获得了若干油污扩散的特征。 关键词：几何形态为曲面的连续介质有限变形理论；海面油污扩散；油污面密度 1 引言 连续介质根据其几何结构的差异可以分为两种，用现代几何学的观点，第一种可以看作 是 Euclid 流形而另外一种则可以看作是欧氏空间中的二维曲面及 Riemann 流形。谢锡麟等 研究了几何形态对应为 Euclid 流形和 Riemann 流形的连续介质的有限变形理论[1]，指出二者 力学及数学力学机制上的差异。几何形态为曲面（Riemann 流形）的连续介质的有限变形理论 [1]按一般有限变形理论[2,3]（实际对应 Euclid 流形）平行发展，主要包括构型构造及曲线坐标 系，变形梯度及其基本性质，变形刻画，输运方程以及守恒律控制方程等。 二维流动一般是指在三维欧氏空间中的平面或曲面上的运动。物理上来讲，如果连续介 质的流动限制在薄层中，那么都可以看作是一种二维流动。流体力学中肥皂膜的流动应该可 以看作是经典的二维流动。而平坦肥皂膜对于之正交的圆柱体的钝体扰流也已经有前人做过 详细的研究。当前对肥皂膜的流动显示一般采用光干涉技术来实现。而肥皂膜流动的涡结构 一般都与膜的面密度，或者曲面的厚度相耦合。密度相对于其平均值变化的最大振幅可以达 到平均值的20%，但是密度变化的机制鲜有人涉足。厚度足够薄的膜自身的有限振幅震动也可 以看作是二维流动。已有前人在此领域做了很多研究，相关的理论将会对纳米技术中膜理论 的应用产生正面的作用。不论是限制在薄层内的流动，还是薄膜自身的有限振幅震动，他们 的几何性质都可以看作是三维欧氏空间中的二维曲面，也就是三维欧氏空间中的二维Riemann 流形。而Euclid 流形与Riemann 流形几何结构上的差别也反映了二维流体所满足的有限变形 理论与普通三维连续介质有限变形理论的差别。 石油及其炼制品在开采、炼制、贮运和使用过程中进入海洋环境而造成的污染是目前一 种世界性的严重的海洋污染。石油入海后即发生一系列复杂变化，包括扩散，蒸发，溶解， 乳化，光化学氧化，微生物氧化，沉降，形成沥青球，以及沿着食物链转移等过程。入海石 *本研究受国家自然科学基金面上项目（11172069），上海市教委 2011 年上海高校本科重点教学改革项目的资 助。通讯作者：谢锡麟 复旦大学 力学与工程科学系；Email: xiexilin@fudan.edu.cn

5二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 油首先在重力、惯性力、摩擦力和表面张力的作用下,在海洋表面迅速扩展成薄膜,进而在 风浪和海流作用下被分割成大小不等的块状或带状油膜,随风漂移扩散,速度约为风速的3% 石油中的氮、硫、氧等非烃组分是表面活性剂,促进石油扩散。本文将只考虑石油油污在海 面的扩散运动。 图1航拍海面油污扩散图 本文基于几何形态为曲面的连续介质有限变形理论,对海面上油污的扩散运动建立力学 模型,并按 Lagrange观点进行数值求解,研究海面对海面上油污的扩散运动在几种不同情况 下的运动进行了模拟 2海面油污扩散的动力学模型 本部分主要基于当前物理构型对应之曲线坐标系显含时间的有限变形理论,具体建立海 面油污扩散的动力学模型 21连续性方程 连续性方程的积分形式为 其中ρ为面密度。根据第一类面输运可得其微分形式为 p+b=0 其中速度的散度6=V·v=V 2.2动量方程 动量守恒方程可以表示为 =E+++F+E 式中,F表示二维流体的表面张力,F表示二维流体的内压力,F表示二维流体的黏滞 力,Fm表示二维流体的外压力,Fm表示二位流体受到的面力

第二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 - 2 - 油首先在重力、惯性力、摩擦力和表面张力的作用下，在海洋表面迅速扩展成薄膜，进而在 风浪和海流作用下被分割成大小不等的块状或带状油膜，随风漂移扩散，速度约为风速的3%。 石油中的氮、硫、氧等非烃组分是表面活性剂，促进石油扩散。本文将只考虑石油油污在海 面的扩散运动。 图 1 航拍海面油污扩散图 本文基于几何形态为曲面的连续介质有限变形理论，对海面上油污的扩散运动建立力学 模型，并按Lagrange观点进行数值求解，研究海面对海面上油污的扩散运动在几种不同情况 下的运动进行了模拟。 2 海面油污扩散的动力学模型 本部分主要基于当前物理构型对应之曲线坐标系显含时间的有限变形理论，具体建立海 面油污扩散的动力学模型。 2.1 连续性方程 连续性方程的积分形式为 0 d dt  d    其中 为面密度。根据第一类面输运可得其微分形式为     0 其中速度的散度      V V 。 2.2 动量方程 动量守恒方程可以表示为 int ext ten pre vis pre sur d Vd F dt   F F F F        式中， Ften 表示二维流体的表面张力， int F pre 表示二维流体的内压力， Fvis 表示二维流体的黏滞 力， ext F pre 表示二维流体的外压力， F sur 表示二位流体受到的面力

5二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 如果表面张力系数y视为常数,则动量守恒方程的微分形式可以表示为 P=H(-p)-Vp+vev+v8p)+hh(+4m+Jm 2.3 Lagrange观点下控制方程 在 Lagrange观点下,连续性方程的分量形式可以表示为 a(、n09+vr-H)=0 如果曲面可以表示为Mnge型(x,y)=y,则动量守恒方程的分量形式在只考虑 L=(x,y,) 重力的情况下可以表示为 川+x￥+2(8=,+82=m,)x+2(g=n乙+g2-0,)+g1==+82=n, =山Vy+Vv"+K-3b 2(Vb)V-bb,'V ap 8-Pg(8=+g=y b x [3VP+(Vb)y+vVP2-2H2-2。y]+H(-p)+4 2 3数值方法 3.1数值徽分 在上面的方程中,某量对坐标x或y的导数对于不同的量,含义是不一样的。如果该量是 曲面量(如gn、b、I等),那么其对坐标的导数即该时刻该曲面量对该时刻曲面坐标的导 数。如果该量是物理量(如V、p等),则其对坐标x或者y的导数需要做如下的理解 (,1)(x2,1)=(5 ax(,/)+g (

第二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 - 3 - 如果表面张力系数  视为常数，则动量守恒方程的微分形式可以表示为 ( ) [ ( ) ( )] sur dV H p n p V V Hn V pn f dt                    2.3 Lagrange 观点下控制方程 在 Lagrange 观点下，连续性方程的分量形式可以表示为 3 ( , , ) ) ( ( ) 0 i i i i x x x t t V H t x V                 如果曲面可以表示为 Monge 型 ( , , ) ( , , ) x x y t y z x y t             ，则动量守恒方程的分量形式在只考虑 重力的情况下可以表示为 1 1 2 1 1 2 1 3 3 1 2 2 2 [ ] 2( ) 2 2( ) ( z ) 3 ( ) g l i j l l l l l l ij tx x tx y ty x ty y tt x tt y l s s l l ls s l l s t s s G s s t s sl l l s x y g z z g z g z g z V x x x z x z x z z g z z V V V b V b V x p g g g z g z K b b x                                           1 2 2 2 2 2 2 3 2 3 2 2 2 1 1 1 3 ( ) 2( 2 ) ( ) 1 2 2 i j tx tt ty ij x y x y x y s t t s s t s t s s G x y z b x x x x z z z z z z g b V b V V H K V H p p z z z z                                                              3 数值方法 3.1 数值微分 在上面的方程中，某量对坐标 或 的导数对于不同的量，含义是不一样的。如果该量是 曲面量（如 ij g 、 ij b 、 l ij 等），那么其对坐标的导数即该时刻该曲面量对该时刻曲面坐标的导 数。如果该量是物理量（如 i V 、  等），则其对坐标 x 或者 y 的导数需要做如下的理解 1 1 2 2 ( ( , ), ) ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) A i A i i i f f f f x t t t x t t x t t x t x x x x                                             

5二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 式中,f=f(x(2,)1)表示该时刻的某物理量,即为初始构型坐标。其中的(x2,1) 可以看作是变换=5(x,D)对应的 Jacobi矩阵的元素,它可以由 202ax 得出。 类似的,该物理量的二阶导数根据链式求导法则要满足 a(x2(52,)=、0 B(52,1)(x2)(x,)+(5, (xs, 1) 藉此模拟中所有的量都可以表示为对初始构型的导数,可以采取等距 Lagrange插值的五点公 式 3.2数值积分 四阶 Adams预测校正方法是一种线性多步法,与作为单步法的 Runge- Kutta方法不同之 处在于,它在计算某节点ln1处的近似值y时用到了前面多个点l1,L1…处的信息 这种方法称为多步法。 四阶 Adams预测校正方法可以表示为 y=y+(551-591+372-9f3) yn=y+,(97m+19-5+f2) 其中∫=∫(1,y2),f1=∫(1)。第一式称为预测公式,第二式称为校正公式,并使用 经典四阶 Runge-Kuta公式作启动值计算。它的实现模式可以表示为 1)P步(预测):j1=y1+(55f-591+372-9-3) 2)E步(计算) f(Lvau h (3)C步(校正):1n=y+n4(9f+191-J1+f2)

第二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 - 4 - 式中， f f x t t  ( ( , ), )    表示该时刻的某物理量，   即为初始构型坐标。其中的 ( , ) A i x t x       可以看作是变换   ( , ) x t     对应的 Jacobi 矩阵的元素，它可以由 1 1 1 2 2 2 2 1 2 2 2 1 1 1 1 1 2 x x x x x x x x                                                                   得出。 类似的，该物理量的二阶导数根据链式求导法则要满足 2 2 2 (x ( , ), ) ( ( ) ( , ) ( , ) , , ) ( , ) A B A i j A B i j A i j f f f t t x x t t x t x x x x x t x t                                              藉此模拟中所有的量都可以表示为对初始构型的导数，可以采取等距 Lagrange 插值的五点公 式。 3.2 数值积分 四阶 Adams 预测校正方法是一种线性多步法，与作为单步法的 Runge-Kutta 方法不同之 处在于，它在计算某节点 i 1 t  处的近似值 i 1 y  时用到了前面多个点 1 , , i i t t  处的信息 1 , , i i y y  ，这种方法称为多步法[4]。 四阶 Adams 预测校正方法可以表示为 1 1 2 3 1 1 2 1 (55 59 37 9 ) 24 (9 19 5 4 ) 2 i i i i i i i i i i i i h y y f f f f h y y f f f f                        其中 ( , ) i i i f f t  y ， 1 1 1 ( , ) i i i f f t y     。第一式称为预测公式，第二式称为校正公式，并使用 经典四阶 Runge-Kutta 公式作启动值计算。它的实现模式可以表示为 (1) P 步（预测）： 1 1 2 3 (55 59 37 9 ) 24 i i i i i i h y y f f f f          (2) E 步（计算）： 1 1 1 ( , ) i i i f f t y     (3) C 步（校正）： 1 1 1 2 (9 19 ) 24 i i i i i i h y y f f f f         

二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 (4)E步(计算):f=f(141,y1),为下一步的预测做准备。 4数值模拟 本文数值模拟基于§2.3所述 Lagrange观点下控制方程,但尚未能考虑油面的表面张力 作用及其内压力修正,油面自身内摩擦力以及油面同油污之间的摩擦力 4.1海面行波运动工况 运动曲面方程:(x,y,1)=0.lsin-x-o,其中λ=1.0,O=4.0。点阵大小 64×64。初始位置:x方向[-1,1;y方向[-1,1。初始速度:均为零。考虑重力作用,忽 略黏度。 图1展示了此工况曲面在t=0.30时的主要的几何量参数(平均曲率H和量H3)的分 布。不同时刻的几何量参数分布相似,只有绝对数值之间的差别。图2-图6列出了曲面上油 污的主要物理量随时间的变化情况 可以看到,由于曲面在y方向具有平移对称性,因此油污的密度场、速率场和速度场都 保持了相应的平移对称性。另外可以注意到,油污的密度在曲面的波谷处较大,而在波峰处 较小;随着时间的推移,油污整体位置在向负x方向拓展,而并没有随行波行进方向前进的趋 势。由图8和图9之间的对比可以知道,油污将在t=060到t=062这段时间内,在曲面的波 峰处密度由零值变为负值。这表明油污将在该处断开,分为几部分的较小的带状油污。可以 预见最终大片的油污将会断裂成为若干块带状碎块,其数量与该区域内波谷的数目相同。 遗憾的是由于此次模拟中采用的方程(主要是质量守恒方程)中并没有描述目标流体断 开的情况,因此模拟只能进行到断开现象发生,而后的结果都将因为密度的初始值为负而导 致密度急剧发散,失去模拟的意义。 (a)平均曲率分布 (b)H3的分布 图1海面行波的几何量分布(t=0.30)

第二十五届全国水动力学研讨会暨第十一届全国水动力学学术会议 - 5 - (4) E 步（计算）： 1 1 1 ( , ) i i i f f t y     ，为下一步的预测做准备。 4 数值模拟 本文数值模拟基于§2.3 所述 Lagrange 观点下控制方程，但尚未能考虑油面的表面张力 作用及其内压力修正，油面自身内摩擦力以及油面同油污之间的摩擦力。 4.1 海面行波运动工况 运动曲面方程： 2 z x y t ( , , ) 0.1sin x t            ，其中  1.0，  4.0 。点阵大小： 6464 。初始位置： x 方向 [ 1,1]  ； y 方向 [ 1,1]  。初始速度：均为零。考虑重力作用，忽 略黏度。 图 1 展示了此工况曲面在 t  0.30 时的主要的几何量参数（平均曲率 H 和量 3 HV ）的分 布。不同时刻的几何量参数分布相似，只有绝对数值之间的差别。图 2-图 6 列出了曲面上油 污的主要物理量随时间的变化情况。 可以看到，由于曲面在 y 方向具有平移对称性，因此油污的密度场、速率场和速度场都 保持了相应的平移对称性。另外可以注意到，油污的密度在曲面的波谷处较大，而在波峰处 较小；随着时间的推移，油污整体位置在向负 x 方向拓展，而并没有随行波行进方向前进的趋 势。由图8和图9之间的对比可以知道，油污将在 t  0.60 到 t  0.62 这段时间内，在曲面的波 峰处密度由零值变为负值。这表明油污将在该处断开，分为几部分的较小的带状油污。可以 预见最终大片的油污将会断裂成为若干块带状碎块，其数量与该区域内波谷的数目相同。 遗憾的是由于此次模拟中采用的方程（主要是质量守恒方程）中并没有描述目标流体断 开的情况，因此模拟只能进行到断开现象发生，而后的结果都将因为密度的初始值为负而导 致密度急剧发散，失去模拟的意义。 (a)平均曲率分布 (b) 3 HV 的分布 图1 海面行波的几何量分布(t=0.30)