横向变形

横向变形（共3篇）

横向变形 篇1

在现代化冷轧宽带钢连续退火生产或热镀锌退火生产过程中,带钢经常会出现图1所示的横向瓢曲行为,其主要特征为带钢沿纵向出现一条或若干条“隆起”,即薄板在部分狭窄区域宽度上发生集体不等量也不同步的弹性或弹塑性面外位移.由于此类屈曲现象常常出现在连退、镀锌等生产线上的退火炉或锌锅高温环境中,因此,现场又称之为“热瓢曲”[1].

与板带材生产线上常见的浪形(如中浪、边浪等)和倾斜浪形[2]宏观表现明显不同,其特征是在带钢纵向出现一条甚至数条“隆起”,而且此种瓢曲现象波长从几十厘米到几十米不等;波宽方向则表现为很窄,只有几厘米到十几厘米宽,此类屈曲一旦产生都无法消除.轻微的横向瓢曲缺陷会造成废次产品数量增加,而严重的瓢曲则造成炉内断带事故并迫使机组停机,对机组生产及产品质量危害极大.但由于带钢瓢曲机理复杂,多年来对其成因一直未能得到理性的、正确的解释,成为严重困扰生产、亟待解决的技术难题.

传统认知认为横向温差、辊凸度等是导致“热瓢曲”的根源[3,4,5,6,7,8,9,10,11,12,13],本文通过分析认为,退火炉内为防止带钢跑偏而设置的锥形导向辊导致了带钢宽向局部张力集中,并诱发纵向不均匀拉伸变形而产生横向压应力,最终造成了带钢局部宽度上的“隆起”现象才是横向瓢曲产生的力学根源.为此,本文以某大型冷轧宽带钢连续退火生产线为研究对象,对带钢瓢曲现象的形成机理、影响因素与改进措施进行深入的研究,并通过解析分析揭示深层次的力学临界条件和应力场的变化规律,为退火炉内带钢工艺张力的设定提供参考依据.

1 横向瓢曲的有限元仿真

考虑到带钢材质各向力学性能并无显著差异,因而采用材料各向同性硬化屈服准则.同时,结合立式退火炉内导向辊的布置特点:带钢绕上下导向辊蛇形前行,相邻上下导向辊及其导辊间带钢就构成了炉内最小的一个结构单元,也就是带钢运行的一个最小周期,多个这样的重复结构单元(如图2)则构成了整个退火炉的内部结构框架.

本文采用了通用有限元软件进行建模仿真,如图3所示.边界约束条件:考虑到对称性,取1/2板宽作为研究对象,板宽中点作对称约束;带钢两端施加均布载荷,一般退火炉中L=18m,假设炉内温度均匀.取退火炉内高温态带钢弹性模量为1.05GPa,泊松比为0.23,屈服应力21.12MPa,材料各向同性.

图4中S50,S100分别表示距O-O切线下方50mm,100mm位置的带钢横截面.从计算结果来看,为防止跑偏采用的锥形导向辊加大了炉内带钢张应力横向分布不均匀性,带钢辊肩位置存在应力集中,此处确实能导致横向瓢曲模态的产生,这也是现场常见的起皱区域.而在辊肩以外对应的带钢区域,张应力急速衰减.本文在下面的解析分析中将采用高斯函数来描述张应力的不均匀分布.

2 横向瓢曲前屈曲弹性应力场

通用的有限元仿真计算工时长,难以形成直观的数学模型来预测炉内瓢曲临界张力,这样既不利于电气控制也难以形成可操作的量化工艺参数.为建立退火炉内不同厚度、宽度、材质带钢的热瓢曲临界张力模型,在此对横向瓢曲行为展开了板内应力场分布和屈曲临界条件的解析研究.根据业内学者的惯用做法,先截取退火炉内无限长带钢在相邻导辊间一段长度作为研究对象,将这一特殊板形缺陷形式的屈曲模态和边界张力条件抽象为数学模型.

2.1 位移函数

基于现场“热瓢曲”浪形模态在宽度方向挠度急剧变化这一基本几何特征,在挠度函数中引入Hermite多项式H(x)作为宽度方向的衰减函数,使其满足条件H(x=0)=1,H(x=±bw)=0,H'(x=±bw)=0,坐标如图5所示.同时,基于浪形屈曲区域具有局部瓢曲的特征,在挠度函数表达式中设定屈曲区域为一待定的局部区域Ω1[-bw,bw;-aw,aw],而非整个方板区域Ω[-b,b;-a,a]在长度方向表现为一个半波,在宽度方向挠度函数需满足如下位移约束条件

在长度方向挠度函数需满足的位移条件为

