标枪气动力数值模拟

标枪气动力数值模拟（精选9篇）

标枪气动力数值模拟 第1篇

标枪气动力数值模拟

利用FLUENT软件对标枪进行了风洞试验的数值模拟,避免了真实实验对标枪的`破坏.得到了标枪的气动力数据,并提出了投掷该标枪的最佳出手攻角.

作 者：刘晓军 史振东 马晋超 陈昭怡 LIU Xiao-jun SHI Zhen-dong MA Jin-chao CHEN Zhao-yi 作者单位：太原理工大学理学院,太原,030024刊 名：科学技术与工程 ISTIC英文刊名：SCIENCE TECHNOLOGY AND ENGINEERING年，卷(期)：20099(11)分类号：O355关键词：标枪 气动力 数值模拟

标枪气动力数值模拟 第2篇

垂荡板水动力的数值模拟

将计算垂荡板的水动力问题简化为二维竖直平板的绕流问题,计算了在各种来流情况下垂荡板的受力,给出了在计算模型下流场后尾涡片的发展图景.与试验结果对比后发现,尽管计算模型和试验时的三维模型差距很大,但只需要调整计算条件,计算结果与试验结果即可较好地吻合.这表明,在类似的`情况下对于不同形状的平板,存在可以简化到二维的可能.

作 者：纪亨腾 范菊 黄祥鹿  作者单位：上海交通大学,船舶与海洋工程学院,上海,30 刊 名：上海交通大学学报  ISTIC EI PKU英文刊名：JOURNAL OF SHANGHAI JIAOTONG UNIVERSITY 年，卷(期)： 37(8) 分类号：O352 U661.1 关键词：垂荡板   水动力   数值模拟   涡  

尾矿坝动力数值模拟与永久变形 第3篇

关键词：尾矿坝,ABAQUS,动力响应,永久变形

尾矿库是尾矿堆存场所, 是金属非金属矿山选矿厂不可缺少的配套设施, 是维持矿山生养的重要环保和安全设施[1]。由于尾矿库为具有高势能的巨大危险源, 威胁下游居民及设施的安全, 带来很多环境和灾害问题[2]。因此, 运行期间必须保证尾矿坝的正常工作。

目前, 国外对于尾矿坝的研究领域主要为尾矿坝的构筑类型、筑坝方式以及渗流、固结、爆破、地震等因素对尾矿坝的影响等方面[3,4]。本文根据实际工程提供的地形CAD等高线图进行山体建模, 实现复杂地形的模拟, 并对坝体中的各种尾矿材料进行适当分区, 根据坝高和地形图进行布尔运算, 获得最终尾矿坝模型。然后对Abaqus进行二次开发, 使用Fortran编程语言, 实现尾矿坝的地震动力响应和地震永久变形等方面的数值模拟, 为实际工程安全评估提供一定的指导。

1 计算理论

目前, 三维尾矿坝的地震反应分析主要为等效线性法。由于每个迭代过程是线性的, 普通用户计算成本能接受且计算效率高, 在大多数情况下, 尤其是应变较小、地震强度小于 (0.3~0.4) g时均能提供较合理的结果[5]。

1.1 本构模型

在等效线性模型中, 认为尾矿材料是粘弹性体[6], 尾矿的非线性和滞后性特征主要是通过等效剪切模量和等效阻尼比来反映。它的基本原理是把材料应力应变关系所体现出的滞回环等价地用倾角和面积相等的椭圆来表示。虽然近似的椭圆和实际的滞回环间存在着一定的误差, 但对于一般阻尼水平较低的地震反应还是可以使用的, 最后得到动剪切模量及阻尼比。本文采用沈珠江提出的形式[7]。

最大动剪切模量计算式为:

模量、阻尼比随应变水平的变化, 迭代收敛标准为:

1.2 永久变形分析

根据坝体各单元的静应力状态和地震过程中的动力反应, 结合动三轴试验的残余应变试验结果, 可以得到坝体各单元在地震过程中积累的残余应变势[5]。本文利用沈珠江提出的残余应变增量模型确定坝体单元各时段内的残余应变增量。

式中, Δεvp、Δγp和Δε1p分别为残余体积应变增量、残余剪应变增量和残余轴向应变增量;γd为动应变幅值;s1为静应力水平;N和ΔN分别为总振动次数及其时段增量;μ为泊松比;c1、c2、c3、c4、c5为试验参数。

由于相邻单元之间相互制约作用, 求出的应变势并不是各单元的真实永久应变。需根据永久变形沿最大剪应力面发展的原理将各单元残余应变进行转换, 其基张量表达式:

式中, dεpij为直角坐标系下各残余应变增量;δij为克罗内克记号;sij为应力偏量张量。

等效静结点力{F}为:

式中, [D]为弹性矩阵;[B]为应变转换矩阵。

2 工程实例分析

2.1 工程概况

陕西省拟建某尾矿坝, 初期坝为碾压堆石坝, 坝顶标高+544 m, 坝高86 m, 坝顶宽4 m, 坝轴线长282 m, 上游边坡1∶1.9, 下游边坡1∶2。初期坝顶以上采用尾矿堆筑子坝, 子坝高度2 m, 子坝坝顶宽3m, 子坝上游坡1∶3, 下游坡1∶4.5, 形成堆积坝平均堆积坝坡1∶50, 计算工况为堆积坝标高+604 m, 干滩长500 m, 浸润线埋深为8 m。

2.2 材料参数与地震波的选取

材料邓肯模型参数见表1。等效线性模型参数及残余变形参数见表2[8]。

注:F=0.01。

根据施工场地近场构造图及区域地震构造图, 计算采用地震局合成的100年2%超越概率、峰值加速度为0.189g、特征周期为0.69 s的人工地震波。3条地震加速度时程曲线如图1所示。

2.3 计算模型与材料分区

根据工程提供的地形CAD等高线图进行山体建模, 然后将拟建坝体与山体进行布尔运算, 获得最终坝体模型, 使用4面体单元对坝体进行网格划分。坝体模型与材料分区如图2、图3所示。

为更好地模拟动力反应, 本文对整个坝体进行动力分析, 并给出最大断面上部分节点的时程曲线 (节点A位于初期坝顶, 节点B位于h/2坝高, 节点C位于堆积坝顶) 。

3 有限元计算结果分析

编写邓肯E-μ模型的Umat子程序, 实现ABAQUS对坝体进行静力计算, 其目的是获得地震开始前的应力状态, 提取围压。之后进行模态分析, 得到第一振型的频率为1.18 Hz。再编写等效线性法的Umat子程序, 实现ABAQUS对坝体的动力计算。顺河向、坝轴向和竖直向分别输入图1中的地震波, 将竖直向地震波峰值调整为水平向。然后进行迭代计算, 直到满足收敛条件。

