基于无网格粒子法的犁体挖掘土壤仿真模型

基于无网格粒子法的犁体挖掘土壤仿真模型

作为一种有效的填补和保护海洋石油管道的运行设备，该装置在海洋气产业中广泛用于传输。根据该设备的工作方法，其重力是衡量设备性能和运营成本的重要指标。因此，对牵引车身的力学过程有要求，可以有效评估设备的可靠性。根据海参崴式钻床的运行环境，有必要建立海参崴式钻床、海流模拟模型等。由于海底犁式挖沟机需要对海底土壤进行有效剪切,因此建立水流作用的海底土壤的仿真模型至关重要.在土力学领域,处理大变形的主要方法有有限元法和离散元法.有限元法处理大变形时,畸变的网格会导致计算不准确或者雅各比矩阵产生负值从而计算失败.而离散元只能计算一定数量的粒子仿真.由于有限元法和离散元法处理大变形问题的限制,采用完全的拉格朗日的无网格数值计算的方法(SPH)应用到该领域.该方法首先是由Gingold等提出,用于解决天文问题.之后,这种方法成功应用于其他计算物理力学的领域,比如材料强度的动态响应,大变形的流体动力学问题和海底爆炸问题.宿崇等采用FEM与SPH耦合算法,模拟了金属切削加工中材料的变形与应力情况.张阿曼等采用修正的SPH法成功的模拟水下冲击波导致的结构破坏的过程.在SPH方法处理土壤大变形问题上,杨望等采用修正库伦摩尔屈服准则,利用ANSYS/LS-DYNA求解器仿真土壤直剪试验.JiangZhong等采用LS-DYNA软件进行旋转犁片模型的实验.丁峻宏等比较自适应网格划分、ALE和SPH解决土壤切削仿真的研究,但是没有具体分析SPH建模方法.本文基于土壤的Drucker-Prager(D-P)屈服准则,采用无网格粒子法,将连续体粒子的差分方程转换成运动方程,再由拉格朗日数值方法进行求解.由于此方法不存在网格畸变问题,可以有效描述土壤的非线性大变形的过程而不需要额外的处理.本文主要研究基于一阶线性水波作用的海底土壤弹塑性模型,建立犁体,海底土壤及海流三者的相互作用模型,通过陆地实验验证ABAQUS动力学仿真软件求解无网格粒子法的可行性,模拟了犁体作用下土壤的运动过程,分析了不同作业环境下水流及海底土壤对犁体牵引力的影响.1曲面后翻抛变形海底犁式挖沟机的犁体主要用来将海底管道下方的土壤切削及翻抛.挖沟机的牵引力主要抵消土壤挖掘力.如图1所示,海底犁式挖沟机犁体部分主要包括主切削面和翻抛曲面.主切削面通过入土角,碎土角及推土角的同时作用挖掘土壤并引导其上升到翻抛曲面后翻转覆盖前一行程.由于土壤的剪切强度是由颗粒间的作用力决定,因此当土壤切削力达到极限平衡状态时土壤开始裂解.在主切削面切削土壤时,土壤可近似看成楔形体,犁体主切削面受到摩擦力及粘聚力、土壤重力、土壤支撑反力的作用.主切削面裂解的土壤以初始速度沿翻抛曲面的土迹线上升,设土壤沿土迹线的局部坐标系为uvw,当土壤上升速度达到0时,土壤则堆积到两侧,则犁体的受力方程:Μ=Ρ+μmg+Q-1(Ft+Fcr+Gu-μ(Fn-Gw)).(1)式中:P为犁体主切削面的切削土壤受力;m为挖沟机质量,μ为土壤与钢板的摩擦系数;Q为uvw转换成全局坐标系xyz的变换矩阵;Ft为翻抛面土壤的切向力,Fn为法向力;Gu,Gw为翻抛面土壤重力在局部坐标系uvw的分量;Fcr为水流产生的阻力.2土壤-地球计数2.1海底模型2.1.1abaqus/opplicit中屈服准则本文的土壤材料采用弹塑性模型,应用弹性理论的弹性应变和塑性增量理论计算的塑性应变.通常假设土壤为理想线性材料,近似成为线弹性且各向同性材料,因此可以采用胡克定律的弹性模量和泊松比描述土壤的本构关系.D-P模型可以用来定义如土壤等非金属摩擦材料的塑性行为.当压力增大,土壤强度变大.这类模型的屈服准则的屈服面基于J2-p子午面.屈服面包括线性函数,双曲函数和指数函数.采用ABAQUS/EXPLICIT进行求解计算,选择线性D-P模型,其屈服准则为f(p,√J2)=√J2-ptanβ-d=0.(2)式中:偏应力J2为√J2=√0.5sijsji=0.5q[1+k-1s-(1-k-1s)(r/q)3]‚s为偏应力张量;主应力为p=σ11+σ22+σ33,σ为应力张量;ks为三轴拉伸强度与三轴压缩强度之比,反映了中主应力对屈服强度的影响,一般取0.778;d为在J2-p应力空间的截距,根据抗压强度σc可以确定d=(1-1/3tanβ)σc.参数β可以通过三轴压缩和拉伸实验得出的摩尔库仑模型的摩擦角φ进行转换,tanβ=6sinφ/(3-sinφ).土壤的弹性阶段和塑性阶段连续,假设土地总应变及其增量分为可恢复的弹性应变和不可恢复的塑性应变,则其弹塑性模量矩阵的一般表达式为dσij=Dijdεij-Dijdλ∂g/∂σij.(3)式中:D为弹性模量矩阵.假设土壤适应流动性规则,塑性势面函数g与屈服面函数f重合,则满足材料的稳定性,弹性模量矩阵为对称阵.按照相关流动法则∂f∂σij=-13tan(βδij)+sij2√J2.(4)dλ可以表示成dε的函数:dλ=(-Κδijtanβ+Gsij2√J2)dεijΚtan2β+G.(5)式中:K为体积变形模量,G为切变模量.2.1.2海底土壤温度与波浪传播关系海底土壤主要受到水波产生的动压力的影响.设海水自由表面方程为yw=η(xw,t).根据一阶线性波理论,如图1所示,设水深为h,波浪为ω=2π的/T正弦波.可以求得一阶近似的速度场为˙xw=Aωchk(y+h)shkhcos(kx-ωt),(6)˙yw=Aωshk(y+h)shkhsin(kx-ωt).(7)根据浅水波概念,水深小于半个波长,因此浅水波动压力分布为Ρwave=ρwgΗcos(kx-ωt)/2ch(khw).(8)式中:A为波幅,ρw为水密度;根据色散关系,波浪的频率与波数之间的关系k=2π/L,L为波长,k为波浪的波数.假设海底土壤的孔隙水的压力与海底表面的压力相同,不考虑土壤的粘度,土壤应变产生的膨胀,流体压力导致的体积膨胀及内部颗粒间产生的接触力,则土壤平衡方程为σij,y+ρby=ρ¨yw.(9)式中:ρ=(1-ns)ρs+npw,ns为孔隙比,ρs为土壤密度,ρw为流体密度;σij,y为土壤沿着y方向受到的应力,¨y为流体在y方向的加速度,by为流体在y方向的体力加速度.2.1.3无网格粒子法求解应力无网格粒子法是基于流体粒子的想法,将计算域离散成有限数量粒子,通过采用与邻近粒子间插值的计算方法,将控制方程的积分形式转化成邻近离散粒子求和的形式,可以得到粒子的密度、速度、应力等物理量,因此不需要网格计算.该方法的控制方程包括粒子a的位置ra,固有质量ma,相关的密度ρa,压力pa,速度向量va,动态粘度μa,及热力学物理量如温度或者热量.其核函数为wh,则场函数方程A(r)为A(r)=∫ΩA(r′)wh(r-r′)dr′+o(h2).(10)式中:Ω为积分域,h为光滑粒子长度.通过应用Monte-Carlo方程,将上式离散成A(r)=∑bmbρb(Ab)wh(r-rb)+o(h2).(11)式中:b是这个域内每个粒子的下标,无穷小体积dr′代替粒子b的体积mb/ρb.最后,忽略高阶o(h2),并且考虑球形核,wh只取决于粒子对的距离,则A可以写成粒子a的位置形式.在标量场下,核函数可微,粒子a的梯度可以表示为∇A(ra)=∑bmbρbAb∇awh(rab).(12)式中:ᐁa为核函数的梯度.上式可以应用到不同向量场的离散,其中核函数通常采用三阶样条曲线,因此可以近似推导出粒子a的速度梯度为∂vai∂rj=∑b=1Νmbρb(vbi-vai)∂wab∂raj.(13)采用无网格粒子法得到的粒子a的各项物理量,将式(13)代入到式(5),则可得无网格粒子法的应力应变方程为dσaijdt=Dij-Dij(-Κδijtanβ+Gsij/2J2).(-tanβδij/3+sij/2J2)(Κtan2β+G)-1dεij.(14)式中:dεij=12(∑b=1Νmbρb(vbi-vai)∂wab∂raj+∑b=1Νmbρb(vbj-vaj)∂wab∂raj).2.2并建立农村气候模型如图1所示,在海底,只考虑x和y方向的水流作用.在不可压缩流体中,在x方向,流体受到的应力为流体在犁体表面产生的总阻力Fcr及海流沿着x方向的流速为σw,x-CDρwu022ds=ρww.(15)式中:CD为阻力系数,u0为海流在x方向的速度,ds为作用面积,w为水波在x方向产生的加速度.根据上述建立的土壤应力应变方程及水流对犁体表面的作用力后,在直角坐标系中建立犁体在海底的动力学方程:σp+σs+[σw,x00]+[0ρpbp0]=ρpU¨p.(16)式中:ρp为犁体密度,bp为犁体体力加速度,U¨p为犁体运动加速度,σp为犁体所受应力,σs为土壤作用在犁体的应力,σw,x为水作用在犁体的应力,忽略流体对犁体沿y方向作用力.2.3y=-hw-d为了分析整体运动过程,需要对各部分加入边界条件.首先,假设海底表面为非渗透性,则当y=-hw-d时,位移usx=usz=0;当x=0和x=xe时,usx=0;当z=0和z=ze时,usz=0.假设水波产生的压力只作用在海底表面,没有流体渗透,则犁体表面受到流体的切向压力为0,即∂p∂τ=0.3模拟模型和验证3.1模型的模拟3.1.1土壤条件模型的建立ABAQUS/CAE中有一系列连续的土壤模型,包括经典塑性模型,Drucker-Prager模型,双曲线模型,摩尔库伦失效模型和复杂模型.这些模型分别应用在各种特殊土壤类型.本文采用线性D-P模型,它经常应用在建立具有弹塑性行为的土壤.在分析过程中,土壤的剪切强度和摩擦角是常数,没有土壤应变硬化.在ABAQUS中建立长10m,宽5m,高5m的土壤层,建立PC3D单元,单元长度为0.2,光滑粒子长度为0.1.实验陆地场地和海底土壤的土工参数如表1所示.3.1.2模型的建立与组合海底犁式挖沟机的犁体结构根据沟型设计,采用水平直元线法,通过设定要求的入土角、碎土角、推土角,可以确定主切削面.再根据设计要求确定导曲线面、导曲线的开度、角度、高度,采用半螺旋线法确定水平直元线角的角度.在SOLIDWORKS建立犁体的三维模型,通过与ABAQUS的接口导入犁体三维离散刚体模型文件.犁体定义为壳体单元,单元类型为R3D4.根据结构钢的材料参数定义材料属性.由于海底犁式挖沟机为对称结构,取右侧犁体分析,根据挖沟作业的实际情况,将犁体模型与土壤模型装配,如图2所示,设定上述的边界条件.本文的接触条件采用面面接触,定义土壤与犁体的摩擦系数为0.5,土壤与犁体采用TIE约束关系.3.1.3海流流速仿真模型海底犁式挖沟机作业水深为20m,因此,根据文献,确定水波参数,如表2所示.结合海水水波的动态压力,对土壤与海流接触表面施加压力.由于水波压力为函数形式,将函数转化成位置函数与时间函数的乘积,分别采用分析场及增益函数定义.海流流速对犁体的作用采用不可压缩流体作用在三维绕流物体上的阻力.如图3所示,为某点在不同波高下的水波的动压力随时间变化曲线.为了进一步仿真海水与犁体的相互作用,可以采用欧拉体建立海水模型.由于欧拉体与PC3D单元不能建立相互作用,因此只能仿真海水对犁体的作用,不能分析海水对泥沙的冲刷作用.3.2不同人工粘度比例因子的确定建立上述仿真模型,采用ABAQUS显示动力学分析软件的无网格粒子法进行迭代计算,通过后处理,可得犁体随着入土到工作稳定沟深时的牵引力大小.在数值仿真时,经常会出现大变形的非物理性振动,导致不稳定.这是在初始步发生粒子冲击,粒子就会极不稳定.为了增加数值的稳定性,在动量方程中需要添加人工粘度,一般在INP文件中加入线性和二阶粘度,可以抑制由计算机响应产生的高频率噪声.人工粘度的表达式为Π=f(αs,cs,βs).(16)式中:cs为土壤的波速,一般为450～600m/s,αs和βs分别为线性和二阶人工粘度的比例因子.对于固体粒子仿真,Libersky等建议将2个参数设定为2.5,Ha等建议2个参数取0.1.为了确定合适的比例因子,分别计算不同人工粘度下的犁体牵引力的对照试验,由图4可以看出,当粒子速度小于声速,αs一定时,βs对仿真影响不大.当βs一定时,αs增大则峰值近似,结果更趋于稳定.因此本文仿真人工粘度比例因子均取2.5.为了验证仿真模型的可靠性,进行海底犁式挖沟机陆地实验.海底犁式挖沟机试验样机的陆地实验包括提供45t牵引力的绞车,恒定运行速度0.06m/s.绞车由在3m沟深中的2t重的钢管固定,陆地实验实况如图5所示.实验场地土工参数如表1所示.挖沟机牵引力由拉力传感器测量,由于陆地实验存在7°向下的牵引角度,通过计算可以对陆地实验数据及犁体牵引力的仿真实验结果进行对比,如图6所示.将试验数据与仿真数据进行3次样条曲线拟合,可以看出仿真曲线与实验曲线的趋势相同,仿真曲线与试验曲线的最大震荡误差为15%.仿真误差在一定程度上提高确定拖船牵引的安全性,满足允许的误差,因此可以证明采用无网格粒子法模拟土壤的犁体动力学仿真的可靠性,为模拟海底环境提供可行办法.仿真曲线出现震荡是由于仿真过程中冲击产生的粒子振动,为了进一步消除震荡,在今后的研究中可引入人工应力方法.4讨论通过与陆地试验结果的比较,验证了土壤仿真模型的可靠性.基于以上结果,可以进一步讨论在海底情况下,土壤及海流对犁体作用力的影响.4.1波高变化对冷却系统的影响为了验证海流对犁体牵引力作用的影响,在陆地模型的基础上,修改为海底土工参数,并加入水波对海底表面的动压力及海流对犁体的绕流阻力.本文只讨论浅水波对犁体牵引力的作用.通过仿真模型的结果可以得出水波对犁体牵引力有波动影响,海流的流速增大,对牵引力产生线性影响,如图7所示.图8中根据计算水波作用的牵引力与无波时的牵引力的差值比,可以得出波高变化对犁体牵引力的影响.当波高为2m时,震荡呈正弦变化,初始阶段波动幅值较大,随着深度的增加,差值比的震荡响应逐渐减小且朝正向发展.当波高增大时,初始阶段波高变化对牵引力震荡影响较大,这是由于波高对海底表面的动压力的增加.随着挖掘土壤深度的增加,土壤的裂解堆积在犁体表面的土壤增多,切削土壤导致粒子振动,从而削弱了由于水波产生的牵引力波动.4.2土地应力分析土工参数对于犁体牵引力的大小起到决定性作用,因此假设水波及海流流速一定时,讨论不同土工参数对其影响作用具有重要意义.表3为不同类型土壤的土工参数.针对表3中每种情况进行仿真,仿真结果如图9所示.由结果可以看出,随着沟深的增加,牵引力呈非线性上升趋势.提高每项土工参数,均会提高牵引力.在初始阶段