根据现场瓢曲带钢的实测与局部宽度拉伸实验结果的瓢曲浪形缺陷形式,将其瓢曲浪形的挠度函数描述成下式形式

式中,,为衰减函数;A为挠度幅值,mm.

2.2 边界张力

取边界张应力分布为高斯函数分布形式,既可以描述锥形辊平直段大小,又能反映载荷集中分布的影响,边界载荷分布为

式中,N0为边界不均匀分布载荷尖峰值,N/mm;kN为带钢的整体张力系数;bN为载荷半值宽,即x=bN时,,其值描述了锥形辊凸度宽度,mm.

从式(4)可以看出边界载荷表达式分两部分描述,第1项表示带钢宽度方向不均匀载荷分布,当半值宽bN很小时,第1项带钢边部张力接近0,而中部张力值为N0;第2项表示施加在带钢上的整体平均张力水平.

2.3 前屈曲应力场求解

为了求解其内部应力场分布,可先假设材料为线弹性.其力学基本方程如下:

①平衡方程

②本构方程

③应变协调方程

④边界条件

对于平面矩形薄板由边界载荷诱导产生的应力场求解,可以引入Airy应力函数(x,y),如式(9),使得其满足上述方程,从而构成了弹性应力场求解的数学提法

各应力分量与Airy应力函数(x,y)存在如下关系

(x,y)是重调和函数,由式(6)、式(7)和式(10)可得到协调方程

因此,当前的求解目标是获得满足式(10)、式(11)和边界条件式(8)的应力函数(x,y),根据Timoshenko最小功原理[14],可知单位厚度矩形板的应变能为

采用变分法求上式的极小值,即可得到关于应力函数(x,y)的协调方程式(11).因此,满足边界条件式(8)并使应变能式(12)为极小值时即为所求的应力函数.为此,取级数形式的应力函数为

在此采用伽辽金法求解应力函数,结合格林积分公式可知

式中Ω为矩形积分区域[-b,b;-a,a],mm.

通过求解上式即可获得弹性应力场函数,即前屈曲应力场分布.

2.4 前屈曲应力场分布

当取半板宽b=1个单位长度,半板长a=2单位长度,板厚h=1个单位厚度,张应力半值宽bN=0.3单位长度,整体张力kN=0,不均匀张力幅值N0=1单位线载荷时,带钢矩形板内的应力场分布情况如图6所示.

3 弹性屈曲临界条件求解

3.1 基本方程

对于理想平板带钢承受平面内边界载荷时,当外载荷比临界载荷低时,平板保持原有的平面形态;当外载荷等于临界载荷时,板既可在原先的平面位置保持平衡,也可以在新的微弯状态下达到平衡.后者是板处于临界状态的标志,因此,建立面内载荷作用下板在微弯状态下的平衡微分方程是获得临界载荷的直接有效途径.在此,采用薄板弯曲理论中的Kirchhoff假定.

根据物理方程、几何方程以及弯曲力和弯曲应力的关系建立弹性曲面微分方程如下

式中Nx=σxh,Ny=σyh,Nxy=τxyh,单位为N/mm.

在此取,,,为N0取单位载荷时的弹性应力场.

由于横向屈曲的屈曲区域具有局部瓢曲特征,因此,假设屈曲区域在一个待定的局部区域Ω1[-bw,bw;-aw,aw],见图5所示.

对于屈曲区域Ω1[-bw,bw:-aw,aw],可以认为这个区域四边为简支.

因此,由式(14)和式(3)构成了求解屈曲载荷N0的数学提法.为获得问题的求解,在此采用伽辽金虚位移原理解法,有

此方程可以采用通用的符号运算软件进行求解.

3.2 屈曲区域求解

屈曲区域的求解可以归结为从已知应力场分布的矩形带钢搜索一个浪形屈曲区域Ω1[-bw,bw;-aw,aw],并使得带钢在此区域内屈曲时获得的载荷最小,即求解使得屈曲载荷N0最小的屈曲区域参数aw,bw,从而求得屈曲临界载荷值.据此,假设最小屈曲载荷对应的屈曲区域为:Ω1[-bcr,bcr;-acr,acr].[-bcr,bcr;-acr,acr].

3.3 算例分析

表1比较了3种长宽比(a/b)下,带钢的屈曲载荷可知,a=500mm,b=500mm时取得的无量纲屈曲临界载荷最小;带钢所受整体张力水平对屈曲临界载荷影响非常显著,张力水平越大,临界载荷越大,带材越难发生屈曲;而载荷半值宽对临界载荷的影响还与平均张应力水平有关,因为载荷半值宽增加时,实际也增加了平均张力水平,因此,载荷半值宽的变化涉及到多变量耦合的影响.