3.1 加速度位移反应分析

在地震作用下尾矿坝坝体节点A、B、C的顺河向加速度时程曲线如图4—图6所示;坝体历程最大加速度云如图7所示。对比分析顺河向最大加速度云图可知, 坝体加速度随高度的增大而增大, 在坝顶中部达到最大, 最大反应值为3.7 m/s2;坝顶振动的水平位移比竖向位移大, 坝顶水平和竖直位移值为0.16 m及0.05 m。C点x、y、z方向位移时程曲线如图8所示。

分析以上变形规律, 可以看出坝体竖向、顺河向最大反应加速度均位于坝体顶部。而且随着坝高最大反应加速度逐渐放大, 有较明显的鞭梢效应。通过监测点的位移时程线可以看出, 尾矿坝振动周期大于输入地震波的周期, 尾矿坝的自振频率与地震波的频率相差较大, 坝体响应并不强烈。

3.2 动应力应变分析

地震过程中, 坝体最大动剪应力云如图9所示, 呈层状分布, 坝体的最大动剪应力出现在最大横剖面的底部, 其值达到0.15 MPa。

地震过程中, 坝体最大动剪应变云如图10所示, 参考点动剪应变时程曲线如图11所示。从图11可以看出, 最大动剪应变随着坝高的增大而增大, 最大值出现在强震阶段。

3.3 永久变形分析

在地震作用下, 尾矿坝结构会产生塑性变形且不可恢复, 造成尾矿坝的永久变形或残余变形。为保证其使用价值和外形美观, 需要对坝体永久变形作定量分析, 为抗震预留坝顶超高提供依据。坝体各方向永久变形云图如图12—图14所示。

结果表明:①垂直水流方向的永久变形主要发生在坝顶附近, 有两端向中间挤压的趋势, 变形沿坝轴线近似对称;②顺水流方向永久变形位移主要以坝顶向下游方向运动为主, 尾矿库沉积滩也有向上游变形情况, 但数值较小;③竖向位移发生在坝顶, 坝体中间位移最大, 震陷率0.24%, 基本符合要求。

4 结论

本文对尾矿库进行了动力计算, 输入地震波峰值为0.189g, 特征周期0.69 s。

①坝体的自振周期大于地震波一阶频率1.18Hz。尾矿坝最大加速度及位移出现在地震荷载历时的强震段, 且节点加速度峰值随着其位置高度的增加而增大。②尾矿坝最大加速度及位移出现在地震荷载历时的强震段, 且节点加速度峰值随着其位置高度的增加而增大。③动剪应力最大值在坝底 (0.15 MPa) , 由坡内向坡外逐渐减小。④垂直水流方向的地震永久变形有两端向中间挤压的趋势, 顺水流方向永久变形位移主要由坝顶向下游方向运动, 尾矿坝沉积滩也有向上游变形, 但数值较小;竖向永久变形较大, 发生在坝顶。总体而言, 尾矿坝坝坡在地震过程中能保持稳定。

参考文献

[1]郭金宏.基于水库与尾矿库的区别[J].黑龙江水利科技, 2012, 40 (11) :72-74.

[2]邓敦毅, 邵树强, 潘建平.自然状态下尾矿坝边坡稳定性的数值模拟研究[J].地下空间与工程学报, 2010, 6 (2) :414-428.

[3]中华人民共和国国家标准编写组.GB 50421—2007有色金属矿山排土场设计规范[S].北京:中国计划出版社, 2007.

[4]Scott MD, Lo RC.Optimal tailings management an Highland Valley Copper[J].CIM Bulletin, 1992, 85 (962) :85-88.

[5]费康张建伟.ABAQUS在岩土工程中的应用[M].北京:中国水利水电出版社, 2010.

[6]沈珠江.一个计算砂土液化变形的等价粘弹性模式[C]//中国土木工程学会第四届土力学及基础工程学术会议论文选集.北京:中国建筑工业出版社, 1986, 199-207.

[7]沈珠江, 徐刚.堆石料的动力变形特性[J].水利水运科学研究, 1996 (2) :143-150.

高速列车气动噪音数值模拟 第4篇

摘要：采用k-ε两方程湍流模型模拟了高速列车在有隔音屏和无隔音屏状态下的粘性流动，将压力变化通过傅立叶变换转化成噪音值，通过相同点的噪音值（dB）比较，表明隔音屏能有效的降低高速列车产生的气动噪音。

关键词：高速列车；气动噪声；fluent；隔音屏

1.引言

随着国家经济的高速发展，传统的交通运输业体现出了越来越大的局限性，为了缓解交通运输的瓶颈，国家开始大力建设高速铁路，越来越多的高速列车驰骋在祖国大地。高速列车一般指时速在200公里以上的火车，目前动车组的平均时速已经超过200公里，而刚刚兴起的大部分高速列车（俗称：高铁）的时速已经可达到350公里每小时。

虽然高速列车能极大的提高运输效率，但同时也带来了一些环境危害，比如说电磁辐射，振动，噪音等等。据2002～2003年对某些小区的调查显示，普通列车运行时产生的噪音达80dB左右[1]。列车产生的噪音源有很多，一部分是机械噪声，比如说列车与铁轨的碰撞，还有一部分是气动噪声，这主要是列车高速运行的时候造成附近压力场的变化而产生的噪声。随着列车时速的提高，气动噪音会大幅增加。据模拟显示，这一部分噪声主要是低频噪声，而低频噪音会对人体健康产生长远的影响[2]，因此需要采取有效措施降低高速列车产生的低频噪音。

计算流体力学（Computational fluid dynamic）是最近几十年兴起的一门新兴学科，该方法是通过建立计算域网格，离散控制方程，利用计算机强大的运算能力迭代求解控制方程，进而得到流场的各种物理量。目前成熟的计算流体力学商业软件已经非常成熟，业界采用比较多的主要有fluent、star-cd、cfx等软件，其中fluent以其强大的适应性和良好的计算精度得到了广泛的应用，对于噪音的模拟，已经有相关文献证明了fluent在这方面的适用性[3][4]。

本文主要利用fluent对高速列车周围的流场进行了模拟，并通过压力场的变化计算出噪音分布，通过比较相同点在有隔音屏和无隔音屏时噪音大小的比较，证明了fluent在计算高速列车气动噪音的适用性和隔音屏降低噪音的有效性。

2.控制方程与边界条件

