非稳态特性范文

非稳态特性范文（精选8篇）

非稳态特性 第1篇

2009年1月21日收到 轮胎实际总是处于非稳态工况,因此为了提高车辆动力学的仿真精度必须建立精确的轮胎非稳态模型。迄今,研究者建立了许多对各种非稳态特性的描述方法[1]。但这些方法在考虑模型对任意工况的适用性、如何调整松弛长度及高速低速特性无缝结合这3个问题时,都或多或少存在一些问题。本文基于当前工作点轮胎物理模型描述,提出了一种全新的轮胎非稳态特性建模方法。

1 当前工作点轮胎非稳态特性物理描述

轮胎稳态模型Fs～S反映的是胎体变形稳定的情况下胎面与路面间相互作用的特性,故也称之为胎面模型,表达了轮胎力与滑移率间的非线性关系,如图1所示。而滑移率又可表示为滑移速度的函数,因此轮胎力与滑移速度间的关系可用阻尼器来等价,该线性阻尼器称为胎面阻尼器Fd～Vs,如图2所示,其阻尼系数C0由胎面模型在当前工作点的特性确定。

当前工作点是指将胎面模型等效为胎面阻尼器时,胎面模型的力与滑移率完全等价于胎面阻尼器的力与滑移速度的点,在胎面模型中表示为点(S0,Fs0),胎面阻尼器中表示为点(Vs 0 ,Fd0 ),且Fs0=Fd0。

胎体的松弛效应是影响轮胎非稳态特性的最重要因素,在本文中,对胎体的建模Fc～u是将胎体考虑为两向拉压弹簧,表示胎体的侧向及纵向平移变形特性。

在将胎面模型等效为胎面阻尼器、胎体特性考虑为两向拉压弹簧后,轮胎非稳态特性的物理描述如图3所示:轮心处的名义侧向或纵向输入Vsn经胎体松弛后产生胎面有效输入Vs,进而产生轮胎侧(纵)向力及回正力矩。

2 非稳态特性建模理论与求解流程

基于轮胎非稳态特性物理描述,当前工作点轮胎非稳态特性用图4所示的物理模型来等价表达,其系统方程为:

(1)式中,u为胎体平移变形,Kc为胎体平移刚度,Vs为胎面阻尼器相对滑移速度。该方程可用于精确描述离开当前工作点后轮胎物理模型的特性,并以此来近似描述轮胎真实特性。

如图1～图2所示,胎面阻尼器当前工作点时的阻尼系数C0可从下面推导中得到:

K0为胎面模型中,连接原点与当前工作点直线的斜率,当轮胎胎面模型位于线性工作区时,K0为轮胎的侧偏或纵滑刚度。顾及滑移率的定义有S0=Vs0/Vr,Vr为印迹更新速度,则从式(2)可知:

定义l0=K0/Kc为当前工作点松弛长度,容易看出当轮胎工作于线性区时,l0就是轮胎的松弛长度。引入l0则C0表示为:

顾及Vs与轮心处的名义输入Vsn及胎体变形速度$\dot{u}$的关系及式(4),从式(1)可推得:

利用式(5)就可以求得离开当前工作点一个时间步长后胎体变形u,进而求得轮胎非稳态力。但正如前文所述,此时求得的非稳态输出都是基于当前工作点衍生的近似系统对真实轮胎特性的预测,代表真实轮胎的胎面模型与线性阻尼器是有明显区别的,如图1～图2所示,因此为了得到轮胎的特性,需要用胎面模型进行修正。图5所示为此非稳态模型的求解流程。

3 TMPT部分非稳态特性仿真结果

TMPT(Tire Model Performance Test)是2002～2007年由国际车辆动力学协会(IAVSD)组织的国际性“轮胎模型性能测试”活动[2]。郭孔辉院士及其课题组开发的UniTire模型参加了此次活动[3]。本文的轮胎非稳态建模方法已应用于UniTire模型,图6～图7给出了TMPT部分非稳态特性仿真结果。

4 结论

(1)所建模型适用任意轮胎非稳态操稳工况。

(2)模型具有非定长松弛长度特性,可描述精细非稳态特性。

(3)所述的建模方法可用一套方程来统一表达低速、高速非稳态问题,无高速低速间不同模型切换问题。

摘要：为合理表达非稳态工况下轮胎力学特性,基于当前工作点轮胎物理模型描述,提出了一种全新的轮胎非稳态特性建模方法,给出了模型求解流程,并展示了TMPT部分非稳态工况仿真结果。所提出的模型具有适用于任意非稳态工况、非恒定松弛长度、高速与低速特性统一表达的特点。

关键词：轮胎模型,非稳态特性非恒定松弛长度高速特性低速特性

参考文献

[1]Pacejka HB.Tyre and vehicle dynamics.Butterworth-Heinemann,2002

[2]Pl chl L P.Tire model performance test(TMPT).Taylor&Francis Group,2007:225

非稳态特性 第2篇

顺序函数法求解二维非稳态热传导逆问题

基于有限控制体积法对二维非稳态热传导问题的.数值模拟,导出了处理二维非稳态热传导逆问题的顺序函数法,并用该方法来对一典型的圆环域边界条件反演问题进行了反演计算.结果表明,顺序函数法是求解二维非稳态热传导逆问题的有效方法,当测量噪声比较小时,顺序函数法能得出较高精度的反演结果;当测量噪声增大时,反演结果仍能较好地再现出精确解在空间和时间方向上的变化趋势,具有一定精度.

作 者：钱炜祺 蔡金狮 作者单位：中国空气动力研究与发展中心总体技术部,四川绵阳,621000刊 名：空气动力学学报 ISTIC EI PKU英文刊名：ACTA AERODYNAMICA SINICA年，卷(期)：200220(3)分类号：V211.3 O357.5+3关键词：二维非稳态热传导逆问题 顺序函数法 热流密度

液晶光阀的双稳态特性研究 第3篇

液晶作为一种应用广泛的特殊功能材料,既具有强烈的各向异性物理特征,同时又像普通流体那样具有流动性,从宏观物理性质来看, 它既具有液体的流动性、粘滞性, 同时又具有晶体的各向异性, 能像晶体一样发生双折射、布拉格反射、衍射及旋光效应等,也能在外场作用下产生电光或磁光效应。同时还具有热效应、压电效应、光化学效应、光生伏特效应等性质。

液晶光阀(LCLV) 是20世纪70年代发展起来的、利用液晶的光学特性制作的一种空间光调制器,它是信息与激光技术领域的关键电子器件之一, 是近年来应用最多的空间光调制器,能广泛应用于光学信息处理、光信息贮存、大屏幕显示、投影电视、全息照相、图像识别、光学计算机等领域[1,2,3]。

本文测试使用的TB3639型液晶光阀使用的是正性向列相液晶, 液晶分子在电场的作用下, 通过液晶分子重新取向后的扭曲和光学双折射现象而形成光的调制[4,5,6],本文主要测量TB3639型液晶光阀在电压驱动下的光学双稳态特性,并讨论其产生光学双稳态特性的非线性效应物理机制。

1 基本原理

当组成物质的原子或分子的空间尺寸远大于它们随时间变化的尺寸,同时其变化过程相对缓慢时,则可以把组成物质的每个原子、分子的微观状态看做连续体,用它来描述宏观的物理特性,因为液晶宏观物理问题同固体和流体中的问题相似,液晶分子的尺寸远小于它们随时间变化的空间尺寸,其变化过程相对缓慢,有关液晶的许多重要物理现象都可以在忽略了单个组成分子行为的情况下,忽略液晶分子的统计起伏,则可以把液晶当作连续的弹性介质来处理。

液晶是液态的物质,不可能产生像固体那样的形变,但是,描述液晶分子排列取向的指向矢,在外场的作用下可以改变方向,去除外场后,由于液晶分子间的相互作用和表面锚定作用,又有恢复到原来取向的趋势,液晶的这种形变表现为外力作用下液晶指向矢的分布发生变化,这种指向矢取向的变化和固体的弹性形变相类似,同样具有展曲、扭曲、弯曲等形变形式,可以参照固体弹性理论,并考虑外界条件的影响,引用相应的弹性常数来解决液晶的物理问题,这就是液晶连续弹性体理论的基础[7]。

TB3639型液晶光阀是扭曲角为90°的TN型液晶光阀,所用液晶是添加了少量手性材料的向列相液晶,根据液晶的连续弹性体理论,中心分子的倾角会随着外加电场或磁场的变化发生较大的改变。图1是扭曲角达到360°时的一种强锚定的液晶盒外加电场与中心分子倾角的关系曲线,图中曲线的变化关系是采用差分迭代方法[8,9]计算出在具有大预倾角、扭曲角为360°的液晶盒中,液晶层中心部分的液晶分子的倾斜角是怎样随电压变化的,如图中实线所示。