4 横向瓢曲的试验分析

为验证带钢宽度方向上不均匀载荷对横向瓢曲的影响,在此借鉴了冲压领域中采用方板对角拉伸实验(Yoshida buckling test,YBT[15])来描述薄板冲压时不均匀拉伸引起的皱曲行为.本文选择矩形方板的局部拉伸实验来近似反映横向瓢曲行为.如图7所示,将载荷施加于钳口夹持端这一局部宽度范围内以实现矩形板的宽向不均匀拉伸,同时,由于受实验设备限制,将所取带钢样本长宽比例缩小,为避免应力尖峰,夹持端与试件连接处用圆弧倒角过渡,试件材料为某厂1420连续退火后的普碳带钢.

表2为常温下板的瓢曲实验结果.从两种长宽比矩形板的局部宽度拉伸实验可以看出,长宽比为1时,矩形带钢的屈曲临界应力比较小,能看到明显的屈曲模态,其板中心面外位移与拉应力的曲线如图8所示.当方板规格为250mm×80mm×0.23mm时,矩形板在所施加载荷达到307 MPa条件下,其获得的面外位移仍不如80 mm×80 mm×0.23 mm方板明显,加载的钳口夹持区域已屈服变形而停止了拉伸,这也说明了其屈曲临界点高于前者.由于试验条件的限制,解析模型的加载条件和设计的钳口加载方式并不能完全吻合,综合各因素,表中解析法和试验结果存在17.3%的相对误差,对工程应用具有很好的指导作用.

5 结论

本文通过有限元仿真分析了退火炉内带钢热瓢曲产生的原因,并对退火炉内带钢热瓢曲行为进行了半解析、半数值研究,获得了其前屈曲应力场、临界屈曲载荷的求解,鉴于篇幅所致,在此未对后屈曲生成路径展开讨论,并运用试验手段验证了解析计算结果的可靠性,为现场抑制和消除热瓢曲缺陷提供了理论参考,文中建立的关于薄带材屈曲行为的求解方法可为其他浪形缺陷形式甚至薄板冲压起皱行为提供一种求解思路.

摘要：薄宽带钢的横向瓢曲是一种产生机理与解决对策都未知的板形缺陷.基于艾利应力函数和Tim★shenko最小功原理建立了横向瓢曲前屈曲变形的力学模型,获得了其前屈曲应力场分布的表达式;运用伽辽金虚位移原理解法获得了其临界屈曲载荷.研究结果可为连退线上板带材生产过程中横向瓢曲的预测和控制提供参考依据.

关键词：热瓢曲,临界应力,屈曲,退火炉,带钢

横向变形 篇2

作为铁路轨道结构中最重要的组成部件，钢轨的受力情况十分复杂，在车轮荷载的作用下，钢轨会产生竖向弯曲等变形。

在以往的研究当中，一般是基于钢轨经典力学分析的方法，将钢轨假定为无限长梁，点支承或连续支承于下部轨枕或地基上，如图1所示。由于为单一长梁，传统的方法仅能够求解钢轨竖向、横向等整体截面的位移及应力，而不能反映由于偏心荷载等作用引起的扭转应力以及钢轨内部局部的应力变化。

对于钢轨截面分析的一般理论求解，正常也仅考虑竖向荷载作用于钢轨中心对称位置的轨头面上，针对钢轨的扭转变形，则主要是考虑由横向作用力引起。张永兴等为了细化钢轨的水平位移，在研究钢轨扭转变形过程中，考虑了竖向荷载偏心的影响，但为减少计算量，钢轨截面简化成方方正正的工字形状，即头部、腰部及底部视为三个矩形，分别计算三部分的扭转及惯性矩，通过计算由扭转引起的横向位移与横向力引起的弯曲位移叠加求得整体的水平位移。事实上，即使在直线区段，车轮对钢轨的竖向载荷也并非作用于钢轨头部中心位置。图2为现场调查的某轨道结构直线段钢轨上的光带示意图，可见轮轨接触位置偏向钢轨内侧，存在相应的偏心，同时，作用力也并非单点作用，而是形成具有一定宽度的轮轨接触面。

传统的简化方式已经越来越不适应高速铁路高精度的要求，且随着计算机技术的不断发展，对钢轨内部应力变化以及轨头局部位置位移的求解成为可能，本文针对直线上钢轨的竖向偏心荷载进行力学分析，采用有限元计算软件将钢轨考虑为三维实体，竖向偏心荷载按单一作用力施加于轨头内侧一定位置，以模拟和实际轨道结构更为近似的受力模型。

2 力学模型