2.1控制方程

在直角坐标系中，三维守恒型NS方程可以写为：

Ut+Ex+Fy+Gz=1Re（Evx+Fvy+Gvz）（1）

式中U=（ρ，ρu，ρv，ρw，e）T，E，F，G为对流项，Ev，Fv，Gv为粘性项，其中ρ，u，v，w，e分别表示气体的密度、x方向的速度分量、y方向的速度分量、z方向的速度分量和单位体积的总内能。对于理想气体p=（γ-1）[e-ρ2（u2+v2+w2）]，对于空气，比热比γ取1.4。

对于粘性系数，本文采用工业上常用的k-ε湍流模型进行计算，其控制方程如下：

紊动能k方程：

ρDkDt=xi（μ+μtσk）kxi+Gk+Gb-ρε（2）

紊动能耗散率ε方程：

ρDεDt=xi（μ+μtσε）εxi+C1εεk（Gk+C3εGb）-C2ερε2k（3）

2.2边界条件

物面边界条件：

对于粘性流动，物面上应该满足无滑移条件，即u=v=w=0 。此外，假设物面为绝热壁条件，即Tn=0，这里n为物面外法线方向。

远场边界条件：

在实际数值计算中，计算区域的大小是有限的，因而在计算区域的远场边界处需要引入无反射边界条件，以保证物体产生的扰动波不被反射回内场这一流动的物理特征，在fluent中的压力远场边界条件采用Riemann不变量关系来处理远场无反射边界条件，其远场边界处的Riemann不变量定义如下：

R-=qn-2aγ-1；R+=qn+2aγ-1（4）

其中qn表示边界上的外法向速度分量，a表示当地的音速，γ 为气体比热比。

3.算例分析

本文采用商业软件fluent模拟了高速列车在无隔音屏和有隔音屏情况下的流场及其产生的气动噪音。

算例一：高速列车在无隔音屏情况下流场和气动噪音的模拟模拟。计算状态为来流速度u=100m/s。本文采用fluent的前处理器gambit生成计算网格，图1为表面网格分布图。考虑到速度100m/s已经可以认为是可压缩流动，因此求解器选择耦合求解器，空气密度设置为理想气体。图2为压力等值线分布图，从图中可以看出列车头部压力变化非常剧烈。为了比较隔音屏对噪音传播的影响，设置了两个监测点，点A位于隔音屏外，距离列车中心10米，高度4米，点B位于隔音屏内，距离列车中心5米，高度4米。根据相对性，计算时将列车设为固定，那么监测点相对于列车就是运动的，运动轨迹为一条直线。图3为该直线上的压力分布。通过将该直线距离与压力的关系除以列车的速度，就得到了监测点时间与压力的关系，将其通过fluent的傅里叶变换，就得到了监测点的噪音值（dB）。

算例二：高速列车在有隔音屏情况下流场和气动噪音的模拟。计算状态与无隔音屏一致。图4为表面网格分布图（为了显示方便，将一侧的隔音屏隐藏）。图5为压力等值线分布图。

图6为监测点压力分布图。通过与图3比较，可以看出有隔音屏的情況下，监测点A和B的压力差更加明显，证明隔音屏能有效的降低压力的传播。图7为监测点噪音值分布图，从图中可以看出，对于频率20HZ左右的低频噪音，加了隔音屏之后，监测点A处的噪音值从80dB降低到了60dB左右，而监测点B处的噪音值始终在90dB左右。由此可见，隔音屏对于降低低频噪音还是有相当效率的。

4.结论

本文利用fluent软件模拟了高速列车在有隔音屏和无隔音屏状态下的气动噪音，通过相同点噪音值（dB）的比较，表明隔音屏能有效的降低列车通过时产生的气动噪音，同时也证明了fluent对于气动噪音计算的适应性。出于简化的目的，本文的计算模型采用的是简化后的列车车头，因此未能捕捉到小的压力脉动，而这一部分压力脉动是产生高频气动噪音的噪音源，因而采用实际列车模型，利用非定常计算捕捉小的压力脉动，计算高频段的气动噪音将是下一步工作的重点。（作者单位：四川航空股份有限公司机务工程部）

参考文献：

[1]魏保祥，耿小娟，唐凤琴. 火车噪声、振动对8个小区环境影响调查，中国环境卫生，2004年第7卷第2期：8～10

[2]孟苏北. 城市住宅区低频噪声对人类健康的危害，中国医学导报，2007年第4卷第35期：17～19

[3]齐海政，刘献栋. 汽车排气消声器仿真方法探究，噪声与振动控制，2006，6：73～76

标枪气动力数值模拟 第5篇

混沌摆系统的动力学分析和数值模拟