从图1中可以看出,电压在约1.8~2.3 V的变化范围内时,中心分子倾斜角会发生较大变化,电压直到大于2.3 V时中心液晶分子的指向矢才陡然偏向90°,降低电压时,1.8 V前均维持90°,直到1.8 V以下才陡然降到0°左右。由此电压在1.8 V至2.3 V之间为双稳态区间。在此区间,液晶的指向矢可以为0°的状态,也可以是90°的状态。这种由向列相液晶Friderikz转变引起的双稳态需要有持续的电压,在零场下是不能实现的。

2 测试方法

本文测试采用的TB3639型液晶光阀样品由汕头超声仪器有限公司提供,其技术参数:液晶层厚度7.425 μm,驱动条件为5 Vop,1/1Duty,1/1Bias,测试其电光特性曲线采用的设备是UV-VIS8500型紫外/可见分光光度计,波长范围190~1 100 nm,测量其双稳态特性的主要 设备为WGD-8A型组合式多功能光栅光谱仪,波长范围200~660 nm,所用光源为氦氖激光器(波长632.8 nm)。测试室温为25°C。 驱动为频率100 Hz的交流电场,电压0~5 V[4,5]。

测试液晶光阀的双稳态特性的装置如图2所示。

3 结果及分析

3.1 电光特性曲线T-V

由UV-VIS8500型紫外/可见分光光度计测试得出TB3639型液晶光阀在波长632.8 nm处的透过率-电压电光特性曲线(T-V)。图3为由低电压到高电压变化时透射光强与电压的变化关系,图4为由高电压到低电压变化时透射光强与电压的变化关系。

通过比较图3和图4可以得出电压在1.8~2.6 V范围内液晶光阀具有双稳态特性,如图5所示,a为电压由低到高的变化曲线,b为电压由高到低的变化曲线。

3.2 双稳态特性测试

我们选取2.2 V(中间)处作为其双稳态特性的测试电压点,分别使电压由0~2.2 V和5~2.2 V变化,分别测试其透射光谱特性,通过图2设计实验测试光路,用He-Ne激光器作为测试光源,通过调节偏光片P1与P2之间的角度来调节光源的光强大小,在光线通过由电脑控制的WGD-8A型组合式多功能光栅光谱仪测试后可以得出下面图6中的液晶光阀双稳态特性曲线。其中a为电压由低到高的变化曲线,b为电压由高到低的变化曲线。

4 结 论

液晶光阀在某一电压范围内存在双稳态特性,同时扭曲角对液晶光阀的双稳态特性影响较大,由理论计算出的360°扭曲液晶盒的双稳态特性比较明显,而我们使用的TN型液晶光阀的双稳态特性不是很明显。主要原因是由于其扭曲角只有90°导致其产生双稳态特性的某一电压范围内的液晶分子指向矢差别不是很大,因此其双稳态特性不是很明显。通过测试液晶光阀的双稳态特性,为我们设计零场双稳态液晶显示器件提供了重要的参考依据。

参考文献

[1]Hebbron M C,Makh S S.Development of gallium arsenide-based spatial light modulators[C].Proc SPIE,1987,825:19-23.

[2]Wu Shintson,Uzi Efron,Hsu Tsung-yuan.Near-infrared-to-visible image conversion using a Si liquid-crystal light valve[J].Optics Letters,1988,13(1):13.

[3]Yu F Y S,Tai A.Spatial filtered pseudocolor holographicimaging[J].J Opt,1978,9:269-273.

[4]黄,郑泽浩,欧阳艳东,等.液晶光栅三基色的电光特性研究[J].光电子技术,2009,29(3):161-163.

[5]黄,欧阳艳东,吴永俊,等.液晶显示器防紫外光谱特性分析[J].光谱学与光谱分析,2004,24(5):637-639.

[6]黄锡珉.液晶显示技术发展轨迹[J].液晶与显示,2003,18(1):1-6.

[7]黄子强.液晶显示原理[M].北京:国防工业出版社,2008:26-37.

[8]Wang Qian,He Sailing.Simulation and comparison study ofliquid crystal directordistributions[J].Acta Physica Sinica,2001,50(5):926-932.

圆柱状多孔材料非稳态组织收缩 第4篇

上述学者在研究物料干燥体积收缩时,多以物料体积收缩与平均干基含水率的变化为基点,探索干燥过程中物料整体体积的变化规律,而很少对物料内部水分扩散引起的组织收缩进行研究。张京平等[4,5]通过实验测得小块果蔬体积收缩率与平均干基含水率变化之间的关系,然后将小块果蔬体积收缩率与平均干基含水率的变化规律应用于大块果蔬各个质点,并推导建立了球形果蔬边界的收缩方程。Białobrzewski I.等[6]采用Possion公式描述的收缩模型对苹果切片干燥过程中局部密度变化进行了研究。此外,Białobrzewski I.等[7] 还利用Possion公式模型对胡萝卜立方块的组织收缩进行了研究。然而,Possion公式模型不能描述颗粒中心组织滞后于外层组织的收缩现象。鉴于此,本实验根据介质连续性定理,提出了一种多孔材料干燥组织收缩递推模型,并利用该模型研究了热风对流干燥香蕉的组织收缩,以期深入理解生物多孔材料干燥过程中的组织收缩机理。

1 数学模型

1.1 传质模型

对于无限长圆柱形生物多孔材料的干燥过程,假设物料内部水分的扩散阻力控制着整个干燥过程,则圆柱形物料内的水分扩散可由Fick定律描述为:

undefined

初始条件及边界条件为:

undefined

式中:m为圆柱形物料的局部干基含水率,kg water/kg dry matter;m0为圆柱形物料初始干基含水率,kg water/kg dry matter;me为圆柱形物料的平衡干基含水率,kg water/kg dry matter;t为时间,s;r为某点到圆柱形物料中心轴线的距离,m;Deff为生物多孔材料的水分扩散系数,m2/s。

圆柱形物料内部各点干基含水率随时间的变化可采用有限差分法由式(1)离散求得。

1.2 组织收缩递推模型

为了方便计算与分析问题,采用外节点法[8]将圆柱形物料颗粒从圆柱中心轴线至圆柱面沿径向均匀划分为N个节点,这样可将圆柱形物料颗粒离散为N个单元:第1单元为半径r=Δr/2的圆柱,第2单元至第N-1单元为厚度为Δr的圆筒,第N单元亦即最外一层单元,为厚度为Δr/2的圆筒。

对于圆柱形的生物多孔材料香蕉,由实验[9]得到:

undefined

式中:V0、V分别为圆柱形生物材料香蕉的初始体积和第n时间步的体积,m3;m0、undefined(n)分别为圆柱形物料香蕉的初始干基含水率和第n时间步的干基平均含水率,kg water/kg dry matter。

对于无限长圆柱形香蕉,忽略轴向收缩,式(3)可以修正为:

undefined

式中:r(n)为第n时间步圆柱形生物材料香蕉的外表面半径,m;R0为圆柱形生物材料香蕉的初始外表面半径,m。

对于无限长圆柱物料,忽略轴向收缩,将前面圆柱形生物多孔材料香蕉体积随干基平均含水率收缩变化规律(式(3))应用于每一单元格或单元体,可以得到:

undefined

式中:r(0,i)、r(0,i-1)分别为第i单元初始时刻外边界和内边界半径,m;r(n,i)、r(n,i-1)分别为第i单元第n时间步外边界和内边界半径,m;m(n,i)为第i单元第n时间步的局部干基含水率,kg water/kg dry matter。

2 一致性证明

由式(4)可以得到圆柱形生物材料香蕉经过第n时间步干燥后从最外层至最内层各单元层的外表面半径。

第N单元(最外层):

undefined

第N-1单元(次外层):

undefined

第2单元:

undefined

第1单元(最内层):

undefined

将以上各式相加,可得:

undefined

由于:

undefined

undefined

所以可以得到式(6):

undefined

式中:h为圆柱形生物材料香蕉的高度,m;ρ为圆柱形生物材料香蕉的初始密度,kg/m3。

S1=rundefined

S2=rundefined-rundefined

SN-1=rundefined-rundefined

SN=rundefined-rundefined

比较式(4)和式(6),由组织收缩递推模型基于圆柱形生物材料香蕉局部干基含水率计算得到的香蕉外径收缩尺寸与由颗粒整体体积收缩模型基于平均干基含水率计算得到的香蕉外径收缩尺寸是一致的。由此可见,利用香蕉体积随平均干基含水率的一维线性收缩变化规律研究香蕉各单元体厚度随单元体局部干基含水率的变化是可行的。

3 计算结果分析

对于初始干基含水率m0=2.914kg water/kg dry matter、直径D=40mm的香蕉,采用由Baini等[9]实验测得的有效扩散系数(Deff=1.0910-9m2/s),沿半径方向划分为11个单元节点,时间步长Δt=600s,利用式(1)、(2),采取有限差分方法,计算获得各单元节点随干燥时间变化的局部干基含水率。然后,将圆柱状生物多孔材料香蕉体积随其平均干基含水率收缩的变化规律应用于每一单元格或单元体,由式(5)可得到各单元外边界半径、外边界位移及厚度变化量随时间的变化。