对于实际铁路轨道，钢轨由扣件弹条扣压，按轨枕间距铺设在轨枕之上，高速铁路无砟轨道结构的钢轨则是按照承轨槽间距铺设于轨道板上，钢轨底部由扣件弹性垫板支撑。因此，钢轨与枕下结构的接触实为面接触状态。考虑上述因素，本文将扣件系统对钢轨的作用视为一定宽度的均布弹簧支撑，扣件以下部分视为全约束，建立钢轨的三维实体力学模型，如图3所示。

钢轨截面受力示意图如图4所示，受竖向偏心荷载F作用，e表示偏心值，则钢轨受到扭矩大小为Mt=F×e。为了分析钢轨受竖向偏心荷载的影响，本文仅针对直线上钢轨竖向力的作用进行研究，而不考虑横向力以及轮轨产生的蠕滑等切向力作用。

为保证合理的计算精度，且由于仅分析荷载作用位置附近钢轨的受力和变形特点，结合经验及模型的求解结果，取13跨钢轨长为计算长度，随着钢轨跨数的再增加，两端的计算结果影响已经非常小，可忽略不计。事实上，对于一般的静力分析，钢轨跨数超过10跨时就可以满足求解要求，同时，为了消除边界效应等对结果的影响，此处取13跨则为一种相对保守的取法。对于弹性垫板刚度，线性均布弹簧按其单支弹簧刚度均匀分配。

计算参数如下:1)钢轨:采用60 kg/m轨，其弹性模量为2.1×105MPa，泊松比为0.3，线膨胀系数11.8×10-6/℃，密度7 830 kg/m3。2)扣件系统:扣件弹性垫板的静刚度为22.5 k N/mm，垫板宽度为150 mm，轨枕间距0.65 m。3)荷载:荷载取为22.5 k N。4)偏心值大小:按无偏心、偏心5 mm、偏心10 mm、偏心13.86 mm，此处最大偏心值结合建模时节点位置选取，并无实际意义。直线上钢轨过大的偏心值一般不会存在，因此仅取到不足15 mm，曲线上钢轨则会存在较大的偏心，但同时由于超高的设置，也得需要考虑横向力作用。

3 计算结果及分析

钢轨采用三维实体建模的方式求解可以直接反映钢轨表面及内部的应力值及变形的大小，本文对作用力范围内的一跨钢轨展开研究，选取几种典型的计算结果进行分析，主要包括偏心荷载作用下钢轨轨距点之间的横向位移变化以及钢轨由于扭转产生的附加正应力值的大小两个主要方面。

图5，图6为选取的钢轨的分析截面及点位的示意图。图6中，b为轨距点，c为下颚点，e为轨腰中心位置，h为轨底中心，其他点则为选取的一些最有可能产生钢轨病害的过渡位置。

3.1 荷载作用截面轨头部位的横向位移

由于钢轨横向位移的变化主要体现为轨距的扩大或变小，因此，只对轨头部分a,b,c三点的横向位移进行分析，在无偏心和不同偏心荷载下的位移结果如图7所示。其中，横轴表示偏心值大小，纵轴表示钢轨横向位移值。

从图7中可以看出，随着偏心值的不断增大，钢轨的横向位移也随之增大，最大值已经达到了2 mm左右，同时，轨头下颚点(c点)也达到了1.5 mm左右，可见竖向偏心荷载对钢轨的横向位移有着不可忽视的影响，会导致钢轨轨头向内侧扭转，从而减小了轨距，增大了轮缘撞击轨角及轨头侧面的风险，不可避免地会引起侧磨、轨角及下颚裂纹等伤损。

3.2 偏心荷载作用跨轨距点(b点)横向位移比较

为了反映出偏心荷载对其作用跨钢轨整体的横向位移影响，对轨距点(b点)在一跨内的横向位移进行了分析，如图8所示。定义作用点坐标为0，则中间跨钢轨位置范围为-325 mm～325 mm，横轴表示其位置坐标。

从图8中可以看出，在偏心荷载的影响下，轨头处的横移量受其影响较大，且一跨范围内的整体横移值变化较小，说明其影响值已经超过了一跨钢轨长度，钢轨轨头整体侧移，产生了一定程度的扭转。

3.3 偏心荷载作用对钢轨扭转的影响分析

偏心荷载产生扭矩，因此对钢轨影响的一个主要体现就是产生了扭转作用，钢轨的扭转会导致截面内部正应力值发生相应的变化，因此，偏心荷载对钢轨扭转的影响可以从钢轨截面正应力角度进行分析。为便于分析，将轨头部分和其他各分析点的计算结果分别处理，如图9，图10所示。