在混沌实验摆动力学方程基础上,分析其分别处于保守和耗散非线性系统的动力学特性. 采用实际物理实验系统的参数进行数值模拟,探讨不同的策动频率、振幅和外加阻尼等参数条件对混沌摆系统动力学行为的影响,且分析这些参数在实验过程中所起的`调节作用. 数值模拟结果表明:混沌实验摆系统动力学性质灵敏地依赖于策动频率值;系统策动振幅必须大于一定阈值是混沌相出现的必要条件;磁阻尼参数对系统混沌相动力学状态影响不明显.

作 者：朱桂萍 王健 ZHU Gui-ping WANG Jian 作者单位：扬州大学,物理科学与技术学院,江苏,扬州,225002刊 名：扬州大学学报(自然科学版) ISTIC PKU英文刊名：JOURNAL OF YANGZHOU UNIVERSITY(NATURAL SCIENCE EDITION)年，卷(期)：11(3)分类号：O415.5关键词：混沌摆 动力学特性 实验参数 数值模拟

标枪气动力数值模拟 第6篇

定常情况下二维叶栅气动性能的数值模拟

以某型航空发动机叶轮机为研究对象,采用了周期性边界条件和Spalart-Allmaras湍流计算模型,通过数值模拟,对定常情况下二维叶栅气动性能进行了数值模拟,结果表明:气流转折角的变化是与攻角相关联的,并且进口静压、叶栅负荷和损失都随着攻角的变化而变化.

作 者：阎成鸿  作者单位：民航飞行技术与飞行安全科研基地,四川广汉,618307 刊 名：中国民航飞行学院学报 英文刊名：JOURNAL OF CIVIL AVIATION FLIGHT UNIVERSITY OF CHINA 年，卷(期)：2008 19(2) 分类号：V231 关键词：航空发动机   叶轮机   二维叶栅   气动性能  

标枪气动力数值模拟 第7篇

根据流固耦合的Taylor平板理论和Cole的水下爆炸经验公式,得到了水下爆炸冲击波作用下平板迎爆面的总压力载荷,并将该压力载荷嵌入到ABAQUS程序的用户自定义载荷子程序VDLOAD中,实现了对ABAQUS的二次开发.作为验证,利用该方法模拟了水下爆炸冲击波作用下固支平板的塑性动力响应,模拟结果与试验结果吻合良好.与通用的.模拟方法相比,本文所采用的模拟方法在保证较高计算精度的前提下,可以显著降低计算成本并简化建模过程,对于工程应用具有一定参考和实用价值.

作 者：贾宪振 胡毅亭 董明荣 许学忠 刘家骢 JIA Xian-zhen HU Yi-ting DONG Ming-rong XU Xue-zhong LIU Jia-cong 作者单位：贾宪振,胡毅亭,刘家骢,JIA Xian-zhen,HU Yi-ting,LIU Jia-cong(南京理工大学,化工学院,江苏,南京,210094)

董明荣,许学忠,DONG Ming-rong,XU Xue-zhong(西北核技术研究所,陕西,西安,710024)

标枪气动力数值模拟 第8篇

离心风机将机械能转化成气体的压力能和动能的流体机械, 在农产品加工、石油、化工、污水处理等众多行业中有着广泛的应用, 并被形象地誉为“国民经济的呼吸系统”。对于单级离心风机, 国内外很多研究者通过数值模拟或者实验的手段, 期望通过研究新的结构或者优化来提高风机性能, 研究对象涉及到叶轮、蜗壳、集流器等。于跃平[1]等分析了不同叶片型线叶轮对风机的最高压力点、最佳工况点、压力曲线形状 及效率的 影响。李建峰[2]等通过对9 - 26型风机加装短叶片, 在相同的进出口条件下, 增加短叶片的风机流量可增加约5% , 出口全压平均增加约10% 。多级离心风机用于压力需求比较大的场合, 但是目前已发表的关于多级离心风机的文章, 数量特别少, 如王企鲲[3]等对多级离心压缩机级间静止部件气动性能的CFD分析。

文中主要从压力损失的角度出发, 计算发现叶轮是发生压力损失的主要零部件, 主要研究了叶轮叶片的进口安装角βb1、出口安装角βb2以及叶片z数量对效率η的影响, 为了考虑3个参数的交互影响, 本文基于响应面方法, 通过二次回归得到了关于效率η的函数式, 找到了最高效率点相对应的3个参数。

1 多级离心的结构

八级的离心风机结构如图1所示。图中省略了轴承、轴等, 气体首先进入前盘, 然后经过第一级的叶轮做功, 气体获得动能和压力能, 经过弯道和导叶, 流出第一级, 进入第二级继续获得更大的压力, 最后由后盘流出整个风机。为得到较大的压力, 此风机叶轮采用的是前向叶轮, 叶片采用单圆弧结构, 叶片数量为16。叶轮的详细参数如表1所示。

2 数值计算方法

2. 1 计算区域及网格

流体流过的典型过程依次是经过前盘、叶轮、弯道、第一级导叶, 然后进入第二级重复此过程, 因每一级的零件参数都相同, 所以某一级的气动性能完全能反应整机的性能, 因此数值计算的模型采用整个风机的第一级, 包括叶轮、叶轮和前盘形成的间隙以及导叶与前盘形成的无叶扩压器、前盘和中盘形成的弯道、导叶等4个部分, 计算模型如图2所示。

对于每个部分的流体域, 既可以由布尔操作形成, 如叶轮的有限元模型, 也可以直接建立流体域的三维模型, 如间隙以及扩压器部分。在用gambit网格划分过程中, 可以对每个部分再进行分区画网格, 以得到较好的网格质量, 整个模型采用六面体网格, 其中叶轮和导叶采用Cooper ( 源面) 方式, 其余的均采用Map或者Submap方式划分。为得到合适的网格大小, 进行了多次网格划分, 通过不断调整网格大小, 直到两次结果相差很小时, 得到满意的网格尺寸。整个计算模型网格数量654645, 所有的网格质量都控制在0. 67以内。布尔操作形成的叶轮流域有限元模型如图3所示。

2. 2 求解器、离散格式及边界条件

选取Standard k - epsilon模型, 近壁面采用标准壁面函数, 压力与速度耦合采用Simple算法, 动量方程、耗散率输送方程以及湍动能方程采用一阶迎风格式。在入口 处给定流 量入口条 件qm=1. 74kg / s, 出口为自由出流, 叶轮流域和动壁面转速为2900r/min。

2. 3 气动性能评价指标

内效率是评价气动性能的重要标志, 可由式 ( 1) 确定:

式中: η—内效率;

p2、p1—分别为出口和进口的全压;

qv—体积流量;

T—叶轮对流体的作用力对旋转轴的力矩;

w—转速, 计算式各个参数均采用国际单位。

经Fluent数值模拟, 得到原风机的结果为:p2= 6355Pa, p1= - 1790Pa, T = 101. 89N·m , 计算得到内效率η = 39. 49% 。

2. 4 各零部件的损失压力计算

在一个流道中, 总机械能头大的一段必然要流向小的一端。若忽略重力的影响, 伯努利方程可以表达流体在流动过程中遵守的能量守恒关系, 即:

式中: p1—流场中某一点的压强;

v1—某一点的速度;

p2、v2—另外一处对应的物理量;

w—外界对叶轮所做的功;

Δw—损失压力。

对于各部分的损失压力, 都可由Fluent报告的数据和已知数据经计算得到, 其中ρ = 1. 16kg /m3。

1) 对叶轮: p1= - 1962Pa, p2= 14957Pa, v1=36. 18m / s, v2= 34. 96m / s, w = 20627Pa。

由式 ( 2) 得到: Δw1= 3788. 34Pa。

2) 对弯道: p1= 8729Pa, p2= 8643Pa, v1= v2=16. 22m / s, w = 0。

由式 ( 2) 得到: Δw2= 83Pa。

3) 对导叶: p1= 7493Pa, p2= 6355Pa, v1=17. 48m / s, v2= 24. 88m / s, w = 0。

由式 ( 2) 得到: Δw1= 956. 19Pa。

由此可知, 叶轮是压力损失最多的零部件, 因此重点研究了叶轮参数对效率的影响, 并优化了叶轮的参数。

3 变量对气动性能的影响分析

3. 1 入口安放角 βb1对效率 η 和全压 P 的影响

在出口安放角βb2= 150°, z = 16不变的情况下, 研究了βb1和η以及βb1和P之间的关系。取5个不同的βb1, 数值模拟得到的结果如图4所示。

由图4可知, βb1对效率的影响不大, 最大值为39. 62% , 最小值为38. 91% , 并且从45° ~ 65°的变化过程中, 存在一个最佳效率点。而对于全压, 呈递增趋势。

3. 2 叶片数量 z 对效率 η 和全压 P 的影响

在叶片入口安放角βb1= 55°, 出口安放角βb2=150°不变的情况下, 效率和全压随叶片数量变化规律如图5所示。

由图5可知, 在入口安放角, 出口安放角一定的情况下, 效率随着叶片数量的增加而减少, 全压随着叶片数量先增大后减少, 对此可以解释为: 随着叶片数量的增加, 漩涡损失减小, 而摩擦损失增大, 当两者损失之和最小时即为最佳叶片数, 因此全压呈现了一个峰值, 而效率除了与全压有关之外, 还与转矩有关, 叶片数的递增, 使转矩不断增加, 两者综合体现出效率呈下降趋势。

3. 3 出口安装角 βb2对效率 η 和全压 P 的影响

在入口安放角βb1= 55°, 叶片数量z = 15不变的情况下, 效率η和全压P随出口安放角的变化规律如图6所示。

由图6可知, 全压随着出口安放角的增大而增大, 而叶轮对流体作用力对旋转轴的转矩T也随着增大, 由式 ( 1) , 两者共同体现出效率呈下降趋势。

4 响应面方法以及样本点的设计

4. 1 响应面方法

响应面法 ( RSM) 是数学方法和统计方法结合的产物, 是用来对所感兴趣的响应受多个变量影响的问题进行建模和分析的, 其最终目的是优化该响应值, 它是一种近似方法, 而且在对变量在大区间变化时, 仍然可以发挥其作用, 因而在食品、工程、生物学、医学等很多领域都得到了广泛的应用。对于本文3个变量的响应面法, 可由式 ( 3) 表达响应值y和变量之间的关系。

式中: x—样本的变量; β—系数, 可通过样本点数据求出。

由于是一种近似, 且是由有限的样本点去表达变量在指定区间的表达方式的方法, 因而会存在拟合误差, 其误差大小用统计学中的决定系数R2来评价。拟合误差越小, R2越接近1。

4. 2 变量编码及样本点的设计

将βb1、βb2、z这3个变量编码至 ( - 1, 1) 产生设计空间, 如表2所示。样本点的选取为Box Behnken设计方法, 其试验点的选取为编码立方体的棱的中点, 它是一种常用的回归设计法, 适用于2 ~ 5个变量的优化实验。样本点以及对应的数值模拟结果如表3所示。

利用非线性拟合的方式得到效率η的二次回归方程为:

且计算R2= 0. 9962, 很接近1, 拟合效果很好。由回归方程 ( 3) 求得的最大值处, x1= 0. 3853, x2=- 0. 5197, x3= - 1, 即βb1= 58. 85°, βb2= 138. 6°, z= 12, 此时效率最大值为0. 4555。

4. 3 优化后结果与数值模拟的对比

以优化以后的参数模型再进行数值模拟分析, 经过后处理报告得到全压和转矩, 由式 ( 1) 计算效率为0. 4539, 原风机效率、响应面优化结果以及数值模拟结果的对比如表4所示。

%

5 结论

通过对3个变量βb1, βb2, z进行优化, 通过Fluent数值模拟, 可得到以下结论:

1) 在没有得到变量和目标 ( 本文为效率) 确切表达式的情况下, 利用响应面方法得到了变量和目标之间的关系, 并求得了效率的最大值。

2) 由Fluent数值结果以及响应面优化结果, 两者误差为0. 16% , 可以说明数值模拟结果的可行性, 为风机的设计提供一定的参考依据。理论上使得优化以后的风机效率比原风机提高5. 9% 。

摘要：用Fluent对多级离心风机的某一级进行数值模拟计算, 由伯努利方程发现, 叶轮是压力损失最大的零部件, 因此重点研究了叶轮叶片进口安放角、出口安放角、叶片数量3个变量对效率的影响。进一步应用响应面模型优化方法, 用二次回归得到了效率关于3个变量的函数关系式, 计算出的函数最大值与数值模拟方法误差为0.16%, 且比原风机的效率提高了5.9%, 为风机的改进设计提供了参考依据。

关键词：多级离心风机,响应面模型,二次回归

参考文献

[1]于跃平, 胡继孙, 陈启明, 等.叶片型线对离心风机气动性能影响试验研究与叶轮流场计算[J].流体机械, 2007, 35 (7) :25-29.

[2]李建峰, 吕俊复, 张海, 等.增加短叶片的9-26型风机流场数值模拟[J].流体机械, 2007, 35 (10) :29-32.

[3]王企鲲, 戴韧, 陈康民.多级离心压缩机级间静止部件气动性能的数值研究[J].动力工程, 2005, 25:102-104.

[4]丁问司, 何祥滨.吸粮机三级离心风机内部流场的数值模拟[J].农业工程学报, 2011, 27 (11) :78-84.

[5]石硕, 张晓菲, 张国成, 等.多级离心风机的气动优化设计[J].山东建筑大学学报, 2012, 27 (2) :242-245.

[6]王松岭, 雷泳, 李春曦, 等.基于Fluent的离心风机内部流场三维数值模拟及节能改造研究[J].动力工程, 2005, 25:87-90.

[7]王永菲, 王国成.响应面法的理论与应用[J].中央民族大学学报, 2005, 14 (3) :236-239.

标枪气动力数值模拟 第9篇

摘要：基于显式动力非线性有限元分析方法，利用ANSYS/LSDYNA软件对强夯问题进行分析，得出了显式动力非线性数值模拟方法的一般步骤.结合某路基强夯实例建立三维立体模型，对碰撞过程进行数值模拟，得到了强夯加固范围及夯后土体的应力场、位移场.通过与现场实测数据的比对，验证了显式动力非线性有限元数值模拟方法在强夯问题中的适用性.在此基础上，研究和探讨了岳阳地区强夯处治粉煤灰路基中的夯锤参数选取问题，分析比较了夯后土体的沉降，结果表明同能级下重锤低落距有更好的加固效果.

关键词：强夯；数值模拟；大变形；粉煤灰路基；参数分析

中图分类号：TU470.3 文献标识码：A

随着工业化、城镇化的推进，土地资源愈发紧张，大量工程不得不建立在粉煤灰堆积地区.而粉煤灰地基的承载力较低，往往不能满足工程需要，因此必须进行处理，以使其满足承载力、稳定性和抗变形能力的要求.强夯法凭借其施工工艺简单、经济性好、加固效果明显等优点＼[1＼]在粉煤灰地基处治中应用＼[2＼]较多.

然而，与广泛的工程应用相比，强夯的理论研究稍显滞后.强夯过程的复杂性决定了其难以用精确的解析方法求解，其加固机理和设计理论方面尚不成熟，现场施工主要依靠经验公式和试验确定＼[3＼]，大大限制了强夯法的发展和应用.针对这一状况，国内外学者采用数值方法，对强夯过程及夯后土体应力场、位移场等进行了较为系统的研究.

钱家欢等＼[4＼]运用加权余量法导出了弹性振动问题的边界积分方程，并将其应用于边界元分析强夯问题.Chow等＼[5＼]基于一维波动方程模拟夯锤和土体之间的相互作用，得出应力波的传播特性.孔令伟等＼[6＼]在考虑夯锤自重的基础上，结合夯锤刚体运动方程和成层弹性地基空间轴对称动力问题的传递矩阵法，导出了强夯的边界接触应力与沉降在变换域中的解析式.然而，上述方法均是基于小变形假定，在工程实际中，夯击区将产生较大破坏，上述方法得出的结果往往不尽如人意.

随着数值计算方法的发展，显式动力非线性有限元理论＼[7＼]日趋成熟，为强夯这类非线性大变形问题的求解提供了一种新的选择.Thilakasiri等＼[8＼]率先采用显式动力非线性有限元软件LSDYNA2D对强夯置换软土进行了数值模拟，得到了夯后土体的应力应变关系，但其采用的仍是二维平面模型.国内研究者于德水＼[9＼]、杨建华＼[10＼]、张建辉和杨培轩＼[11＼]分别利用ANSYS/LSDYNA分析了强夯加固湿陷性黄土、碎石土、风成砂的过程，对强夯法的加固机理及其加固范围进行了探索.

可见在常见土类的强夯研究中，显式动力非线性数值方法已经占有一席之地，取得了一定的成果.相较而言，其在粉煤灰地基中的应用目前尚属空白.

鉴于此，本文首先对显式动力非线性有限元数值模拟的一般方法进行探讨，并通过实例验证其在强夯中的适用性，得到夯后土体应力、应变场及夯沉量.在此基础上，结合岳阳某粉煤灰路基强夯工程，对工程中的夯锤参数方案选择问题进行研究.

1显式动力非线性有限元分析方法

1.1显式动力非线性有限元方法简介

在动力学中，显式算法与隐式算法是两个相对的概念.显式算法主要包括central difference method （即中心差分法），是ANSYS/LSDYNA中的主要求解方法，用于分析大变形、瞬态问题、非线性动力问题.

结构系统的通用运动学方程为：

M+C+KU=Rt.（1）

式中：U，Rt，M，C，K分别为结构位移、荷载、质量、阻尼、刚度矩阵.

假定0，t1，t2，…，tn 时刻的节点位移、速度及加速度已知，现求解t+Δt时刻的结构响应.中心差分法对速度、加速度采用中心差分代替，即为：

t=1Δt2（Ut-Δt-2Ut+Ut+Δt），（2）

t=12Δt（Ut+Δt-Ut-Δt）.（3）

将式（2）和式（3）代入式（1）中，整理可得：

Ut+Δt=t.（4）

式中：

=1Δt2M+12ΔtC；（5）

t=Rt-（K-2Δt2M）Ut-（1Δt2M-

12ΔtC）Ut-Δt.（6）

式（5）和式（6）分别被称为有效质量矩阵、有效载荷矩阵.中心差分法在求解t+Δt瞬时的位移Ut+Δt时，只需t+Δt时刻以前的状态变量Ut和Ut-Δt，然后计算出有效质量矩阵、有效荷载矩阵，即可求出Ut+Δt，故称此法为显式算法＼[12＼].

显式算法的优点是它既没有收敛性问题，也无需求解联立方程组，其缺点是时间步长受到数值积分稳定性的限制，不能超过系统的临界时间步长.由于强夯是瞬态非线性过程，从解的精度考虑，时间步长也不能太大，这就在很大程度上弥补了显式算法的缺陷.

1.2强夯模型的建立

1.2.1数值模拟简化

强夯过程复杂，影响因素多.为简化计算，在建立强夯碰撞模型时作如下假设：

1） 路基填土土体均质、各向同性；

2） 与地基土相比，夯锤在冲击过程中被认为是刚体，变形忽略不计；

3） 不考虑夯击过程中产生的热能和声能能量损失；

4） 锤体冲击过程是瞬态大变形问题，孔隙水压力的影响很小，故不考虑孔隙水压力.

受篇幅限制，本文只模拟强夯过程中单点夯击的第一击.

1.2.2材料本构模型

夯锤采用刚体材料即LSDYNA3D中的020RIGID.土体的材料模型是本文的一个难点，目前还没有一个可以很好地反映高能量冲击下土体应力应变特征的本构关系.传统的MohrCoulomb模型应用于大变形问题时，计算结果往往不能收敛.彭建兵等＼[13＼]用DruckerPrager本构模型（简称DP模型）研究冲击荷载下黄土的动力响应问题，取得了较为理想的结果.本文也采用DP模型作为土体本构模型，并通过与实测数据的比对来验证其在强夯数值模拟中的适用性.

1.2.3建模及网格划分

模型包括锤体和土体两部分，本模型可按轴对称进行简化，夯锤和土体均取1/4建模，采用3D solid164八节点六面体实体单元.

由于强夯是瞬态问题，因此夯击影响到的土体深度和宽度是有限的.通过加入适当的边界条件即可模拟无限大土体的情况.对于土体模型深度和宽度的考虑可以根据经验公式和现场试验确定.

为了既能保证计算结果的准确性又尽量缩小模型规模，离夯击中心点较近处，网格较密，离中心点较远处，网格较疏.

强夯碰撞计算时间取0.2 s.时间步长设定为0.4 s，通过在LSDYNA中定义*TIMESTEP下的参数TSSFAC完成.

为消除沙漏模式的变形积累，在LSDYNA中添加沙漏控制卡*HOURGLASS，选择4号模式.

1.2.4接触面、边界条件的建立及荷载的施加

接触时将发生穿透，故采用侵蚀接触选项，通过在LSDYNA3D中定义关键字*CONTACT_ERODING_SURFACE_TO_SURFACE完成.设置夯锤接触面为Master segment，土体接触面为Slave segment.

土体约束分为两部分.在两个对称面（即xz面和yz面）分别约束y方向和x方向位移，底面约束z方向位移，土体模型外侧两个面设置为无反射边界条件*NON_REFLECTING.

本文不模拟夯锤的自由落体，只需对夯锤施加一个初速度即可.重力加速度通过在LSDYNA中定义*Load中的关键字*GRAVITY_PART完成.

1.3工程实例验证

为验证显式动力非线性有限元数值模拟方法的可行性，对某高速公路填土路堤强夯实例建模分析如下.夯锤为铁锤，土体物理力学指标通过试验获得，具体数据见表1.

夯击能为1 200 kN·m，夯锤落距H=10 m，夯锤初始接触速度v=14 m/s.土体水平方向建模6 m，竖向建模10 m，划分539 000个实体单元.夯锤半径1 m，高0.5 m，划分为81个实体单元，模型建立如图1所示.

1.3.1土体的应力、应变云图

图2所示为土体竖向应力云图.如图所示，当夯击发生时，动应力波逐渐从锤底向四周扩散，应力波在土体中传播，使得土体原有结构发生改变，孔隙体积减小.竖向应力等值线云图上小下大，成一梨形分布.这与文献＼[14＼]的实测动应力分布规律相一致.

图3所示为土体竖向位移云图.夯锤碰撞土体使土体发生变形的过程中，土体的位移主要发生在锤底的位置，并向两侧发展.竖向变形随着深度的增加而减小，深度越大，衰减速度越快.土体达到最大位移后会发生一定回弹.同时，土体的竖向位移随水平位移的增大而减小，且衰减的速度很快.

1.3.2不同深度处动应力随时间的变化

图4所示为不同深度处动应力时程曲线.如图所示，在夯击过程中土体的应力波为一尖峰，没有明显的第二应力波.当z=0时，峰值最大，达到2.251 MPa.峰值持续时间很短，且随深度的增加衰减很快.当z=8 m时，应力波峰值为0.158 MPa，与z=0时相比衰减明显.同时，土体动应力沿深度依次达到峰值，反映了应力波在土体中的传播性质.

1.3.3土体竖向位移随深度的变化

图5和图6所示为中心点竖向位移随深度的变化情况.如图所示，随着深度的增加，土体竖向位移减小.当达到最大位移后，土体会发生一定的回弹，稳定后达到最终夯沉量.位于夯锤下方的土体位移最大，回弹也较大，最终夯沉量为0.150 m.当深度达到6 m时，竖向位移仅为3.7 mm.可见，随着深度的增加，土体的竖向位移迅速衰减.由图5可知，土体达到最大位移的时间随深度的增加而推后，这是由于在碰撞中产生了应力波，应力波从夯锤底部向下传播，使得夯锤下土体的位移依次达到最大值.

1.3.4土体竖向位移随水平距离的变化

图7所示为深度0 m， 3 m， 6 m处土体竖向位移随水平距离的变化情况.如图所示，同一深度处，离中心点的水平距离越远，竖向位移越小.同一水平距离处，深度越大，竖向位移越小.当水平距离超过3 m时，可以发现图中3条曲线纵坐标均趋近于0，说明该能级强夯水平影响范围在3 m左右.同时可以观察到，夯击后地表离中心点2 m左右处土体会发生一定的隆起，这与实际情况是一致的.

图8所示为地表夯沉曲线.如图所示，工程实测地表夯沉量为0.131 m＼[15＼]，数值模拟结果为0.150 m，误差为1.9 cm，相对误差14.5%.从整体曲线趋势来看，均反映了随着水平距离的增加沉降衰减以及2 m左右土体发生隆起的现象.

图8所示为地表夯沉曲线.如图所示，工程实测地表夯沉量为0.131 m＼[15＼]，数值模拟结果为0.150 m，误差为1.9 cm，相对误差14.5%.从整体曲线趋势来看，均反映了随着水平距离的增加沉降衰减以及2 m左右土体发生隆起的现象.

2.2强夯处治粉煤灰路基模型的建立

方案1采用轻锤高落距方案.夯锤半径1 m，高0.5 m，落距20 m.夯锤划分为81个实体单元，土体水平方向建模6 m，竖向建模10 m，划分539 000个实体单元.夯锤初始接触速度为v=19.796 m/s.夯击中心处网格加密，边界条件等设定同前，模型如图9所示.

方案2采用重锤低落距方案.夯锤半径1 m，高1 m，落距10 m.夯锤划分为162个实体单元，土体水平方向建模6 m，竖向建模10 m，划分539 000个实体单元.夯锤初始接触速度为v=14 m/s.夯击中心处网格加密，边界条件等设定同前，模型如图10所示.

2.3计算结果比对及分析

就同种土而言，比较强夯加固效果最直观的方法就是比较夯沉量.夯沉量越大，压实度越高，土体工程性质的改善就越为明显.本文将以此为切入点来比较两种方案的加固效果.

2.3.1土体竖向位移随深度的变化

图11和图12所示为中心点不同深度处竖向位移时程曲线.如图所示，z=0时，方案1土体最大竖向位移为0.423 m，到达时间为0.085 s；方案2土体最大竖向位移为0.581 m，到达时间为0.123 s.z=2时，方案1土体最大竖向位移为0.150 m，到达时间为0.091 s；方案2土体最大竖向位移为0.231 m，到达时间为0.138 s.其余深度情况类似，在此不一一赘述.在相同深度处，方案2的最大竖向位移均大于方案1.同时，方案2的曲线更加平缓，到达最大位移的时间较长，这意味着锤体进人土体后相对减速慢，土体所受到的振动持续时间更久，加固时间更长.

图13所示为竖向位移随深度变化情况.如图所示，在相同深度处，重锤低落距方案土体最终沉降要大于轻锤高落距方案土体最终沉降.同时，两条曲线均显示了土体竖向位移随深度增加而衰减的规律.当深度达到8 m时，两种方案的沉降量均已很小，曲线衰减已经很慢.

图13中心点竖向位移随深度的变化曲线

Fig.13The middle settlementdepth curves

2.3.2地表沉降随水平距离的变化

图14所示为地表夯沉曲线.如图所示，轻锤高落距方案地表夯沉量为0.382 m，而重锤低落距方案地表夯沉量为0.552 m，重锤夯沉量比轻锤夯沉量大40%以上.相同能级下，重锤夯击的冲量更大，在锤体进人土体后，相对减速慢，相应的单击夯沉量大，土体压密加固效果好于轻锤高落距，可以产生比轻锤夯击更大的夯坑.同时，轻锤夯击时，周围土体的隆起现象较重锤夯击明显，这在工程中是较为不利的.

3结论

本文对显式动力非线性有限元数值方法的一般步骤进行了探讨，通过实例验证了其在强夯问题中的适用性.在此基础上，探讨了强夯加固粉煤灰中的夯锤参数选取问题，并得出以下结论：

1）数值模拟结果表明，土体竖向变形随着深度和水平位移的增加而减小.锤径2 m且能级为1 200 kN·m时，强夯加固高填土路基的有效影响深度约为6 m，夯点间距不宜超过3 m.

2） 数值模拟结果与工程实测结果基本一致，说明了DP模型可以较好地描述土体在强夯冲击过程中的特性，同时也说明了显式动力非线性有限元方法在强夯问题中的适用性，可以为强夯设计和理论研究提供一定的依据.

3） 通过对强夯处治粉煤灰方案的建模研究发现，在夯击能均为2 400 kN·m且锤径相同的情况下，重锤夯击的单击夯沉量大，加固效果要好于轻锤夯击.同能级下，增大夯锤重量，同时相应减小落距，可以有效减少单点的总夯击次数，提高强夯机具的工作效率，实现缩短工期、降低工程成本的目的.

参考文献

[1]赵明华. 土力学与基础工程＼[M＼]. 2版. 武汉： 武汉理工大学出版社， 2003： 329-331.

ZHAO Minghua. Soil mechanics and foundation engineering＼[M＼]. 2nd ed. Wuhan： Wuhan University of Technology Press， 2003： 329-331.（In Chinese）

[2]王起刚. 强夯法在粉煤灰地基中的应用研究＼[D＼]. 青岛： 中国海洋大学， 2003： 17-26.

WANG Qigang. Study on dynamic consolidation method in coal subground practice＼[D＼]. Qingdao： Ocean University of China， 2003：17-26. （In Chinese）

[3]龚晓南. 地基处理手册 ＼[M＼]. 2版. 北京： 中国建筑工业出版社， 2000： 277-296.

GONG Xiaonan. Manual of ground treatment＼[M＼]. 2nd ed. Beijing： China Architecture and Building Press， 2000： 277-296. （In Chinese）

[4]钱家欢， 帅方生. 边界元法在地基强夯加固中的应用＼[J＼]. 中国科学：A辑， 1987， 29（3）： 329-336.

QIAN Jiahuan， SHUAI Fangsheng. Application of BEM to dynamic consolidation＼[J＼]. Scientia Sinica： Series A， 1987， 29（3）： 329-336. （In Chinese）

[5]CHOW Y K， YONG D M， YONG Y， et al. Dynamic compaction analysis＼[J＼]. Journal of Geotechnical and Geoenvironmental Engineering， ASCE， 1992， 118（8）： 1141-1157.

[6]孔令伟， 袁建新. 强夯的边界接触应力与沉降特性研究＼[J＼]. 岩土工程学报，1998， 20（2）： 86-92.

KONG Lingwei， YUAN Jianxin. Study on surface contact stress and settlement properties during dynamic consolidation＼[J＼]. Chinese Journal of Geotechnical Engineering， 1998， 20（2）： 86-92. （In Chinese）

[7]李裕春， 时党勇， 赵远. ANSYS10.0/LSDYNA理论基础与工程实践＼[M＼]. 北京： 中国水利水电出版社， 2006：11-23.

LI Yuchun， SHI Dangyong， ZHAO Yuan. Basic theory and engineering practice of ANSYS10.0/LSDYNA＼[M＼]. Beijing： Water Conservancy and Hydropower Press， 2006：11-23. （In Chinese）

[8]THILAKASIRI H S， GUNARATNE M， MULLINS G， et al. Implementation aid for dynamic replacement of organic soils with sand＼[J＼]. Journal of Geotechnical and Geoenvironmental Engineering， 2001， 127（1）： 25-35.

[9]于德水. 强夯法处理湿陷性黄土路基的试验研究及数值模拟＼[D＼]. 北京： 北京科技大学， 2005： 48-76.

YU Deshui. The experimental study and numerical simulation of the collapse loss ground improved with dynamic consolidation method＼[D＼]. Beijing： University of Science and Technology Beijing， 2005： 48-76.（In Chinese）

[10]杨建华. 碎石土高路堤的强夯处理与沉降特性研究＼[D＼]. 武汉： 武汉理工大学， 2008： 75-97.

YANG Jianhua. Settlement property and dynamic compaction improvement parameters study for high crushed stone embankment＼[D＼]. Wuhan： Wuhan University of Technology， 2008： 75-97.（In Chinese）

[11]张建辉， 杨培轩. 强夯法处理风成砂地基的数值模拟分析＼[J＼]. 河北大学学报：自然科学版， 2013， 33（2）： 127-133.

ZHANG Jianhui， YANG Peixuan. Numerical simulation analysis of treating aeolian sand foundation using dynamic consolidation method＼[J＼]. Journal of Hebei University： Natural Science， 2013， 33（2）： 127-133.（In Chinese）

[12]何文， 钟志华. 显式有限元技术在车身薄壁梁结构件数值模拟中的应用＼[J＼]. 湖南大学学报：自然科学版， 2005， 32（1）： 1-5.

HE Wen， ZHONG ZhiHua. Application of explicit finite element technology to numerical simulation of body thinwalled structure＼[J＼]. Journal of Hunan University： Natural Sciences， 2005， 32（1）：1-5.（In Chinese）

[13]彭建兵， 陈立伟， 邓亚虹， 等. 车辆动荷载作用下黄土暗穴对路基稳定性影响的数值分析＼[J＼]. 中国公路学报， 2006， 19（4）： 17-22.

PENG Jianbing， CHEN Liwei， DENG Yahong， et al. Numerical analysis of subgrade stability influenced by hidden holes in loess under dynamic load of vehicles＼[J＼]. China Journal of Highway and Transport， 2006， 19（4）： 17-22.（In Chinese）

[14]何长明， 邹金锋， 李亮. 强夯动应力的量测及现场试验研究＼[J＼]. 岩土工程学报， 2007， 29（4）： 628-632.

HE Changming， ZOU Jinfeng， LI Liang. Field tests on measurement of dynamic stress of dynamic compaction＼[J＼]. Chinese Journal of Geotechnical Engineering， 2007， 29（4）： 628-632. （In Chinese）

[15]郭乃正， 邹金锋，李亮. 大颗粒红砂岩高填方路基强夯加固理论与试验研究＼[J＼]. 中南大学学报：自然科学版， 2008， 39（1）： 186-189.