图1为香蕉各节点的局部干基含水率随干燥时间的变化曲线。从图1中可以看到:①在干燥过程中,香蕉内部各点的局部干基含水率呈梯度分布,外层各点的局部干基含水率小于内层各点;②香蕉内部水分向外扩散是逐步进行的,其主要驱动力为干基含水率梯度。由于香蕉表面水分蒸发,香蕉外层与表面形成干基含水率梯度分布。在干基含水率梯度的推动下,香蕉外层水分向表面扩散,而随着外层水分的扩散,外层与次外层局部干基含水率梯度产生,次外层的水分在干基含水率梯度的驱动下向外层迁移,如同外层单元和次外层单元,随着香蕉内部距离表面较近的各单元水分向外扩散,该单元与其近邻内层单元也形成了干基含水率梯度,进而推动近邻内层单元内水分向外层扩散,在图1中表现为内层单元局部干基含水率的变化滞后于外层单元干基含水率的变化。同时还应注意到,随着干燥时间的延长,尤其是到了干燥的后期,香蕉内部水分大量散失,各单元层之间的干基含水率梯度逐渐减小直至为零,干燥过程结束。

图2和图3分别为香蕉各单元外边界半径与各单元外边界位移随干燥时间的变化曲线。

从图2和图3可以看到,香蕉最外层单元外边界半径首先收缩,并且在整个干燥过程中其外边界位移(外边界半径收缩幅度)最大。其它内层各单元外边界半径依次滞后收缩,且外边界位移依次减小。这是因为在干燥时香蕉所含水分是由外层至内层依次扩散蒸发离开香蕉的,香蕉外层的水分首先扩散,内层的水分滞后于外层逐渐扩散离开。

图4和图5分别为干燥过程中圆柱形香蕉各单元厚度变化量与各单元层厚度随干燥时间的变化曲线。图4表明,外层单元厚度的变化量由0快速增至极大值,然后快速回落,而内层单元厚度的变化量相对于外层单元,其由0到极大值的增加及回落速度都要缓慢很多。原因在于干燥过程中,外层单元之间干基含水率梯度大于内层,使外层水分扩散速度大于内层(图1)。此外,从图4还可以看到,外层单元厚度变化量的极大值大于内层,与外层单元体积大于内层相关联。

从图5可以看出,香蕉最外层单元在干燥开始阶段,其厚度从初始厚度降到极小值,然后又缓慢增加到某一稳定值。这一方面是由定值边界条件引起的,另一方面是由筒形物体的属性决定的,筒形物体体积保持不变,其内径收缩,则筒形物体厚度变大。另外,从图5中还可以看到,除了最外层单元外,其它各单元随着单元体内水分的扩散迁移,单元体厚度逐渐收缩变薄,其中距离香蕉外表面较近的几个单元的厚度随着干燥时间的延长达到极小值,之后又随着干燥时间的延长缓慢增加。这是由于香蕉各离散单元在干燥过程中,除了轴心第1单元体外,其余各单元体厚度的变化同时受2个因素的影响:①本单元体内水分的扩散引起的单元体体积收缩;②内层单元体体积收缩引起的本单元体内边界半径的收缩。图5中外层几个单元随着干燥时间的延长单元厚度变薄,主要是由于本单元体内水分的扩散导致本单元体收缩起主要作用;单元体厚度达到极小值后又变厚,主要原因是此时单元体干基含水率已经很低,随着干燥时间的延长,水分扩散引起的单元体体积变化很小,此时内层单元由于水分扩散导致体积收缩所引起的本单元内边界半径收缩起主导作用。

4 结论

(1)基于连续介质理论,建立了生物多孔材料干燥组织非稳态收缩递推模型,该模型可以很好地预测圆柱状生物多孔材料香蕉的组织收缩规律。

(2)将干燥过程中物料整体体积收缩随物料平均干基含水率变化的一维线性关系应用于物料各离散单元,计算分析了生物多孔材料干燥过程中组织的非稳态收缩规律,并从理论上证明了将整体体积随物料平均干基含水率变化的一维线性关系应用于物料各离散单元所采用的递推模型与整体收缩模型时,物料外边界尺寸收缩是一致的。

(3)计算结果表明,物料颗粒在干燥过程中随着水分的对外扩散和介质干基含水率降低,收缩并不是由外到内同步收缩的,而是先失去水分的外部先收缩,失水多的外部收缩幅度大于失水少的内部,没有失水的部位不收缩。

摘要：基于连续介质理论,建立了多孔材料干燥组织非稳态收缩递推模型。将物料整体体积随平均干基含水率变化的一维线性关系应用于物料各离散单元所采用的递推模型与整体收缩模型,证明了物料外边界尺寸收缩的一致性。此外,还利用此模型计算了热风对流干燥条件下香蕉组织随干燥时间的变化。结果表明,香蕉在干燥过程中随着水分的对外迁移,水分减少,香蕉组织收缩不是由外至内同步收缩,而是先失水的外部先收缩,不失水的部位不收缩。

关键词：香蕉,非稳态收缩,递推模型,对流干燥

参考文献

[1] Saravacos G D,Charm S E.A study of the mechanism offruit and vegetable dehydration[J].Food Techn,1962,16(1):78

[2] Perez M G R,Calvelo A.Modelling the thermal conductivi-ty of cooked meat[J].J Food Sci,1984,49(1):152

[3] Talla A,Puiggali J R,Jomaa W,et al.Shrinkage and den-sity evolution during of tropical fruits:Application to bana-na[J].J Food Eng,2004,64(1):103

[4]Zhang Jingping(张京平),et al.Derivation of boundaryshrinkage formula for round fruit during drying process(干燥过程中球形果蔬边界的收缩方程)[J].Trans ChineseSoc Agricult Machinery(农业机械学报),2005,36(9):71

[5]Zhang Jingping(张京平),et al.Preliminary study on theshrinkage equation of boundary of spheral fruits and vegeta-bles(干燥过程中球形果蔬边界的收缩方程)[J].TransChinese Soc Agricult Eng(农业工程学报),2002,18(6):141

[6] Biaobrzewski I.Simulation of changes in the density of anapple slab during drying[J].Int Commun Heat Mass Trans-fer,2006,33(7):880

[7] Biaobrzewski I,Zielińska M,Mujumdar A S,et al.Heatand mass transfer during drying of a bed of shrinkage parti-cles simulation for carrot cubes dried in a spout-fluidized-beddrier[J].Int J Heat Mass Transfer,2008,51(19-20):4704

[8]陶文铨.数值传热学[M].西安:西安交通大学出版社,1988

动力机械非稳态过程的计算方法探索 第5篇

对于动力机械而言, 在非稳定工况情形下, 如突然启动、断电停机、快速调节工况、快速变频等, 描述它性能的每一个参数在短时间内都将发生剧烈的变化。一般情形下, 为了减少研究强度, 通常仅对研究对象开展研究。若对单一研究对象进行内部流动计算, 在指定边界条件时, 则必须知道转速-时间曲线和流量-时间曲线, 或者转速-时间曲线和进出口压力-时间曲线, 而流量与压力的变化又与管路系统特性密切相关。因此, 要准确模拟出内部流场, 需要通过实验确定上述关系。在各式各样的瞬态操作过程中, 若存在转速变化, 则可通过基于用户自定义函数功能来实现, 内部旋转区域作为一个整体旋转。旋转的转子和静止的定子在计算过程中存在相对运动, 两个区域通过滑移面联接。对于瞬态过程中的非定常计算, 每个时间步长计算一次交迭面;对于每个滑移面上的通量按照重叠面的面积比例计算得到, 从而实现两个区域的通量平衡。下面以一台低比转速离心泵的启动过程为例, 来探索数值模拟方法的可行性。

1 模型与方法

1.1 计算模型

本文中的计算模型与参考文献中的计算模型完全相同[1], 模型泵额定流量Q=6 m3/h, 扬程H=8 m, 转速n=1 450 r/min。计算域如图1所示, 计算域分为3部分, 分别是进口段, 叶轮旋转区域和蜗壳静止区域。叶轮旋转区域作为一个整体旋转。

1.2 计算方法

在快速启动过程中, 转速-时间曲线基本上取决于动力源的启动特性, 并且该曲线关系近似满足于指数关系[2]:

其中:nmax=1 450 r/min;t0=0.15 s。转速变化关系通过编写用户自定义函数施加到叶轮上。