从图9中可以看出，在偏心荷载的作用下，轨头区钢轨正应力呈逐渐增加趋势，相对于无偏心时的正应力值，应力增大值即由于扭转产生，为扭转正应力。由于轨头部分受压，初始正应力值为负，扭转作用增加了正应力值，导致压应力值逐渐变小，轨头下颚c点增大为正值，由受压变成受拉状态。结合钢轨产生向内侧的横向位移，轨距变小增加了车轮对钢轨的侧面撞击，在两种因素的影响下，使得钢轨侧磨更为严重，同时由于出现拉应力而增加了下颚处产生裂纹的可能。

图10反映出钢轨扭转会导致轨腰中部以下钢轨截面的正应力值减小，即扭转正应力为负值，但钢轨偏心产生扭曲影响的影响值较小，轨底中心处的应力值大小基本不变，仅在轨腰中心附近出现一定的拉压应力波动。

MPa

偏心荷载引起的钢轨截面的正应力增加值即是由于钢轨的扭转产生，扭转正应力如表1所示。

从表1中也可以看出，偏心荷载产生了钢轨的扭转，且轨头部分的扭转正应力值较大，其中下颚点处扭转正应力达到最大值，且为拉应力。

4 结语

通过对直线上钢轨偏心荷载作用下轨头位置横向位移及截面应力值的求解，可以得出以下结论:

1)车轮竖向偏心荷载会导致钢轨产生横向位移及扭转变形，对于越来越高的精度要求，即使直线上的钢轨，也需要考虑竖向荷载偏心大小的影响。

2)偏心荷载导致轨距减小，增加了车轮对钢轨内侧撞击的概率，同时，钢轨由于扭转产生了正的扭转应力，导致轨头部分的钢轨截面正应力增加，在下颚处甚至出现了拉应力。两者的叠加会导致轨角内侧钢轨磨耗更加剧烈，甚至出现裂纹。

3)钢轨扭转会产生一系列的问题，若考虑横向力等作用，扭转应力又会重新分布，对钢轨扭转问题还值得进一步的研究。

摘要：采用有限元计算软件,建立了钢轨的三维实体力学模型,分析了竖向荷载在不同偏心值下钢轨的横向位移及扭转应力,结果表明,钢轨的扭转应力会对钢轨产生不利的影响,精度要求较高时应考虑钢轨的扭转变形和应力。

关键词：钢轨,三维实体,偏心荷载,扭转应力

参考文献

[1]李成辉.轨道[M].成都:西南交通大学出版社,2005.

横向变形 篇3

一般曲线坐标系下的二维水沙模型,当涉及河床变形及岸坡崩退时,崩岸的宽度与计算网格的尺寸存在一定的差距。如要精确考虑崩岸尺寸,崩岸处网格必须尽量小,但对于整个河道而言,会导致计算速度过慢。为解决实际存在的问题,国内外一些学者进行了移动坐标下可变网格的研究,并取得了较为理想的成果。日本Nagata[1]等学者首先提出了移动坐标下模拟河岸侵蚀过程的数值方法。Jang和Shimizu[2]采用移动坐标下的平面二维水沙数学模型对宽浅型的水流形态,河床变形及河岸的侧向冲蚀变形进行了数值模拟研究。杨建明[3]采用可变网格技术模拟挑流冲刷时因河床变形引起的流场变化。

河流弯道水流呈现明显的三维特性,平面二维水沙模型无法很好地模拟由于弯道横向输沙引起的凹岸冲刷、凸岸淤积的现象。三维数学模型的研究虽然取得了一些成果,但许多技术性问题尚待解决[4];另外一种途径就是对现有的平面二维数学模型进行修正[5],现有的研究成果在文献[5]中作了比较详细的阐述。本文在文献[7]建立的曲线平面二维数学模型基础之上,引入文献[8]中移动坐标,并采取文献[5]的研究方法,考虑悬移质的横向输沙,对河岸岸坡崩塌的模拟,引入夏军强[6]改进后的黏性土河岸变形计算方法。通过上述这些方法对已有平面二维数学模型进行修正。

1 数学模型

1.1 基本方程

移动坐标下,方程中的沿时间偏导项为:

$\frac{\partial }{\partial t}=\frac{\partial }{\partial \tau }+\frac{\partial }{\partial \xi }\frac{\partial \xi }{\partial \tau }+\frac{\partial }{\partial \eta }\frac{\partial \eta }{\partial \tau }(1)$

因此在移动一般曲线坐标下,平面二维水沙数学模型的基本方程如下[8]:

连续方程:

$J\frac{\partial h}{\partial \tau }+\frac{\partial (U+\xi {}_t)h}{\partial \xi }+\frac{\partial (V+\eta {}_t)h}{\partial \eta }=Jq(2)$

ξ方向运动方程:

