孤立冰块流拖曳系数的优化辨识

孤立冰块流拖曳系数的优化辨识 孤立冰块流拖曳系数的优化辨识 第18卷第2期 2006年6月 极地研究 CHINESEJOURNALOFPOLARRESEARCH Vo1.18.No.2 June2006 ^—^—^,,,, ;研究{ ,,,^,,,,,,,, 孤立冰块流拖曳系数的优化辨识 白乙拉李志军卢鹏张丽敏 (大连理工大学应用数学系,大连116024; 大连理工大学海岸和近海工程国家重点实验室,大连116024) 提要流拖曳系数是海冰动力学模型中的重要参数,实际海冰运动时受到水流的底面剪切作用 和侧面正压力作用,相应的拖曳系数分别为摩拖曳系数和形拖曳系数.利用实验室波流水槽设 备,对长方平底,长方粗糙底和圆饼型淡水冰试样进行拖曳运动试验,获得26组冰样运动数据. 根据这些试验数据,利用动量法原理和数学优化方法构造了摩拖曳系数和形拖曳系数的辨识模 型,将摩拖曳系数和形拖曳系数两项分开.在此基础上,分析了冰样粗糙度指标(支承长度率和 冰厚均方根偏差)对这两个系数的影响;发现了摩拖曳系数同冰样水平尺寸和垂直尺寸之比的 关系,该关系能将冰山和浮冰的摩拖曳系数归结到同一条曲线上. 关键词浮冰拖曳系数粗糙度优化辨识 1引言 极地海冰系统是全球最重要的冷源,它直接影响全球尺度的大气环流,大洋环流和气 候变异,因而它成为全球气候变化的驱动器和响应器之一,近年来全球气候变化的加剧推 动了对极地海冰问题的研究.极地海冰数值模拟的研究从2O世纪6O年代开始,主要包 括:海冰生消漂移模型,海冰本构方程,海冰与大气一海洋的相互作用及海冰的动力破坏 等.其中海冰的动力学模型主要描述的是海冰与大气,海洋间的动力相互作用,模型中涉 及到众多的海冰形态参数,而这些参数的准确性直接影响到动力学模型的精度.拖曳系 数是动力学模型中的重要参数,包括风拖曳系数和流拖曳系数,它们直接影响风,流对冰 的拖曳力.所以前人已经在这方面进行了大量研究(WamserandMartinson,1993),目前 常用的拖曳系数确定方法可分为三类:涡动法,剖面法和动量法.涡动法是综合考虑大 气,海水和海冰的温度,速度等物理要素来计算风流的拖曳系数(GuestandDavidson, 1991);剖面法是通过冰上风速和冰下流速剖面来计算冰面粗糙长度,进而得到拖曳系数 (Andreas,1993);动量法则是根据海冰漂移的动力学方程来确定拖曳系数(Ma~inson, [收稿Et期] 【基金项目] [作者简介] 2005年11月收到来稿,2006年3月收到修改稿. 国家自然科学基金重点项目(40233032),科技部科技基础性工作和社会公益研究专项 (2003DEB5J057)和国家自然科学基金(40376006)资助. 白乙拉,男,1961年出生.在读博士,主要从事控制论与海冰数值模式研究. 第2期白乙拉等:孤立冰块流拖曳系数的优化辨识 1990).利用上述方法,已经获得一些海域海冰拖曳系数的经验值,并应用到了动力学模 型中(Hibler,1979;Wueta1.,1997).不过存在的问题是:海冰的形态变化对拖曳系数的 影响没有定量解释,而是给出一个定值,但这个定值的适应性仍缺少定论,所以目前仍然 有学者致力于寻找气一冰和流一冰界面粗糙度同拖曳系数的关系(PrinsenbergandPeterson, 2002;Steinereta1.,1999). 另外,冰山作为海冰的特殊形态,其拖曳系数差异很大,主要是由于水平和垂直尺度 的比值范围较大而引起的.如果冰山的拖曳系数也像大面积浮冰那样,用一个综合考虑 了摩阻力和形阻力影响的拖曳系数来表达就有缺陷.本工作就是以实验室流对冰的拖曳 试验数据为基础,再利用数学中的参数辨识方法,从中分离出摩拖曳系数和形拖曳系数, 从而探讨尺寸效应和冰层底面粗糙度效应.文中重点介绍试验和辨识方法,结果,简单讨 论所得拖曳系数同尺寸和粗糙度之间的关系. 2水槽试验和关键参数提取 试验在大连理工大学海岸和近海工程国家重点实验室的波流水槽中进行,水槽长 22m,宽0.45m,深0.6m.试验之前对造流系统进行标定,在40cm水深条件下,水面以下 5cm和15cm处的流速基本相同,本试验水面以下15cm之内的流速为6.75cm/s. 试验选用不同厚度和形状的淡水冰块,用CCD采集它们在波流水槽中被水流拖曳运 动过程,以此计算其速度和加速度,试验装置如图1所示. 试验时CCD距离水面70cm,有效拍摄宽度范围大约25cm,冰块最大入水深度13cm; CCD采集频率25Hz.每个冰块试验3次,均在有效拍摄区间内从静止开始采集图像.分 析时为了凸现每帧之间的差距,每隔4帧取一次图像,即分析用的时间间隔为0.2s.为了 标定和图像处理的方便,在冰块上表面都贴有边长已知的正三角形标志物,有2cm,5cm, 8cm和10cm4种.图像中每个像素代表0.68587×0.68587mm.图2是CCD采集的某 一 时刻的冰块图像. 图1试验装置 Fig.1.Experimentequipments 图2CCD采集冰在水槽中某一时刻的状态 Fig.2.ACCDviewoficefloeinthetank 试验中冰块在水槽内不发生旋转,从图像中直接获取冰块上标志三角形最前面顶点 l10极地研究第18卷 的像素坐标位置,用来代表该时刻冰块的位置,所以在计算位移时也仅考虑沿水流方向的 位移.由冰块标志点的位移除以时间间隔即为冰块在该加速时段的平均速度. 由于冰块基本是匀加速直线运动,其位移随时问的变化关系为 |s=口+(1) 根据冰块的位移/时间关系拟合出一个二次曲线,其二次项系数就是加速度的一半, 这样求出的加速度能较好地代表整个加速时段的平均值,如图3所示.拟合得到的二次 曲线的一次项系数是初速度,表现了在释放冰块时其他因素的干扰,并不能使冰块达到理 想的"静止"状态,但它对于提取加速度不会造成影响. 5 4 旨3 2 l 0 024T/s68l0 图3某次试验加速阶段的位移/时间拟合曲线 Fig.3.Movementwithuniformacceleration 3流拖曳系数的优化辨识方法和结果 水流对冰块的拖曳力分为摩阻力和形阻力两项,可表示为 .1 = Fl+F2=Ma=pClSl(一)+c2|s2(一)(2)二 其中:F为水流拖曳力,为冰块质量,口为冰块运动加速度,F,,F为拖曳力中的摩阻力 项和形阻力项,p为水的密度,|s为冰块下底面面积,|s为冰块水下部分在垂直水流方向 上的投影面积(侧面面积),为水流速度,为冰块运动速度.公式(2)中的,,p |s,|s是在试验中已知的测试数据,口及由试验数据提取,c与c是待辨识的摩拖曳 系数和形拖曳系数. 设每个冰块重复试验?次,根据每次采集的数据,利用前述方法计算出本次试验中 冰块运动的平均速度和加速度口i,iE,:={1,2,…,?}.对于每个试验组次的?个样 本,理论上每个样本中的摩拖曳系数c和形拖曳系数c是不变的,目的是要确定这两个 参数.为此,对每个试验组次的?个样本,要在c及c的允许范围内,使式(2)左端项 和右端项之差的平方和达到最小,这是一个带约束集条件的函数求极值问题.将上述思 想用数学语言描述如下: 考虑摩拖曳系数c和形拖曳系数c的实际物理意义,参数c和c必有上下界,令 第2期白乙拉等:孤立冰块流拖曳系数的优化辨识111 :={(c,c)lCC?C,Cz?C?C,其中c…c,cz,c为确定的实数}CR 为有界闭集,定义U:={(c,c)I(c,c)?WCR}为辨识参数(c,c)的允许集. 对任意(c,c)?U,定义目标函数为 ,v1 -,(C,C2)=?[M?口一(c.Js+寺cz.Jsz).(一.).p](3)l=l一 这样以-,(c,c)为性能指标,辨识参数(c,c)?Ua的优化辨识模型记为 OPT:minJ(C1,C2) s.t(C1,C2)?Uad 由于对任意(c,c)?,J(C,C)是上的连续可微函数,而且是中的 有界闭集,因此优化辨识模型OPT的最优解(c,c)?U存在且唯一.对此辨识参数 (c,c)的非线性优化问题,本文采用遗传算法求解,其主要运算步骤如下: Step1:设置进化代计数器,;设置最大进化代数To;随机生成L个个体q(c,c) 作为初始群体Q(o). Step2:对每一个体,利用目标函数-,(c,c)计算群体Q(t)中各个个体q.的适应 度. Step3:按轮盘赌选择算子进行选择. Step4:依交叉概率P.?[0.5,0.9],对浮点数编码表示的个体进行算术交叉,交叉系 数从(0,1)区间随机选取. Step5:依变异概率P?[0.001,0.08]进行变异,变异时与进化代数有关.群体Q (t)经过选择,交叉,变异运算之后得到下一代群体P(t+1). Step6:若,?To,则t+--(t+1),转到Step2;若t>To,则以进化过程中所得到的具有最 大适应度的个体作为最优解,终止计算. 每个冰块做三次试验,可以得到三组口和值,也就是每个冰块有三个样本(即N: 3).对于每个冰块的三个样本,为使(2)式左右两端项之差的平方和最小,在允许集 内,按遗传算法的规则依次调整c和c,能使(3)式达到最小值的c,c即为所求.实 际计算时,根据前人的结果(岳前进等,1999;季顺迎等,2003),可以给出参数c和c 的大致取值范围为1.0X10一---9.9X10,,最大进化代数取值为1000,种群个体数 取值为8O. 本试验累积3O组次,通过质量控制得到26组可用数据,其中平底长方体冰块l3个, 带底纹长方体冰块9个,圆饼型冰块4个.经辨识计算得到上述26组有效数据所对应的 摩拖曳系数和形拖曳系数值,试验条件和辨识结果归纳为表1. 4结果分析与讨论 为了寻求拖曳系数同试样尺寸和粗糙度的关系,试验使用不同形状,不同厚度和不同 底纹的冰样进行试验.为此,选择了能够描述表面形态,又能方便的从现场观测中得到, 并同拖曳力有直接关系的两个粗糙度指标.第一个为冰底面的支承长度率R(c),其物 l12极地研究第18卷 理意义是在给定水平位置C上轮廓的实体长度Ml(C)与评定长度z的比率.平底冰 试样的Rmr(c)为1,粗糙底面冰试样的(c)均大于1;第二个为冰厚的均方根偏差R. 表1试验条件和辨识结果 Table1.Resultsofidentification 4.1支承长度率效应 支承长度率越大,表面起伏越大.从摩擦力的物理意义看,表面起伏越大,摩拖曳系 数越大.将辨识的摩拖曳系数(C)和粗糙度指标支承长度率尺mr(c)关联,得到图4.从 图4可以发现摩拖曳系数随冰底面的支承长度率增加而增加.但由于冰样水平尺度和垂 直尺度比值的效应,数据有一定的分散性,但总的趋势是存在的. 007 006 005 004 003 002 001 000 图4摩拖曳系数随粗糙度变化的关系 Fig.4.Skindragcoefficientvs.surfaceroughheSS 第2期白乙拉等:孤立冰块流拖曳系数的优化辨识ll3 4.2厚度均方差效应 形拖曳系数是冰侧面受水流作用力的体现,所以形拖曳力直接与水下冰厚度有关,冰 底面起伏在其侧面上的投影也引起水下厚度的变化,所以评价形拖曳系数,选择水下平均 冰厚度和均方差之和作为评价指标(尺),该指标在现场也容易从实测冰厚剖面资料中获 得.对于实验室试验,R=0.917×H+尺,H为平均冰厚,尺为实测冰厚的均方差,是粗 糙度的指标之一.形拖曳系数同该指标的试验关系见图5,它反映形拖曳系数同评价指 标尺之间存在非线性关系.由于尺是有量纲参数,使用时注意实验室物理模拟的比尺 关系. 4.3试样水平和垂直尺寸效应 对于海上存在的大面积浮冰,由于其水平尺度远远大于冰厚度,即使其形拖曳系数很 大,但形拖曳力也相对摩拖曳力小许多个量级.对于冰山,情况就比较复杂,这时同其水 平尺寸和水下尺寸的比例有直接关系.为此,我们引入一个新指标,即冰样长度与水下厚 度的比值,用它来寻求尺M(c)=1时摩拖曳系数的分布规律,结果见图6. 0 0 0 0 0 0 0.06 0.05 0.04 一 ru0.03 0.02 0.0l 0 图5形拖曳系数与水下厚度的关系图6摩拖曳系数同冰样长度与水下厚度比值的关系 Fig.5.FormdragcoefficientV8.drownFig.6.SkindragcoefficientV8.1ength/drown depthoficefloedepthoficefloe 由图可见,摩拖曳系数随冰样长度与水下厚度比值的增大而减小,图中曲线能够统一 冰山和浮冰摩拖曳系数的取值问题.对于海冰,从小浮冰到冰原对应长度/水下厚度比值 可以在5—1oo00之间变化,根据图6的长方平底拟合曲线[C=0.2352(L/R)'洲,r = 0.908],海上浮冰的摩拖曳系数相对很小,在1.01×10,一12.1×10之间,目前一般 海冰动力学模型取5.5×10,,就在上述范围之内(Hibler,1986).从图中还可看出,粗 糙底面比平底面的摩拖曳系数略高一些,但差异不大,具体可以通过图4解决.对冰山而 言,其长度/水下厚度比值在1.5—100的范围,所以冰山的摩拖曳系数相对浮冰而言,有 较大变化范围4.87×10,--0.1114.它们既可以在10.量级,也可以在10量级,这取 决于冰山的水平尺度和水下厚度的比值.目前在冰山动力学模型方面,有些文献将冰山 的摩拖曳系数设为0.002,然后根据实测资料拟合形拖曳系数,这样就忽视了冰山尺度效 应.仍然产生不同冰山观测结果,产生不同形拖曳系数.有些文献因冰山的摩阻力相对 形阻力较小,摩阻力项被忽略.从本试验的结果看,当冰山较小时,这种忽略会产生错误. 过去因为受观测能力的限制,获得的冰山运动数据均针对大冰山,所以当时的处理是正确 l14极地研究第18卷 的.现在利用船舶雷达跟踪到南极小冰山的运动(LiandLu,2005),其尺寸量级就能用 到图6结果.将图6的结果发展成为包含两个粗糙度指标的统计关系是必要的,但还需 更多的试验数据和不同条件验证,这是未来工作的方向. 5讨论 在实验室水槽中进行了水流对长方平底,长方粗底及圆饼三种不同形状淡水冰试样 的拖曳试验,试验获得的冰块试样运动过程数据,体现了流对试样摩拖曳力和形拖曳力的 联合作用.这种联合作用类似海洋原型情况,但问题比现场原型简单,即缺少风的作用. 利用数学中的优化辨识方法,将摩拖曳力和形拖曳力区分开,进而得到冰块运动的摩拖曳 系数和形拖曳系数.根据材料表面粗糙度指标的概念,引入两个具有物理意义的粗糙度 指标,得到流拖曳系数同冰粗糙度指标的关系.在此基础上,还得到摩拖曳系数同冰样水 平尺度和垂直尺度之比的关系.该关系可以统一解释从冰山到小浮冰块,大浮冰原的摩 拖曳系数相互之间的关系和取值问题. 参考文献 岳前进,季顺迎,毕祥军(1999):渤海海冰数值模式中计算参数的确定,海洋预报,16(3),97—103. 季顺迎,王瑞学,毕祥军,岳前进(2003):海冰拖曳系数的确定方法研究,冰川冻土,25(2)增刊,299--303. AndreasEL(1993):RatlshnessofWeddleSeaiceandestimateofair—icedragcoefficient,JournalofGeophysicalResearch, 98(C7),12439--12452. GuestPSandDavidsonKL(1991):TheaerodynamicrouIghnessofdifferenttypesofseaice,J ournalofGeophysicalResearch, 96(C3),4709—_4721. HiblerWD(1979):Adynamic— thermodynamicseaicemodel,JournalofOceanography,9,815—846. HiblerWD(1986):Icedynamics,In:TheGeophysicsofSeaIce,Ed.byUntemteinerN,PlenumPress,NewYork(NATO ASIB146).577一4O. UzJandLuP(2005):MeasurementofIcebergDraftwithMarineRadar,6thInternationalSymposiumonTestandMe~ure— ment,Dalian,China,2972--2975. Ma~insonDG(1990):IcedriftandmomenmmexchangeinwinterAntarcticpackice,JournalofGeophysicalResearch,94 (C5),1741—1755. PrinsenbergSandPetemonIK(2002):Variationsinair— icedragcoefficientduetoicesurface1"011ghhess,InternationalJo111- nalofOffshoreandPolarEngineering,12(2),1--5. SteinerN.Ha~erMandLemkeP(1999):Sea— icelOrdlg~xneBBanddragcoefficientsinadynamic—thermodynamic$ea—icemodel fortheArctic.Tellus.51A.964—_978. WamserCandMartinsonDG(1993):DragcoefficientforwinterAntarcticpackice,JournalofGeophysicalResearch.98 (C7),12432—12437. WuHDeta1.(1997):Numeficaisimulationfordynamicalprocessofseaice,ActaOceanologicaSinica,16(2),3O3—325. 第2期白乙拉等:孤立冰块流拖曳系数的优化辨识ll5 OPTIMIZATIONANDIDENTlFlCATlONONCURRENT DRAGCOEFFICIENTOFISOLATEDlCEFLOE BaiYila',LiZhijun,LuPeng2andZhangLimin (DepartmentofAppliedMathematics,DalianUniversityofTechnology,Dalian116024,China; StateKeyLab.rat.l了 ofCoastalandOffshoreEngineering,DalianUniversityofTeehn.1.gy,Dalian116024,China ) Abstract Currentdragcoefficientisoneofthemostimportantparametersinseaicedynamicsmod. e1.Inactualmotion,bothshearstressontheundersurfaceandpressureOilprofileareexerted onice,correspondingtoskindragcoefficientandformdragcoefficient, respectively.Byusing awave.currenttankinlaboratory,dragmovementexperimentswereperformedforfreshwater icesampleswithdifferentshapes,includingcuboidwithsmoothunderside.cuboidwithrough undersideanddiskwithsmoothunderside, and26datasetswereobtained.Basedontheseda. ta,anoptimizationmodeltodistinguishskindragcoefficientandformdragcoefficientwases- tablished.Thentheeffectoficeundersurfaceroughnessonthesetwodragcoefficientswasana. 1yzed,theresultrevealedtherelationshipbetweenskindragcoefficientandhorizontal/,vertical sizeoficesamples,whichcanmaketheskindragcoefficientsoffloeandiceberginasame CUn,e. KeywordsFloe,dragcoefficient,roughness,optimization,identificadon.