已有文献表明, 流量变化近似满足三次曲线形式。结合文献[3], 在此处采用人为给定的形式, 最大稳定流量设定为6 m3/h。假定的流量变化关系式如图2中的流量曲线所示。根据质量守恒和几何关系, 可换算成速度, 即在泵进口处采用速度进口边界条件, 在泵出口处采用自由出流条件。进口速度变化关系通过用户自定义函数施加到泵进口处。湍流模型采用RNG k-ε湍流模型封闭时均雷诺方程。计算中的滑移网格技术已被广泛用来计算稳定工况的性能预测[4]。瞬态项的离散采用一阶隐式格式, 压力和速度的耦合采用SIMPLE算法实现, 对流项采用一阶迎风格式离散, 扩散项采用中心差分格式离散。时间步长取0.000 1 s, 从静止到达稳定转速的时间约为1 s。在每个时间步长内取最大迭代次数为200次 (实际上每个时间步长内迭代不到200次便可收敛) , 以保证在每个时间步长内都绝对收敛。收敛残差为0.000 1。计算介质为清水, 密度ρ=1 000 kg/m3, 动力黏度μ=0.001 Pas, 考虑重力影响, 计算域网格总数为508 792。

2 计算结果

图2是计算得到的启动过程离心泵扬程变化的时间历程, 其中转速和流量曲线是外部指定的边界条件与初始条件。可以看到, 随着转速的增加, 扬程波动值也随之增加。最大波动值约为1.5 m, 相对于设计扬程, 可见波动时非常大的。

图3是启动过程中离心泵中截面上总压变化的时间历程。在0.04 s, 由于转速相对较低, 故叶轮搅动水体做功, 使得此时内部压力场变化并不明显。随后随着转速的不断上升, 压力梯度变得更为明显。在任意时刻, 压力逐渐演变为从进口处的低压转变为叶轮出口处的高压。并且中心压力随着转速的增加而迅速降低。在0.6 s之后, 叶轮部分压力变化并不明显, 变化主要发生在蜗壳内部。从压力场的演化趋势与扬程外特性曲线的变化趋势结合分析可知, 两者具有较好的对应关系。

3 结论

此算例采用滑移网格技术及用户自定义函数相结合的方式对单一泵进行了启动过程的三维非定常流动数值预测, 从内流场和外特性的结果看, 数值模拟方法是可行的。但由于转速和流量的变化历程是人为给定的, 尤其是流量条件的给定更具不确定性, 因此无法获得准确可靠的非定常流场演化结果。根据前述分析, 要获得准确的流量、压力曲线需要借助于实验结果, 才能数值模拟得到准确的内流场演化特征。

参考文献

[1]Zhang Yuliang, Li Yi, Cui Baoling, et al.Numerical simulation and analysis of solid-liquid two-phase flow in centrifugal pump[J].Chinese Journal of Mechanical Engineering, 2013, 26 (1) :53-60.

[2]Tsukamoto H, Ohashi H.Transient characteristics of a centrifugal pump during startup period[J].ASME Journal of Fluids Engineering, 1982, 104 (1) :6-13.

[3]平仕良, 吴大转, 王乐勤.离心式水泵快速开启过程的瞬态效应分析[J].浙江大学学报:工学版, 2007, 41 (5) :814-817.

非稳态特性 第6篇

近几年,随着大区电网互联,实际系统中发生了一些由周期性负荷或扰动引发的功率振荡现象,这类低频振荡难以采用传统的负阻尼理论进行解释,因而引起了越来越多研究人员的关注,并开展了电力系统强迫功率振荡方面的研究[1,2,3,4,5,6,7,8,9,10]。这些研究表明,持续的周期性小扰动会引发电力系统强迫振荡,当扰动频率接近系统固有频率时,会引起系统共振,导致大幅度的功率振荡,或称为共振机理的低频振荡[3,4,5]。它具有起振快、起振后保持等幅同步振荡和失去振荡源后振荡很快衰减等特点。文献[3]计及周期性扰动的多机系统,通过对非齐次线性方程组解的推导分析,提出了当扰动频率接近系统固有振荡频率时将发生大幅度强迫功率振荡的结论。文献[4]基于单机无穷大系统分析,提出了电力系统强迫功率振荡的基础理论,并分析了其主要的影响因素。文献[5]通过对河北南网安保线低频功率振荡的分析,认为它属于强迫振荡,并通过仿真指出强迫振荡可由发电机的轴系、励磁器、调速器之一的周期扰动引发。文献[6]从能量的角度分析了电力系统强迫振荡共振稳态时内外能量的转换关系和特征。在强迫振荡扰动源方面,文献[7]基于电力系统线性化模型的频率响应,研究了弱阻尼情况下,低速柴油发电机引起的强迫振荡问题及其影响因素。文献[8]发现汽轮机热力系统中存在激发电力系统强迫振荡的扰动源。文献[9]研究了汽轮机功率变化的原因,发现调节汽门扰动频率与电力系统固有振荡频率一致或接近时,均可能引起电力系统发生共振机理的低频振荡。文献[10]从机理上研究了原动机功率与负荷持续周期性小扰动造成电网功率振荡的区别,分析说明了两者的不同性质。

目前,传统的负阻尼低频振荡分析方法已经比较成熟,通过对电力系统线性化方程的模态分析可以获得系统的全部振荡模式,对特征向量的分析,又可以获得参与各振荡模式的强相关机组信息[11,12,13]。而电力系统强迫功率振荡方面,已有的研究很多都是基于单机无穷大系统下的分析,在多机系统情况下,强迫功率振荡发生共振的条件、振荡分布特点与单机无穷大系统又有很大的不同,目前这方面研究还有待于进一步深入。

本文采用复模态叠加方法对多机系统强迫功率振荡进行分析,研究其发生共振的条件及其影响振荡幅度的主要因素。通过比较弱阻尼情况下多机系统强迫功率振荡稳态响应与传统的负阻尼低频振荡响应之间的相似点和不同点,有助于加深对多机系统强迫功率振荡的认识和理解。

1 多机系统复模态模型

在包含n台发电机的多机系统中,由于主要关注振荡频率在0.1 Hz～2.5 Hz之间的机电模式,而它们主要由转子运动方程决定,所以可以只考虑转子摇摆方程[14]。发电机采用经典二阶模型(假定E′恒定),其线性化的状态方程为:

$\dot{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}\mathit{u}(1)$

为了消除方程间的耦合,采用模态坐标,可令x=Φz,代入可得系统新的坐标方程为:

$\dot{\mathit{z}}=\mathit{\Phi }{}^{-1}\mathit{A}\mathit{\Phi }\mathit{z}+\mathit{\Phi }{}^{-1}\mathit{B}\mathit{u}=\mathit{\Lambda }\mathit{z}+\mathit{\Phi }{}^{-1}\mathit{B}\mathit{u}(2)$

式中:Φ为A的右特征向量矩阵。

由于状态方程只计入转子摇摆方程,当n机系统正常稳定运行时,通常含有反映机电振荡模式的n-1对共轭复特征根,其相应特征向量也是共轭复数向量,右特征向量矩阵Φ=[Φ1,Φ2,,Φn-1,Φ*1,Φ*2,Φ*n-1],左特征向量矩阵Ψ=Φ-1=[ΨT1,ΨT2,,ΨTn-1,Ψ*T1,Ψ*T2,,Ψ*Tn-1]T,所以采用复模态方法[15],模态坐标也可以表示为共轭复数的形式,即z=[z1,z2,,zn-1,z*1,z*2,,z*n-1]T。此时系统时域响应可以表示为n-1个机电振荡模式的叠加形式,即$\mathit{x}(t)={\displaystyle \sum _{r=1}^{n-1}(}\mathit{\Phi }{}_{r}z{}_r+\mathit{\Phi }{}_r^{*}z{}_r^{*})$,而每个振荡模式由2个独立的共轭模态坐标组成。

当采用传统方法分析负阻尼低频振荡时,一般忽略机械功率变化。此时,系统第r阶振荡模式解耦方程为:

$\{\begin{array}{l}\dot{z}{}_r=\lambda {}_{r}z{}_r\\ \dot{z}{}_r^{*}=\lambda {}_r^{*}z{}_r^{*}\end{array}(3)$

系统负阻尼自由振荡为n-1对共轭模态响应的叠加,其时域解为:

$\mathit{x}(t)={\displaystyle \sum _{r=1}^{n-1}(}\mathit{\Phi }{}_r\mathit{\Psi }{}_r^{\text{Τ}}\mathit{x}{}_0\text{e}{}^{\lambda {}_{r}t}+\mathit{\Phi }{}_r^{*}\mathit{\Psi }{}_r^{*\text{Τ}}\mathit{x}{}_0\text{e}{}^{\lambda {}_r^{*}t})(4)$

第i个状态量可表示为:

$x{}_i(t)={\displaystyle \sum _{r=1}^{n-1}(}\phi {}_{ir}z{}_{r0}\text{e}{}^{\lambda {}_{r}t}+\phi {}_{ir}^{*}z{}_{r0}^{*}\text{e}{}^{\lambda {}_r^{*}t})(5)$

上式各值均为复数,可表示为φir=|φir|ejγir,zr0=|zr0|ejθr,λr=-αr+jωdr,则时域解为:

$\begin{array}{l}x{}_i(t)={\displaystyle \sum _{r=1}^{n-1}|}\phi {}_{ir}||z{}_{r0}|\text{e}{}^{-\alpha {}_{r}t}(\text{e}{}^{\text{j}(\omega {}_{dr}t+\gamma {}_{ir}+\theta {}_r)}+\\ \text{e}{}^{-\text{j}(\omega {}_{dr}t+\gamma {}_{ir}+\theta {}_r)})=2{\displaystyle \sum _{r=1}^{n-1}|}\phi {}_{ir}||z{}_{r0}|\cdot \\ \text{e}{}^{-\alpha {}_{r}t}\cos (\omega {}_{\text{d}r}t+\gamma {}_{ir}+\theta {}_r)(6)\end{array}$

此时,系统表现为零输入响应。当系统阻尼为负时,系统将增幅振荡;当阻尼为0时,系统将维持等幅振荡;当阻尼为正时,系统振荡将逐渐衰减。

2 多机系统强迫功率振荡

下面分析多机系统强迫功率振荡的情况。采用模态坐标后,系统解耦方程可以写为:

$\left[\begin{array}{c}\dot{\mathit{z}}\\ \dot{\mathit{z}}{}^{*}\end{array}\right]=\left[\begin{array}{cc}\mathit{\Lambda }& \\ & \mathit{\Lambda }{}^{*}\end{array}\right]\left[\begin{array}{c}\mathit{z}\\ \mathit{z}{}^{*}\end{array}\right]+\mathit{\Phi }{}^{-1}\text{Δ}\mathit{{\rm P}}{}_{\text{Τ}}(7)$

式中:ΔPT=ΔPTmsin ωt为机械功率扰动。

因此,系统第r阶振荡模式解耦方程为:

$\{\begin{array}{l}\dot{z}{}_r=\lambda {}_{r}z{}_r+\mathit{\Psi }{}_r\text{Δ}\mathit{{\rm P}}{}_{\text{Τ}}\\ \dot{z}{}_r^{*}=\lambda {}_r^{*}z{}_r^{*}+\mathit{\Psi }{}_r^{*}\text{Δ}\mathit{{\rm P}}{}_{\text{Τ}}\end{array}(8)$

其解由通解和特解2部分组成,由于系统阻尼的存在,更关注自由振荡消失后,系统强迫振荡的稳态响应,此时系统各变量都以正弦形式变化,因此可采用复相量法[15]求解,即取$\text{Δ}\mathit{{\rm P}}{}_{\text{Τ}}=\text{Δ}\widetilde{\mathit{{\rm P}}}{}_{\text{Τ}\text{m}}\text{e}{}^{\text{j}\omega t}$参与运算,获得结果$\mathit{x}(t)=\tilde{\mathit{B}}\text{e}{}^{\text{j}\omega t}$,然后取虚部作为最终时域响应,其中Δ$\widetilde{\mathit{{\rm P}}}$Tm和$\tilde{\mathit{B}}$均为复相量。下面公式中的变量也表示为复相量形式。则解耦方程的解为:

$\{\begin{array}{l}z{}_r=\frac{\tilde{\mathit{\Psi }}{}_r\text{Δ}\widetilde{\mathit{{\rm P}}}{}_{\text{Τ}\text{m}}}{\text{j}\omega -\tilde{\lambda }{}_r}\text{e}{}^{\text{j}\omega t}\\ z{}_r^{*}=\frac{\tilde{\mathit{\Psi }}{}_r^{*}\text{Δ}\widetilde{\mathit{{\rm P}}}{}_{\text{Τ}\text{m}}}{\text{j}\omega -\tilde{\lambda }{}_r^{*}}\text{e}{}^{\text{j}\omega t}\end{array}(9)$

系统强迫功率振荡稳态解也为n-1对共轭模态响应的叠加,其解为:

$\begin{array}{l}\mathit{x}(t)={\displaystyle \sum _{r=1}^{n-1}\tilde{\mathit{\Phi }}}{}_{r}z{}_r+\tilde{\mathit{\Phi }}{}_r^{*}z{}_r^{*}=\\ {\displaystyle \sum _{r=1}^{n-1}[}\frac{\text{j}\omega (\tilde{\mathit{\Phi }}{}_r\tilde{\mathit{\Psi }}{}_r+\tilde{\mathit{\Phi }}{}_r^{*}\tilde{\mathit{\Psi }}{}_r^{*})-}{(\text{j}\omega -\tilde{\lambda }{}_r)\cdot }\\ \leftarrow \frac{(\tilde{\mathit{\Phi }}{}_r\tilde{\mathit{\Psi }}{}_r\tilde{\lambda }{}_r^{*}+\tilde{\mathit{\Phi }}{}_r^{*}\tilde{\mathit{\Psi }}{}_r^{*}\tilde{\lambda }{}_r)}{(\text{j}\omega -\tilde{\lambda }{}_r^{*})}\text{Δ}\widetilde{\mathit{{\rm P}}}{}_{\text{Τ}\text{m}}\text{e}{}^{\text{j}\omega t}](10)\end{array}$

定义ωnr=|λr|为第r阶振荡模式的固有频率,阻尼比ζr=αr/ωnr,阻尼振荡频率$\omega {}_{\text{d}r}=\sqrt{1-\zeta {}_r^2}\omega {}_{nr}$,取$\mathit{A}{}_r=-(\tilde{\mathit{\Phi }}{}_r\tilde{\mathit{\Psi }}{}_r\tilde{\lambda }{}_r^{*}+\tilde{\mathit{\Phi }}{}_r^{*}\tilde{\mathit{\Psi }}{}_r^{*}\tilde{\lambda }{}_r),\mathit{B}{}_r=(\tilde{\mathit{\Phi }}{}_r\tilde{\mathit{\Psi }}{}_r+\tilde{\mathit{\Phi }}{}_r^{*}\tilde{\mathit{\Psi }}{}_r^{*})$,则式(10)可以整理为:

$\mathit{x}(t)={\displaystyle \sum _{r=1}^{n-1}\left(\frac{\mathit{A}{}_r+\text{j}\omega \mathit{B}{}_r}{\omega {}_{nr}^2-\omega {}^2+\text{j}2\zeta {}_r\omega {}_{nr}\omega }\text{Δ}\widetilde{\mathit{{\rm P}}}{}_{\text{Τ}\text{m}}\text{e}{}^{\text{j}\omega t}\right)}(11)$

为了方便分析,假设系统中某一台机组上存在简谐机械功率扰动,相位为0,其在Δ$\widetilde{\mathit{{\rm P}}}$Tm中对应元素为$\tilde{p}$l,则$\tilde{\mathit{\Psi }}{}_r\text{Δ}\widetilde{\mathit{{\rm P}}}{}_{\text{Τ}\text{m}}=\tilde{\psi }{}_{rl}\tilde{p}{}_l,\tilde{\mathit{\Psi }}{}_r^{*}\text{Δ}\widetilde{\mathit{{\rm P}}}{}_{\text{Τ}\text{m}}=\tilde{\psi }{}_{rl}^{*}\tilde{p}{}_l$。

则第i个状态变量的解为:

$\begin{array}{l}x{}_i(t)={\displaystyle \sum _{r=1}^{n-1}[}\frac{\text{j}\omega (\tilde{\phi }{}_{ir}\tilde{\psi }{}_{rl}+\tilde{\phi }{}_{ir}^{*}\tilde{\psi }{}_{rl}^{*})-}{\omega {}_{\text{n}r}^2-\omega {}^2+}\\ \leftarrow \frac{(\tilde{\lambda }{}_r^{*}\tilde{\phi }{}_{ir}\tilde{\psi }{}_{rl}+\tilde{\lambda }{}_r\tilde{\phi }{}_{ir}^{*}\tilde{\psi }{}_{rl}^{*})}{\text{j}2\zeta {}_r\omega {}_{\text{n}r}\omega }\tilde{p}{}_l\text{e}{}^{\text{j}\omega t}]=\end{array}$

${\displaystyle \sum _{r=1}^{n-1}}$$\frac{\frac{1}{\omega {}_{\text{n}r}^2}(a+\text{j}\omega b)}{1-\left(\frac{\omega }{\omega {}_{\text{n}r}}\right){}^2+\text{j}2\zeta {}_r\frac{\omega }{\omega {}_{\text{n}r}}}\tilde{p}{}_l\text{e}{}^{\text{j}\omega t}](12)$

式中:$a=-(\tilde{\lambda }{}_r^{*}\tilde{\phi }{}_{ir}\tilde{\psi }{}_{rl}+\tilde{\lambda }{}_r\tilde{\phi }{}_{ir}^{*}\tilde{\psi }{}_{rl}^{*}),b=\tilde{\phi }{}_{ir}\tilde{\psi }{}_{rl}+\tilde{\phi }{}_{ir}^{*}\tilde{\psi }{}_{rl}^{*}$都为实数。

取v=ω/ωnr为频率比,则第r阶振荡模式时域响应为:

$x{}_i^r(t)=\text{Ι}\text{m}[\frac{\frac{a}{\omega {}_{\text{n}r}^2}+\text{j}\frac{v{}_{r}b}{\omega {}_{\text{n}r}}}{1-v{}_r^2+\text{j}2\zeta {}_{r}v{}_r}\tilde{p}{}_l\text{e}{}^{\text{j}\omega t}〗(13)$

取x${}_i^r$(t)=B${}_{il}^r$sin(ωt-ϕ${}_{il}^r$),则

$\begin{array}{l}B{}_{il}^r=\sqrt{\frac{\left(\frac{a}{\omega {}_{\text{n}r}^2}\right){}^2+\left(\frac{v{}_{r}b}{\omega {}_{\text{n}r}}\right){}^2}{(1-v{}_r^2){}^2+(2\zeta {}_{r}v{}_r){}^2}}|\tilde{p}{}_l|(14)\\ \text{ϕ}{}_{il}^r=\arctan \frac{2\zeta {}_{r}v{}_{r}a-(1-v{}_r^2)v{}_{r}b\omega {}_{\text{n}r}}{a(1-v{}_r^2)+2\zeta {}_{r}v{}_r^{2}b\omega {}_{\text{n}r}}(15)\end{array}$

对式(14)求导,可得存在最大幅值的频率比为:

$v{}_r^2=\sqrt{1+2(1-2\zeta {}_r^2)\left(\frac{a}{b\omega {}_{\text{n}r}}\right){}^2+\left(\frac{a}{b\omega {}_{\text{n}r}}\right){}^4}-\left(\frac{a}{b\omega {}_{\text{n}r}}\right){}^2(16)$

由式(16)可见,系统共振频率总落在下列频带范围内:1-2ζ${}_r^2$v${}_r^2$1。当该模式阻尼比ζr比较小时,vr≈1,即可以认为系统的共振频率等于模式固有振荡频率。当扰动频率ω与系统第r阶模式的固有频率ωnr相近时,第r阶振荡模式的稳态响应振幅比其他阶模式响应大得多,可以近似将其认为是系统总响应。为了便于分析,将复特征值和复特征向量元素表示为$\tilde{\phi }{}_{ir}=|\tilde{\phi }{}_{ir}|\text{e}{}^{\text{j}\gamma {}_{ir}},\tilde{\psi }{}_{rl}=|\tilde{\psi }{}_{rl}|\text{e}{}^{\text{j}\sigma {}_{rl}},\tilde{\lambda }{}_r=\omega {}_{\text{n}r}\text{e}{}^{\text{j}\text{ϕ}{}_r}$形式,则第i个状态变量时域解为:

xi(t)≈Im$\frac{|\tilde{\phi }{}_{ir}||\tilde{\psi }{}_{rl}||\tilde{p}{}_l|}{2\zeta {}_r\omega {}_{\text{n}r}}$

$\text{e}{}^{\text{j}(\gamma {}_{ir}+\sigma {}_{rl})}+\text{e}{}^{-\text{j}(\gamma {}_{ir}+\sigma {}_{rl})}+\text{e}{}^{\text{j}(\gamma {}_{ir}+\sigma {}_{rl}-\text{ϕ}{}_r+\frac{\text{π}}{2})}+$

$\text{e}{}^{-\text{j}(\gamma {}_{ir}+\sigma {}_{rl}-\text{ϕ}{}_r-\frac{\text{π}}{2})}$ejωnrt (17)

由于多机系统分析中,主要关注弱阻尼模式,其阻尼比ζ较小,即特征值实部远小于虚部,ϕr≈π/2。此时,式(17)可以整理为:

$x{}_i(t)\approx \frac{|\tilde{\phi }{}_{ir}||\tilde{\psi }{}_{rl}||\tilde{p}{}_l|}{\zeta {}_r\omega {}_{\text{n}r}}\sin (\omega {}_{\text{n}r}t+\gamma {}_{ir}+\sigma {}_{rl})(18)$

由此得到了弱阻尼情况下,多机系统强迫功率振荡共振时的状态变量稳态时域响应。

3 多机系统强迫功率振荡稳态响应特性

多机系统强迫功率振荡稳态响应特性如下:

1)多机系统强迫功率振荡稳态响应可以看做是各阶振荡模式稳态响应的叠加。假设第r阶振荡模式为弱阻尼的情况下,近似认为当扰动频率ω接近该振荡模式固有频率ωnr时,系统发生共振,此时系统响应可以用第r阶振荡模式稳态响应近似代替。当系统存在多个弱阻尼模式时,则存在多个可能引发系统大幅度共振的扰动频率。

2)在系统共振频带范围内,强迫功率振荡稳态响应振幅不仅与扰动幅值、系统阻尼大小有关,而且与振荡源所在机组参与该阶振荡模式的程度有很大关系,若机组参与该模式较小,即使发生共振,振幅也不大。因此往往更加关心系统阻尼比较小、扰动源机组参与程度较高的模式。

3)多机系统负阻尼低频振荡与强迫功率振荡稳态响应表现形式不同,负阻尼振荡响应由n-1个不同阻尼振荡频率ωdr的模式叠加,每一模式响应含有不同的e-αrt衰减项,一般通过Prony方法分析实际振荡曲线时,可以获得多个振荡频率及阻尼比等信息。而强迫功率振荡稳态响应由多个同频率等幅振荡的振荡模式稳态响应叠加而成,通过Prony方法只能获得单一振荡频率,难以获得阻尼比等具体单个模式信息。