$\begin{array}{l}J\frac{\partial hU}{\partial \tau }+\frac{\partial (U+\xi {}_t){\rm M}}{\partial \xi }+\frac{\partial (V+\eta {}_1){\rm M}}{\partial \eta }=\\ -gh(\frac{\partial {\rm Z}}{\partial \xi }\xi {}_x+\frac{\partial {\rm Z}}{\partial \eta }\eta {}_x)J+\frac{\partial }{\partial \xi }[DJ(q{}_{11}{\rm M}{}_{\xi }+q{}_{12}{\rm M}{}_{\eta })]+\\ \frac{\partial }{\partial \eta }[DJ(q{}_{12}{\rm M}{}_{\xi }+q{}_{22}{\rm M}{}_{\eta })]-\frac{gn{}^2{\rm M}J\sqrt{u{}^2+v{}^2}}{h{}^{4/3}}+Jqu{}_0(3)\end{array}$

η方向运动方程:

$\begin{array}{l}J\frac{\partial hV}{\partial \tau }+\frac{\partial (U+\xi {}_t){\rm N}}{\partial \xi }+\frac{\partial (V+\eta {}_1){\rm M}}{\partial \eta }=\\ -gh(\frac{\partial {\rm Z}}{\partial \xi }\xi {}_x+\frac{\partial {\rm Z}}{\partial \eta }\eta {}_y)J+\frac{\partial }{\partial \xi }[DJ(q{}_{11}{\rm N}{}_{\xi }+q{}_{12}{\rm N}{}_{\eta })]+\\ \frac{\partial }{\partial \eta }[DJ(q{}_{12}{\rm N}{}_{\xi }+q{}_{22}{\rm N}{}_{\eta })]-\frac{gn{}^2{\rm N}J\sqrt{u{}^2+v{}^2}}{h{}^{4/3}}+Jqv{}_0(4)\end{array}$

泥沙连续性方程:

$\begin{array}{l}J\frac{\partial hS{}_k}{\partial \tau }+\frac{\partial (U+\xi {}_t)hS{}_{{\rm K}}}{\partial \xi }+\frac{\partial (V+\eta {}_t)hS{}_{{\rm K}}}{\partial \eta }=\\ \frac{\partial }{\partial \xi }[DJ(q{}_{11}hS{}_{{\rm K}\xi }+q{}_{12}hS{}_{{\rm K}\eta })]+\\ \frac{\partial }{\partial \eta }[DJ(q{}_{12}hS{}_{{\rm K}\xi }+q{}_{22}hS{}_{{\rm K}\eta })]-J\alpha W{}_{{\rm K}}(S{}_{{\rm K}}-S{}_{*{\rm K}})+Jqs{}_0(5)\end{array}$

河床变形方程:

${\gamma }^{\prime }\frac{\partial {\rm Z}{}_b}{\partial \tau }={\displaystyle \sum \alpha }\omega {}_{{\rm K}}(S{}_{{\rm K}}-S{}_{*{\rm K}})(6)$

式中:x、y为物理平面坐标;ξ、η为贴体坐标下的坐标;D为紊动黏性系数;n为Manning糙率系数;u、v分别为x 和y方向的流速,m/s;J为雅可比变换系数,即J=xξyη-xηyξ;U和V分别为曲线坐标中流速在ξ、η方向的分量,即U=ξxu+ξyv, V=ηxu+ηyv;ξt=-xtξx-ytξy,ηt=-xtηx-ytηy;h为水深,m;z为水位,m;q为单位面积上水流的源汇强度,m/s;M=uh, N=vh;q11=ξ${}_x^2$+ξ${}_y^2$,q12=ξxηx+ξyηy;q22=η${}_x^2$+η${}_y^2$;Mξ、Mη、Nξ、Nη、Skξ、Skη表示偏导数;ρ为水体密度,g/cm3;γ′为泥沙干密度,kg/m3;α为恢复饱和系数;ε为泥沙扩散系数;ω为泥沙颗粒沉速,m/s;ωK、SK、S*K分别对应分组泥沙颗粒沉速、含沙量和挟沙力。

本文引用文献[5]的研究成果,即在泥沙连续方程(5)中加入横向输沙项,根据图1得到修正后的泥沙连续性方程如下:

$\begin{array}{l}J\frac{\partial hS{}_{{\rm K}}}{\partial \tau }+\frac{\partial (U+\xi {}_t)hS{}_{{\rm K}}+Q{}_{n{\rm K}\xi }}{\partial \xi }+\frac{\partial (V+\eta {}_t)hS{}_{{\rm K}}-Q{}_{n{\rm K}\eta }}{\partial \eta }=\\ \frac{\partial }{\partial \xi }[DJ(q{}_{11}hS{}_{{\rm K}\xi }+q{}_{12}hS{}_{{\rm K}\eta })]+\\ \frac{\partial }{\partial \eta }[DJ(q{}_{12}hS{}_{{\rm K}\xi }+q{}_{22}hS{}_{{\rm K}\eta })]-J\alpha W{}_{{\rm K}}(S{}_{{\rm K}}-S{}_{*{\rm K}})+Jqs{}_0(7)\end{array}$

