北极航道的开通加速了各国在极区的产业投入，不过强冰况环境可能给海工结构物的结构安全带来威胁[1]．美国船级社在冰区规范中指出：当进行冰固作用分析时不仅要考虑安全系数范围内的设计载荷，更要考虑意外过载情况可能导致的结构非线性响应[2]．为了更好地分析结构物变形和海冰破坏的内部耦合关联，建立一套能实现冰体与结构物实时动态耦合作用的数值预报方法具有十分重要的工程意义．目前，国内外学者开展了一系列冰体与结构物耦合数值预报方法研究．文献[3]基于拉格朗日流固耦合模型分析了液柱掉落悬臂冰梁动态过程；文献[4]利用离散元(discrete element method，DEM)方法建立了不同形状海冰和垂直立柱接触数值预报模型，但是无法进行立柱结构响应预报分析；文献[5]基于DEM方法构造了海冰离散单元模型；文献[6]利用连续光滑流体动力学(smoothed particle hydrodynamics，SPH)方法开展了船冰流耦合作用预报研究；文献[7]基于近场动力学(peridynamics，PD)方法开展了船艏与碎冰接触作用预报研究，分析了船冰作用过程冰载荷时变特征；文献[8]以Abaqus为二次开发平台，基于断裂力学和损伤力学理论提出了一种冰与结构物耦合预报模型；文献[9]基于黏聚单元法(cohesive element method，CEM)开展了冰桨接触作用预报，但是研究中没有进一步考虑螺旋桨的结构响应；文献[10]基于近场动力学和有限元方法耦合建立了冰桨接触预报模型，但是模型中忽略了结构变形本身对计算结果的影响；文献[11]利用光滑粒子流体动力学-有限元法(SPH-FEM)建立了新型抗冰结构响应模型，但模型在模拟冰体动力破坏上还不够完善．通过上述研究可以发现：针对冰固耦合问题，目前学者多关注在冰载荷和海冰破坏的预报，较少考虑强冰况下结构物可能产生的非线性响应，更是鲜少关注结构物非线性变形特征给冰固耦合的反馈影响．本研究通过自编程手段，利用常规状态型近场动力学方法建立海冰破坏模型，基于非线性动力学理论建立结构物响应预报模型，并提出一种交界面冰体物质点与结构物面元信息的交互传递方法，最终建立冰体与结构物非线性动力耦合预报模型．基于该模型开展了冰柱高速冲击悬臂梁案例分析，计算结果验证了数值方法的有效性．1 基础理论1.1　状态型近场动力学理论在近场动力学理论中，研究对象将离散成一系列包含材料密度、弹性模量等物质信息的物质点．在该计算域中，物质点只与其领域范围内物质点存在力的作用，为了更好地表示不同物质点之间的本构关系，将两个物质点之间的相对位置和相对位移矢量定义为ξ=x'-x，η=u(x',t)-u(x,t)，其中：x和x'分别为物质点变形前后的位置；u为t时刻该物质点的位移矢量[12]．近场动力学中物质点运动须满足质量守恒和动量守恒，而在常规态基理论中，变形后力密度矢量T|x|的大小不相等，方向与变形矢量方向Y|x|平行．即T|x|×Y|x|=0，其本构模型为T=t1M1，式中：M1为单位方向状态矢量，M1=(ξ+η)/|(ξ+η)|；t1为力标量状态．近场动力学最核心的方程是物质点的运动方程，该理论的公式推导也都是围绕推导与求解运动方程进行的，运动方程表示为ρu¨(x,t)=∫Hx[T(x,t)x'-x-T(x',t)x-x']dVx'+b(x,t),式中：ρ为材料密度；u¨为物质点加速度矢量；b为外力密度；Vx'和Hx分别为近场动力学邻域的体积和邻域半径．物质点对作用会随着材料变形而变化，通过伸长率s来表示这种变化，即s=(|ξ+η|-|ξ|)/|ξ|.当伸长率超过某一阀值时，物质点对发生断裂，两端作用点失去力作用，该阀值定义为临界伸长率s0，引入判定函数μ来判断点对是否发生断裂，定义为μ(x,ξ)=1    (ss0);0    (其他).1.2　非线性结构力学结构经过有限元法离散为若干单元后，结构动力学方程归纳为Mq¨t+Ctq˙t+Ktqt=Rt，(1)式中：M为总体质量矩阵，不随时间变化；Ct，Kt和Rt分别为t时刻结构总体的阻尼矩阵、刚度矩阵和外载荷；q¨t，q˙t和qt分别为节点在t时刻的加速度、速度和位移．本研究采用纽马克法．依据纽马克法的基本假设，在t至(t+Δt)时间段内，加速度、速度和位移的关系为q˙t+Δt=q˙t+[(1-α)q¨t+αq¨t+Δt]Δt；qt+Δt=qt+q˙tΔt+[(1/2-β)q¨t+βq¨t+Δt]Δt2，式中：Δt为时间增量步；α和β为不依赖频率的常数．经过推导，非线性结构动力平衡方程可归纳为K¯t+ΔtkΔΔqt+Δtk+1=R¯t+Δtk (2)式中：K¯t+Δtk为t+Δt时间增量步下第k次迭代的总刚矩阵；ΔΔqt+Δtk+1为t+Δt时间增量步下第k+1次迭代的位移增量；R¯t+Δtk为t+Δt时间增量步下第k次迭代的外载．K¯t+Δtk和R¯t+Δtk具体展开为：K¯t+Δtk=4Δt2M+2ΔtCt+Δtk+Kt+Δtk；R¯t+Δtk=M4Δtq˙t+q¨-4Δt2Δqt+Δtk+Ct+Δtkq˙t-2ΔtΔqt+Δtk-Kt+Δtkqt+Δtk+Rt+Δt,式中上标k为该时间增量步下的第k次迭代．2 基础模型验证2.1　非线性结构力学计算模型以悬臂梁受恒定载荷作用下的阻尼振动算例来验证本研究构建的非线性结构动力学模型的有效性．该悬臂梁左端为固定支座，右端为自由端，在梁上表面施加沿x轴负方向的均布载荷P=10 kN/m2．悬臂梁的长度L1=0.05 m，宽度B=0.01 m，厚度H=0.001 m．材料的弹性模量E=210 GPa，密度ρ=7 600 kg/m3，泊松比υ=0.33．首先针对该模型开展网格收敛性分析，选择一套适合的网格尺寸．沿着悬臂梁长度方向设置了三组不同尺度的网格尺度L1/20，L1/30和L1/40，其中L1为悬臂梁长度．图1给出了不同网格尺度下悬臂梁非线性振动预报结果，可以看出：悬臂梁最大位移发生在梁体的自由端处，这和悬臂梁的变形特性相吻合．图2给出了悬臂梁不同网格尺度下非线性振动位移时历曲线，其中U为加载方向上变形的位移．可以看出：悬臂梁在载荷作用下迅速发生最大程度的变形，然后在阻尼作用下往复振动，振幅迅速衰减，最后停留在平衡位置．当网格尺度为L1/30和L1/40时，无论是振动的频率还是最后稳定的平衡位置符合度均较好，表明当网格尺度小于L1/30时，预报结果基本收敛，因此为了保证计算精度的同时尽可能减少计算量，后面分析中选择的网格尺度为L1/30．10.13245/j.hust.210820.F001图1不同网格尺度下悬臂梁非线性振动预报结果10.13245/j.hust.210820.F002图2悬臂梁不同网格尺度下非线性振动位移时历曲线然后开展了不同载荷强度下悬臂梁线性和非线性结构动力学程序预报结果对比，用以说明在强载作用下采用非线性计算模型对保证计算结果准确性的重要影响．须要注意的是，强弱载荷是针对悬臂梁的结构响应提出的，弱载工况表明可以用线弹性、小变形假设来计算结构响应，而强载则须要考虑到结构响应中的非线性特性．案例中计算参数设置和上述研究一致，选择网格尺度L1/30．针对弱载工况，当t0.0 s时，在梁上表面施加p=-10 kN/m2的均布载荷；针对强载工况，当t0.0 s时，在梁上表面施加p=-100 kN/m2的均布载荷．图3和4分别给出了不同大小载荷作用下采用线性、非线性模型预报的时历位移曲线和变形云图．10.13245/j.hust.210820.F003图3悬臂梁线性、非线性时历位移曲线10.13245/j.hust.210820.F004图4悬臂梁线性、非线性变形云图图3可以看出：在10 kN/m2载荷作用下，非线性预报结果和线性预报结果符合度很高，这是因为当载荷较小时，结构发生的形变也较小．当载荷增大为100 kN/m2时，对于线性和非线性模型来说无论是振动曲线时变特征还是最后稳定收敛结果均存在较大差异，线性模型最后稳定收敛于0.022 m，而非线性模型收敛于0.019 m，结果误差15.7%，非线性计算稳态结果小于线性计算结果．非线性模型振动收敛速度也明显要快于线性模型的收敛速度，以上分析表明忽略结构非线性会减弱悬臂梁结构的刚度，与文献[13]中的结果相符合．2.2　海冰破坏近场动力学计算模型为了验证近场动力学方法能合理的模拟海冰破坏过程，本研究参考文献[14]开展的柱状冰体冲击试验，利用状态型近场动力学方法构建相应的数值计算模型．针对试验结果进行同工况下的数值预报对比计算．试验中柱状冰体直径D=0.117 m，长度L2=0.421 m，冲击速度V=60 m/s．近场动力学中邻域半径δ=3.015dx，泊松比υ=0.33，时间步长dt=3.6×10-6 s，弹性模量E=9.31 GPa，极限键伸长率s=4.9×10-4[15]．图5给出了不同网格尺度冰柱冰柱撞击时历载荷对比，图中f为冲击载荷．可以看出：网格尺度对冰柱冲击载荷的幅值影响不大，对载荷曲线特征特别是下降的斜率存在较大影响，预报和试验结果中载荷曲线都表现出先迅速增加到幅值后相对缓慢减小的特征，这是因为在初始接触时刻，冰柱撞击到刚性板上发生较快的弹性变形，载荷迅速增大．图5还反映出在对载荷的预报中无论是冲击幅值还是整体载荷曲线和试验结果的符合度较高，证明了本研究建立冰体破坏预报模型的准确性．10.13245/j.hust.210820.F005图5不同网格尺度下冰柱撞击时历载荷对比3 非线性耦合预报算法及算例分析3.1　非线性耦合预报算法基于前期建立的非线性结构动力预报方法和海冰破坏动力方法，提出了一种计算冰与悬臂梁结构物非线性动力耦合预报方法，给出了冰柱高速度撞击悬臂梁的计算算例．总载荷的计算公式为Fi=14∑finp|np|，式中：Fi为第i块面元的总载荷；fi为物质点冲击面元的载荷大小；np为面元的法向量．3.2　算例分析为了验证本研究建立的冰与悬臂梁结构物非线性动力耦合预报模型的准确性，以冰柱体冲击悬臂梁过程模拟来开展基准算例分析，悬臂梁的长度L1=0.112 m，宽度B=0.024 m，厚度H=0.004 m，材料的弹性模量E=210 GPa，密度ρ=7 600 kg/m3，泊松比υ=0.33．图6为冰柱撞击悬臂梁实时动态耦合过程，图中Q为冰柱的破坏剧烈程度．可以看出：预报模型结果能较为合理地反映冰柱和悬臂梁作用的过程．冰柱高速冲击作用到悬臂梁自由端附近，悬臂梁受冲击载荷作用逐渐发生变形．图6(c)中，当裂纹刚扩展到圆柱体半腰区域时，上端面区域率先发生破坏，是典型的层裂破坏[16]．同时对比图6(b)、(c)和(d)可以发现：图6(b)和(c)中，由于冰体刚和结构物接触作用不久，因此结构还能保证完整，冰柱基本是垂直于悬臂梁上表面的．图6(b)中悬臂梁上表面变形较小，冰柱基本呈竖直状态，而图6(c)中悬臂梁上表面有一定的变形，冰柱朝着悬臂梁自由端倾覆，图6(d)中冰柱又沿着悬臂梁固定端倾覆，这是因为该时刻下圆柱冰体结构不再完整，悬臂梁发生的形变导致圆柱左下端受到严重破坏，而右下端冰柱结构相对完整，因此在弯矩作用下圆柱体朝着悬臂梁固定端倾覆．10.13245/j.hust.210820.F006图6冰柱撞击悬臂梁实时动态耦合过程图7为冰柱撞击悬臂梁过程模拟，为了更好地10.13245/j.hust.210820.F007图7冰柱撞击悬臂梁过程模拟理解这一过程，图7(c)给出了整个过程的概念图．从图7可以看出：冰柱在冲击过程中先是冰柱底部紧贴着变形的悬臂梁上表面呈一定倾角竖立，柱体表面出现斜向上扩展的张开型裂纹．随着冲击过程的进行，裂纹沿着原方向进一步扩展并逐渐增宽，最后可以看到冰体底部结构会沿着裂纹扩展路径破坏裂解成碎冰块．当裂纹扩展到冰柱上半区域时就不再向上延伸，这是因为底部结构的完全破坏，冰柱主体结构几乎不再和壁面发生作用，裂纹失去了进一步扩展的冲击能量，裂纹停止扩展．4 结论本研究针对强冰况下结构物产生的非线性变形响应问题，建立了一种冰与悬臂梁结构物非线性动力耦合预报方法，可以为强冰况下极地海工结构物的响应预报和强度校核提供方法参考．通过研究主要得出以下结论．a. 建立了冰与悬臂梁结构物非线性动力耦合预报方法，该方法对于梁体变形反馈给柱体破坏的影响模拟的较为准确，预报方法具有可行性．b. 冰柱冲击悬臂梁过程中，悬臂梁的变形使得冰柱表面产生斜向上的裂纹，裂纹随着冲击进行逐渐扩展延伸并进一步拓宽，最后冰柱底部沿着裂纹扩展路径破坏裂解成碎冰块．c. 冰柱的姿态会随着梁的变形情况和柱体结构损伤情况的变化而改变，前期柱体贴合着变形的梁上表面向梁自由端倾覆，后期冰柱左下端区域结构瓦解，柱体向梁固定端倾覆．本研究主要工作是建立冰与结构物耦合预报方法，但目前只涉及了形状较为简单的悬臂梁和柱状冰，后续须开展更多复杂工况计算和试验验证工作．