4)多机系统负阻尼低频振荡,某个振荡模式下系统内不同机组Δδi之间的幅值比和相位差是固定的,可以由右特征向量元素衡量,与系统初始条件和时间t无关。而多机系统强迫功率振荡稳态响应在弱阻尼共振情况下,不同机组Δδi之间的相对振幅大小和相位差同样可以近似采用右特征向量元素衡量,与振荡源所在位置及时间t无关。这一点可以看出两者的振荡分布是相似的。

5)多机系统强迫功率振荡弱阻尼共振情况下,各状态量相位呈现出一种较为复杂的关系,没有显现单机无穷大系统共振情况下转子角响应相对于机械功率扰动存在π/2相位滞后的简单特点,但是各机组Δδi的相位差满足右特征向量分布。

4 算例分析

采用2区4机系统,详细参数见文献[16],系统结构、振荡模式及转子相关右特征向量见附录A。发电机采用5阶模型,负荷以恒阻抗表示。此时系统共有3个振荡模式,固有振荡频率分别为0.578 3 Hz,1.164 5 Hz和1.196 4 Hz。阻尼比ζ分别为2.95%,10.54%和10.51%。其中,模式1为机组G1和G2相对机组G3和G4的区间振荡模式,模式2为G1与G2间的局部振荡模式,模式3为G3与G4间的局部振荡模式。

假设机械功率扰动位于G1上,扰动幅值为0.03(标幺值),此时系统幅频特性如图1所示。

当扰动频率接近模式1固有频率0.578 3 Hz时,系统发生共振,由于该模式系统阻尼比较小,且4台机组都参与该模式,所以振荡幅值都很大,G3参与程度最大,振荡幅值也最大。当扰动频率接近模式2固有频率1.164 5 Hz时,系统也会发生共振,该模式系统阻尼比较大,共振振幅比起模式1小很多。模式2的G1和G2参与程度大,所以G1和G2振幅相对较大,而G3和G4几乎没有参与该振荡模式,振幅变化不明显。由于G1对模式3的参与程度很小,所以当G1存在强迫振荡扰动源时,难以激发模式3的共振。

图2给出了扰动源位于G1时系统的相频特性,当扰动频率ω接近于系统模式1和模式2固有频率时,扰动源所在G1的Δδ相位接近滞后外施扰动π/2,各机组间Δδ相位差与该共振模式右特征向量各对应元素的相位差基本一致。

当G2存在强迫机械功率扰动时,激发的共振模式、共振时各机组Δδ之间振幅比和相位差与G1存在扰动源时比较相似。

图3和图4给出了G3存在强迫振荡扰动源时的幅频特性和相频特性曲线。

随着扰动频率的变化,存在模式1和模式3的共振,而由于G3参与模式2的程度很小,所以难以激发模式2的共振。同时,由于G3参与模式1的程度比G1大,所以相同的扰动下,模式1共振的振幅比图1中大一些。共振情况下,机组间Δδ幅值比和相位差按照发生共振模式的振荡振型分布。当G4存在强迫振荡扰动源时,与G3存在扰动源时情况比较相似,只是激发各模式的绝对振幅不同。