式中:QnKξ、QnKη曲线坐标下沿ξ、η方向的横向输沙率。

1.2 模型辅助方程及相关问题的处理

(1)水流挟沙力公式。

采用与实际吻合较好的张瑞瑾公式:

$S{}_{*}=k\left(\frac{\overline{u}{}^3}{gh\overline{\omega }}\right){}^m(8)$

式中:$\overline{u}$为垂线平均流速,m/s;g为重力加速度,m/s2;h为水深,m;$\overline{w}$为泥沙平均沉速,m/s;k、m为系数。

其中,$S={\displaystyle \sum _{{\rm K}=1}^{n}S}{}_{{\rm K}},S{}_{*}={\displaystyle \sum _{{\rm K}=1}^{n}S}{}_{*{\rm K}}$

(2)横向输沙率公式。

采用张红武公式[9],具体形式如下:

$Q{}_{n{\rm K}}=\frac{75.86\overline{u}hS{}_{{\rm K}}}{{\rm N}{}_{0}r}\exp \left(0.0931\frac{\omega }{\kappa u{}_{*}}\arctan \sqrt{\frac{1}{\eta }-1}\right)f(\eta )(9)$

其中:

$\begin{array}{l}f(\eta )=[(1+5.75\frac{g}{C{}^2})\eta {}^{1.875}-0.88\eta {}^{2.14}+(0.034-\\ 12.5\frac{g}{C{}^2})\eta {}^{0.857}+4.72\frac{g}{C{}^2}-0.088]\end{array}$

${\rm N}{}_0={\displaystyle \underset{0}{\overset{1}{\int }}\eta }{}^{\frac{1}{7}}\exp \left(0.0931\frac{\omega }{\kappa u{}_{*}}\arctan \sqrt{\frac{1}{\eta }-1}\right)\text{d}\eta$

式中:r为流线曲率半径;C为Chezy系数;u*为摩阻流速。

(3)分组挟沙力及挟沙力级配。

先用挟沙力公式求出非均匀悬移质泥沙的总挟沙力,然后用挟沙力级配乘以总挟沙力,得到分组挟沙力公式:

S*K=PK*S*

挟沙力级配与泥沙来源的级配有关。水中泥沙有两个来源,一是上游来沙,而是床面补给,挟沙力级配采用如下公式:

${\rm P}{}_{{\rm K}}=\frac{S{}_{*{\rm K}}^{´}+S{}_{{\rm K}}}{{\displaystyle \sum _{{\rm K}\ge {\rm K}{}_d}(}S{}_{*{\rm K}}^{´}+S{}_{{\rm K}})}$

(4)床沙级配采用韦直林模式。

将河床组成概化为表层的交换层、中层的过渡层及底层的泥沙冲刷极限层。规定在每一个计算时段内,各层间的界面都固定不变,泥沙交换限制在表层进行,中层和底层暂时不受影响。在时段末,根据床面的冲刷或淤积往下或往上移动表层和中层,保持这两层的厚度不变,而令底层厚度随冲淤厚度的大小而变化。

1.3 方程求解

模型采用水流泥沙非耦合方法求解,计算物理量布置于同位(非交错)网格,微分方程的数值离散采用有限体积法(控制容积法),离散方程的求解基于Pantankar压力校正法原理[10],具体计算步骤如下。①二维水流计算。②根据水深判断各断面是否属于滩岸,并对岸坡进行崩塌稳定性分析。计算土体崩塌结果,包括单位堤防崩塌土体的宽度及崩塌体积。横向移动岸坡崩塌处的计算网格,计算下一步时泥沙侧向输入项及河床地形高程。③悬移质泥沙扩散方程的求解。④根据河床变形方程计算河床变形,并修正河床高程。⑤进行下一时段的计算。

2 算例及结果分析

由于目前缺少河岸冲刷过程以及弯道内悬移质泥沙不平衡输移的实测资料,故以一概化的90°弯道的河床变形作为模拟对象,模拟凹岸的冲刷后退,对结果进行定性分析。

2.1 模型概化