新英格兰10机39节点系统分析见附录B。从结果来看,在多机系统中,是否存在弱阻尼的振荡模式,以及扰动源所在机组参与该振荡模式的程度对于共振振幅大小的影响非常大。对于弱阻尼局部振荡模式,只有在参与该振荡模式的发电机本身受到接近该模式固有频率的强迫功率扰动时,才会引起较大的振荡。而对于弱阻尼区间振荡模式的系统,由于参与的机组较多,更容易发生强迫功率振荡共振的情况,参与程度接近的几台机组中无论哪一台存在区间模式固有频率的强迫扰动,都会激发大幅度的功率振荡,虽然振荡源位置不同,但系统振荡分布却非常相似。

5 结语

本文采用复模态叠加方法推导了多机系统强迫功率振荡的稳态响应特性,分析其发生共振的条件以及主要影响因素。分析表明多机系统在某一振荡模式阻尼比较小的情况下,可以近似认为参与机组扰动频率接近该模式固有振荡频率时,会发生大幅度共振。共振幅值的大小不仅与外施扰动大小、系统阻尼有关,而且与扰动源所在机组参与该共振振荡模式的程度有很大关系。电力系统强迫功率振荡由于存在具体的振荡源,与负阻尼低频振荡起因不同,但是其各状态量共振稳态响应幅值和相位仍然近似按照右特征向量分布,这一点与负阻尼低频振荡非常相似。本文的推导基于线性理论,而多机电力系统毕竟是一个非线性系统,当系统非线性特性表现比较强的时候,本文的结论存在局限性。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。

非稳态特性 第7篇

1 旅行车与原型车的稳态回转特性的对比分析

原型车是一款具有独特悬架结构的旅行轿车,其悬架系统不同于一般车型,前部为麦弗逊支撑杆摆臂式独立悬架,后部同样为独立悬架,但带有横置的钢板弹簧,使该车具有较高的稳定性和抗侧倾能力,图1为以该型车为原型车所开发的旅行车。

汽车稳态回转特性试验通过中性转向点的侧向加速度值an、不足转向度U和车箱侧倾度KQ三项指标进行评价计分。其中,稳态回转试验的中性转向点的侧向加速度值an为汽车操纵稳定性稳态回转特性中具有否决权的指标[1]。图2为旅行车的稳态回转试验曲线,表1为侧向加速度值an、不足转向度U和车箱侧倾度KQ三项指标的评分值[2]。

对比分析图表可知,原型车的中性转向点的侧向加速度值an、不足转向度U和车箱侧倾度KQ三项指标的评分值都合格;旅行车的不足转向度U和车箱侧倾度KQ两项指标的评分值都合格,但中性转向点的侧向加速度值an评分值不合格,由于其为操纵稳定性中具有否决权的指标,因此旅行车操纵稳定性评价不合格。

旅行车和原型车采用相同的悬架结构,但其稳态回转特性却明显不同,从图2转弯半径比与侧向加速度的关系曲线可知,原型车表现为明显不足转向特性,而新开发的旅行车的稳态回转试验曲线在车身侧向加速度为0.4~0.6g的范围内存在一个较为明显中性转向点,因此对新开发的旅行车的稳态回转试验数据进行拟合,稳态回转特性总体表现为具有过度转向趋势,中性转向点的侧向加速度为4.6 m/s2,评分值为56.6分。

2 旅行车稳态回转特性的改进试验情况

为了改善旅行车稳态回转特性,针对旅行车稳态回转特性出现中性转向点的侧向加速度过低的问题进行了一系列试验研究,并在一定程度上改善了旅行车的操纵稳定性。上述研究工作从试验角度分析了车辆各总成、部件等设计参数对车辆操纵稳定性的影响,为旅行车操纵稳定性的改进设计提供了一定的试验依据,这些前期研究工作主要包括以下内容[2]。

(1)增加前悬架横向稳定杆直径。增加横向稳定杆直径将导致前桥侧倾角刚度的变大,使汽车前、后桥侧倾角刚度比值变大,将增加汽车的不足转向趋势。试验结果表明:进一步增加横向稳定杆直径,并不能明显改善旅行车稳态回转特性。

(2)为旅行车的后悬架加装横向稳定杆。试验结果表明:加装横向稳定杆后旅行车的稳态回转特性进一步恶化,这也从试验角度验证了汽车前、后桥侧倾角刚度比值变小,将增加汽车过度转向趋势的结论。

(3)降低前轮气压和增高后轮气压,希望降低前轮侧偏刚度和增加后轮侧偏刚度。试验结果表明:轮胎气压的改变对稳态回转特性的改善不明显。通过轮胎专项试验可知,尽管旅行车轮胎气压不同(前200k Pa,后220k Pa),轮胎侧偏刚度变化不大。

(4)降低钢板弹簧刚度。降低钢板弹簧刚度将导致后桥侧倾角刚度的变小,使汽车前、后桥侧倾角刚度比值变大,将增加汽车的不足转向趋势,但会对车箱侧倾度评价指标产生不利影响,同时造成车身后部下沉,承载能力下降,离地间隙过小,通过性差。

(5)将后悬架的单摆臂结构改为平行四边形的双摆臂结构。试验结果表明:改进方案对旅行车的稳态回转特性产生不利影响,同时使后悬架结构复杂,因此该方案不可行。

(6)前轮采用大的负前束。试验结果表明:前轮采用负前束将改善旅行车稳态回转特性,但会加剧轮胎磨损,该方案不具备工程适用性。

(7)改变汽车轴荷分配,使前轴轴荷变大,后轴轴荷变小,有利于改善旅行车稳态回转特性,但该方案不具备工程适用性。

3 旅行车稳态回转特性的影响因素分析

旅行车和原型车的悬架结构型式相同,但旅行车稳态回转特性却存在中性转向点的侧向加速度过低的问题。对比分析旅行车和原型车在整车质心位置、悬架侧倾角刚度、转向轮附加转角等方面的差异,探讨这些因素可能对稳态回转特性产生的影响。

3.1 整车质心差异对稳态回转特性的影响

表2为旅行车和原型车的轴荷分配及质心高度。根据汽车操纵稳定性的两自由度模型[3],稳定性因子:当K>0时,汽车的稳态响应为不足转向;当K=0时,汽车的稳态响应为中性转向;当K<0时,汽车的稳态响应为过多转向。假设前后轮的轮胎侧偏特性相同(k1=k2),旅行车整车质心距前轴距离a接近质心距后轴距离b,K=0,因此旅行车的稳态回转特性是中性转向,而原型车整车质心距前轴距离a小于质心距后轴距离b,K>0,故原型车的稳态回转特性是不足转向。

与原型车相比,旅行车整车质心略有提高,假设两车整车质量相同,在同样的侧向加速度下,旅行车所受的惯性侧倾力矩比原型车大,车身侧倾角变大,使车辆的操纵稳定性变差。

综上所述,整车质心位置的改变使旅行车的稳态回转特性相对原型车变差。

3.2 前后悬架侧倾角刚度差异对稳态回转特性的影响

旅行车和原型车的前悬架结构型式相同,都是麦弗逊式悬架装置(带横向稳定杆),但螺旋弹簧刚度和自由长度不同。旅行车和原型车的后悬架结构型式相同,都是单摆臂独立悬架带双向作用筒式液力减振器,但钢板板簧的刚度和弧高不同。旅行车的螺旋弹簧刚度为16.5±0.6(N/mm),螺旋弹簧自由长度为460(mm),钢板弹簧刚度为44.4(N/mm),钢板弹簧弧高为185(mm);原型车的螺旋弹簧刚度为15.3±0.6(N/mm),螺旋弹簧自由长度为445(mm),钢板弹簧刚度为27.0(N/mm),钢板弹簧弧高为160(mm)。

从保证汽车具有不足转向特性的目的出发,要求前轴的侧偏角大于后轴的侧偏角,从而要求前悬架的侧倾角刚度大于后悬架的侧倾角刚度。由于旅行车前悬架螺旋弹簧刚度比原型车只是略有提高,而钢板弹簧刚度却比原型车大幅度提高,因此,如果忽略其他相关因素的影响,可以认为由于旅行车前悬架螺旋弹簧刚度和钢板弹簧刚度的改变,使得前后悬架侧倾角刚度变小,不利于操纵稳定性的改善,有增加旅行车过度转向趋势。

3.3 侧倾时转向系与前悬架运动干涉对操纵稳定性的影响

汽车前悬架和转向系杆系干涉引起的附加转角主要表现在两方面[4]:一方面转向轮依靠前悬架导向杆系来保证车轮与车身的相对运动关系;另一方面,转向机与转向轮之间的传动杆件也使车轮与车身之间有一定的相对运动关系。这两种运动关系是由两套独立的机构来实现的,因此,要想做到完全一致几乎是不可能的,这种不一致的效果就造成车身侧倾和垂直位移时的运动干涉引起转向轮的附加转角,改变汽车的稳态回转特性,这种干涉既可能表现为增加汽车的不足转向量,也可能导致过度转向。目前旅行车的转向系断开点对转向摆臂球铰的垂直距离(近似为车轮中心与通过摆臂与车身连接的紧固螺栓的中心平面的垂直距离差)比原型车大23mm,这改变了原型车的原有的转向系与悬架运动学关系,在稳态回转试验过程中给转向轮施加一个更大的附加转角,增加车轮前束,加剧过度转向趋势。