假定有一90°概化弯道,顺直进口段350 m,90°弯道550 m,顺直出口段800 m,河宽100 m。初始河岸高程为8.0 m,出口河底高程为0.0 m,河床初始比降0.000 2,初始糙率0.020。模拟河段共划分网格数为58×21,其中ξ方向主槽平均网格大小约为30 m,η方向主槽平均网格大小约为5 m。如图2、3所示。边界条件:进口断面给定流量600 m3/s,出口水位7.0 m。

计算流场如图2。在弯道顺直进口段,纵向平均流速沿河宽基本呈对称分布。进入弯道,凸岸流速稍增。至弯顶以下,出现相反的调整,流速渐渐移向凹岸。出弯水流,在相当长的距离内,最大流速仍靠近外侧河岸。计算结果与张红武[11]弯道模型流速平面分布结论一致。说明本文中所采用的平面二维水流数学模型计算方法及相关处理是可行的。

2.2 二维泥沙计算

床沙平均粒径为0.087 mm,河岸土体的特性为:容重γs为18 kN/m3、内摩擦角Φ为14°、凝聚力C为20 kN/m2,堤岸土体的起动切应力τc为1.2 N/m2。边界条件用恒定流计算。进口断面给定流量Q=600 m3/s ,出口水位Z为7.0 m。进口来流中悬移质含量为0 kg/m3;挟沙力系数k=0.15,m=0.92;恢复饱和系数冲刷时取α=1.0,淤积时取α=0.25。本模型模拟了凸岸边滩形成后,弯道水流和泥沙特性。模型概化计算地形如图4、5所示。

(1)冲淤变化。图6给出典型断面CS1不同时刻形态图。从图中可以看出,由于横向输沙的作用,导致凹岸冲刷后退,凸岸淤积。随着时间的推移,凹岸不断冲刷,导致河床岸坡失稳坍塌,河道展宽。图7(a)～(d)为T=0 d,10 d,20 d,30 d时刻地形云图。

(2)岸坡崩塌。图8(a)～(b)给出始末时刻崩塌变化下的网格移动的相对位置。从图中看出,凹岸在横向环流的冲刷下逐渐失稳坍塌,最大崩塌后退距离发生在顶冲点附近位置。从图9各时刻崩塌后退位置示意图可看出,随着主槽冲深下切及凹岸的冲刷后退,崩塌位置也逐渐移向下游。

上述计算结果与现有的研究结论是吻合的,说明扩展后的模型具有一定的可行性,能够反映弯道环流横向输沙以及由此产生的凹岸的崩塌后退。移动网格技术的应用,解决了崩岸宽度与计算网格的尺寸存在的矛盾,克服了采用固定网格进行计算时难以精确调整模型计算边界的困难。

下一步将在本文的基础上,结合具体河段的实测资料,将修正后的二维水沙模型应用于实际河道,使其能够较好地模拟弯道水流的水沙运动规律及河岸的崩塌现象。

3 结 语

(1)采用移动网格技术,能够较好的解决崩岸宽度与计算网格的尺寸存在的矛盾,并能精确的调整模型计算边界。在保证计算速度的前提下,提高计算精度。

(2)扩展后的平面二维水沙模型引入横向输沙项,近似考虑了弯道环流所引起的横向输沙。使之更好的反映实际河道的冲淤变化。

(3)采用修正后的模型,对一概化弯道进行计算,计算结果与现有的认识一致,证明该模型具有一定的可靠性和实用性。

参考文献

[1]N Nagata,T Hosoda,Y Muramoto.Numerical Analysis of River图8网格位移图图9各时刻崩塌后退位置示意图Channel Processes with Bank Erosion[J].Journal of HydraulicEngineering,2000,126(4):243-252.

[2]Chang-Lae Jang,Yasuyuki Shimizu.Numerical Simulation of Rel-atively Wide,Shallow Channels with Erodible Banks[J].Journalof Hydraulic Engineering,2005,131(7):565-575.

[3]杨建明,吴建华.动网格技术数值模拟挑流冲刷过程[J].水动力学研究与进展A辑,2001,16(2):6.

[4]刘玉玲,刘哲.弯道水流数值模拟[J].应用力学学报,2007,24(2).

[5]钟德钰,张红武.考虑环流横向输沙及河岸边形的平面二维扩展数学模型[J].水利学报,2004,(7).

[6]夏军强.河岸冲刷机理研究及数值模拟[D].北京:清华大学,2002.

[7]谢作涛,张小峰.一般曲线坐标系平面二维水沙数学模型研究与应用[J].泥沙研究,2005,(6):12.

[8]袁晶.移动坐标下平面二维水沙数学模型研究及应用[D].武汉:武汉大学,2006.

[9]张红武,吕昕.弯道水力学[M].北京:水利电力出版社,1993.

[10]帕坦卡S V.传热与流体流动的数值计算[M].张政译.北京:科学出版社,1984.