4 改善旅行车稳态回转特性的建议

从旅行车和原型车的操纵稳定性影响因素的探讨可知,整车质心位置、前后悬架参数差异以及悬架的运动干涉都可能导致旅行车和原型车在稳态回转特性方面存在差异。针对旅行车稳态回转特性存在的问题,应利用机械系统动力学分析软件建立旅行车的整车操纵稳定性计算机仿真模型,实现对操纵稳定性的动态仿真研究,选取工程上切实可行的参数为设计变量,进行汽车稳态回转特性的计算机仿真分析,确定各参数对汽车稳态回试验的中性转向点的侧向加速度值、不足转向度、车箱侧倾度评价指标的影响,提出工程实用性强的改进设计方案,并进行试验验证,从根本上解决中性转向点的侧向加速度过低的问题,改善旅行车的稳态回转特性。

参考文献

[1]QC/T480-1999.汽车操纵稳定性指标限值与评价方法[S].

[2]某型旅行车车操纵稳定性试验研究总结报告[Z].南京:某汽车公司研究所,2000.

[3]余志生.汽车理论[M].北京:机械工业出版社,2001.

非稳态特性 第8篇

关键词：爆破掘进,粉尘,运移规律,非稳态

0 引言

金属矿山爆破掘进中产生的大量粉尘严重威胁井下工人的健康。由于爆破中粉尘的运动是一个非稳态过程, 受爆破危险性制约, 目前对爆破后粉尘运移分布规律还不够明确, 掌握的除尘技术参数不够确切, 使除尘效果不理想[1,2]。本文以甘肃某金属矿1100爆破掘进工作面为例, 结合现场实测与数值模拟, 研究粉尘运移规律。

1 数学模型与几何模型的建立

1.1 基本控制方程与湍流模型

根据井下实际条件, 将井下空气视为不可压缩粘性流体, 式 (1) 、 (2) 、 (3) 、 (4) 分别为连续方程、动量方程、能量方程和气体状态方程[3]。

式中, ρ为气体密度, kg/m3;xi、xj为笛卡尔坐标系各方向上坐标, m;t为时间, s;ui、uj为气体在笛卡尔坐标系各方向上速度, m/s;p为有效压力, Pa;μ为动力粘性系数, Pas;T为风 (烟) 流温度, K;h为风流的焓值, J;vsj为气流速度, m/s;qr为烟流辐射热, W;gi、gj为i、j方向重力加速度, m/s2;R为通用气体常数, 取8.314, J/ (kgK) 。

湍流流动是一种复杂非线性流动, 基于Reynolds平均法建立RNG k-ε湍流模型, 较好地处理高应变率及轨迹弯曲较大的流动[4,5,6], 式 (5) 、 (6) 分别为k方程与ε方程。

μt0 (μt) 为漩涡修改前 (后) 湍流粘性系数, Pas;μeff为扩散系数;Ω为考虑漩涡的默认估计值。

1.2 粉尘颗粒运动方程

爆破空间内粉尘的体积分数小于10%, 所以基于欧拉-拉格朗日法采用DPM离散相模型用于模拟。粒子所受质量力、升力与热泳力相对于重力数量级较小, 可忽略不计, 颗粒方程为[7]:

式中:CD为曳力系数;up为颗粒运动速度, m/s;u为流体相速度, m/s;ρp为颗粒密度, kg/m3;dp为颗粒直径, m。

通过气流的瞬时速度来跟踪颗粒的随机轨道, 同时考虑颗粒的湍流扩散[8,9,10]。颗粒轨道控制方程为:

式中:τp为颗粒松弛时间, s;u为气流瞬时速度, m/s。

1.3 网格划分与模型检验

实测发现距工作面45m以外的风速基本不变, 掘进面的作业强度集中在工作面附近, 因此选取60m巷道做模拟。巷道宽5.1 m, 高4.6 m, 半圆拱形。风筒直径为420mm, 距工作面21.8 m。距工作面25m处有一硐室, 硐室另配直径490mm的风筒, 据此建立几何模型如图1所示。

采用非结构化四面体单元划分网格, 为改善模拟精度, 在流动激变区做网格调整, 保证计算效率同时提高计算精度[11,12,13,14,15]。为检验模拟结果的可靠性, 将呼吸带高度处 (距地面1.5m) 模拟风速与实测值做对比, 结果如图2所示, 模拟值与实测值基本一致, 因此认为RNG k-ε湍流模型选择合理, 模拟结果可靠。

2 爆破粉尘运移规律CFD模拟

2.1 气相流场边界条件设定

根据实际条件与计算结果, 对气相流场主要边界条件设定如表1所示。

2.2 气相流场风速分布规律分析

图3为工作面y方向风速分布图, 最高风速出现在x轴 (10~20) m, 距离工作面30m以内平均风速较高, 30m后远离涡流, 平均风速较低。风速沿x方向呈先上升后下降趋势, 这是由于爆破安全距离超出了风筒有效射程。掘进空间存在大面积空气滞留区, 距工作面35m外超出射流的有效范围, 平均风速小于0.06 m/s, 因此应采用措施提高供风量。

2.3 爆破粉尘扩散规律分析

如图4所示, 为爆破后第30min巷道垂向、纵向粉尘扩散情况。粉尘浓度在垂高 (0.5~4) m内呈逐渐降低趋势, 这是由于大颗粒粉尘受风流扰动的上浮效应不明显, 只在较低垂高内运动, 而呼吸性粉尘受风流扰动强烈, 均匀分布于空间。粉尘浓度沿工作面纵向60m内先升高后降低, 浓度升高段在 (10~20) m, 该段为风筒形成的涡流区, 风流变化剧烈, 粉尘难以沉降。

2.4 爆破粉尘沉降规律分析

采用点粒子形式, 对爆破后50min不同粒径粉尘沉降状态进行模拟。如图5所示, (90~200) μm粉尘几乎完全沉降, 风筒底部距工作面 (10~20) m处沉降量相对较少。如图6所示, (10~90) μm粒径的粉尘沉降量大于悬浮量, 粒径在15μm以下的长时间悬浮于工作面40m内, 15μm以上的几乎完全沉降。如图7所示, 10μm以下的粉尘受气流扰动较大, 长时间、大范围分布于沿程空间, 悬浮高度达到呼吸带处, 给人体造成危害, 须采取措施对空气进行净化。

将粉尘按“<10μm、10μm~30μm、30μm~50μm、50μm~70μm、70μm~90μm、90μm~200μm (用145μm代表) ”划分, 统计60m内粉尘沉降量, 结果如图8所示。粉尘沉降量随粒径递增而逐渐减少, 在60m范围内, 沉降量呈“减少-增高-减少-增高”趋势。在15m~20m段, 10μm以下粉尘沉降量较多, 10μm~200μm粉尘沉降量较少, 这是由于风筒射流形成的涡流效应造成二次扬尘。

2.5 爆破粉尘浓度随时间变化规律分析

图9为爆破后不同时刻巷道沿程粉尘浓度情况。t=1min时, 掘进空间内大面积区域粉尘浓度较高, 巷道顶板附近浓度低于内部浓度。t=5min时, 沿程各点粉尘浓度急剧下降, 随垂高增加粉尘浓度递减, 巷道底板附近浓度远高于1.5m以上空间的浓度, 10m~20m内空间粉尘浓度较高。t= (10~50) min, 粉尘浓度降低趋势变缓, 40m外的浓度基本不变且明显低于40m以内的浓度。

爆破初期 (1~5) min内, 大颗粒粉尘迅速沉降, 小颗粒粉尘长时间悬浮于空气中, 浓度下降趋势缓慢。由于风筒出口形成涡流卷吸控制区, 使距工作面 (10~20) m粉尘不易沉降。高浓度区域没有沿着巷道扩散运移的趋势, 长时间停留在距工作面40m内。现有的通风排尘无法使粉尘浓度降到标准规定下, 须采用其它措施控制粉尘, 保证正常作业。

3 结论

1) 巷道内风速沿纵向先升高而后降低, 数值模拟结论经实验验证可靠。距工作面 (10~20) m范围内为空气涡流区, 风速较高, 约30m以后, 风速在0.08 m/s左右, 不满足通风标准, 须采取措施加以改进。

2) 粉尘纵向及垂直方向扩散规律:随垂高在 (0.5~4) m递增, 粉尘浓度逐渐降低。距工作面 (0~60) m内, 粉尘浓度呈先升高而后降低趋势, 局部升高区域为风筒形成涡流控制的激变区域。

3) 粉尘沉降规律:距工作面60m内, (90~200) μm粉尘几乎完全沉降;10μm以下粉尘不易沉降, 沿巷道纵向呈“减少-增高-减少-增高”趋势; (10~200) μm粉尘的沉降量沿巷道纵向呈“减少-增高-减少”趋势。