模拟气动力

2024-07-26

模拟气动力（精选十篇）

模拟气动力篇1

现代体育竞技中的标枪项目有不少专家研究过其气动力, 在进行标枪气动力测量过程中, 蔡国华[1]、王倩[2]采用直立天平安装方式研究了标枪的气动力特性。这种方法对标枪的气动力的测量数据有直接影响, 虽然周继和, 石玉琴[3]用水平天平安装方法, 但他们实验方法共同的缺点是:将天平放置在标枪内腔。这样做过实验的标枪结构被破坏, 虽然人们掌握了它的气动力特性, 然而这些标枪实验结果也将不可使用。

实验方法亟待改进使我们考虑到可以用数值仿真来消除实验的弊端。起始于20世纪70年代的计算流体力学CFD与数值方法相结合, 对真实流场的模拟已日趋完善, 已经成为现代各种工程流体计算中不可或缺的一部分。一方面数值计算具有耗时短、经济、可重复的特点, 可大大加快研究进度;另一方面采用数值计算可以获得许多在实验中难以获得的数据, 为实验提供有力地支撑, 从而有力地推动研究向更深层次发展。

虽然国内在空气动力学数值计算起步较晚, 但从20世纪90年代末起进展迅速。大都来自自己所编写的程序, 或是其他商用CFD软件。而FLUENT作为一款商用CFD软件, 其数值模拟的应用已成为工程流体的主导。本文以FLUENT软件为平台, 以1:1的标枪模型[4,5]为研究对象, 用无黏流模型对其进行数值模拟计算, 求解标枪的气动力, 分析标枪的最佳出手攻角。

1 模型建立和网格划分

根据张绪树标枪模型的几何参数[4], 利用GAMBIT的多点输入法进行几何构建, 见图1。使用圆柱形计算区域, 见图2。空气场的半径为3 m长为6 m。由于标枪的长细比达到72, 所以在标枪近壁处采用四面体网格划分, 远壁处采用六面体网格划分。所选的空气场为Invalid无黏流模型, 空气为标准大气{T=288.16 K;ρ=1.225 kg/ m3}。

2 数值模拟计算

2.1 求解器和离散

采用分离显式求解器, SIMPLEC压力速度耦合格式, 压力项采用PRESTO (pressure staggering option) 交错压力格式, 动量方程采用Second Order Upwind二阶迎风格式。

2.2 边界条件

本模型标枪为固定边界, 而以改变入口风载的角度来考虑标枪的攻角。因此求解区域前端和圆柱桶壁入口边界为速度边界条件, 求解域后端为压力出口边界条件P=101 325 Pa。

2.3 计算结果修正及分析

在风洞实验仿真时, 计算出了攻角-30°—30°, 每隔5°时的收敛值, 从数值计算结果中提取所需气动特性, 并对数据进行曲线拟合, 与前人的实验数据进行对比, 趋势基本一致。而与王倩[6]的数据更为吻合。

2.3.1 计算结果修正

风洞实验洞壁干扰的修正方法有经典映像法、壁压信息奇点法、单参数线化壁压信息法等方法[7]。本文选用经典映像法对洞壁干扰进行了修正。

升力 $F_{L} = C_{L_{c}} \frac{1}{2} ρ v^{2} S$ ;

阻力 $F_{d} = C_{d_{c}} \frac{1}{2} ρ v^{2} S$ ;

俯仰力矩 $Μ = C_{m}_{c} \frac{1}{2} ρ v^{2} S$ ;

阻力系数CDc=CDu (1-2ε) +ΔCD+ΔCDW;

升力系数CLc=CLu (1-2ε) +ΔCL;

俯仰力矩系数Cmc=Cmu (1-2ε) +ΔCm;

以上各式的下标c表示修正后的值, 下标u表示未修正的测量值。右端各量可以用下列公式计算:

(1) 阻塞干扰计算

阻塞干扰因子ε可分解为两部分ε=εS+εW。εW是模型尾流阻塞干扰因子, εS是模型实体阻塞因子。其中 $ε_{S} = Κ τ \frac{V}{A^{3 / 2}}$ 。式中V—模型体积;A—试验段横截面积;K—与模型相对厚度或细长比有关的模型形状系数, 可由文献图中查得;τ—取决于风洞试验段截面形状及模型翼展与试验段宽度之比的系数, 可由文献图中查得。有时可近似认为K和τ的乘积等于E, 旋成体取E=0.96。小迎角流线体模型尾流阻塞干扰因子εW的计算公式 $ε_{W} = \frac{S}{4 A} C_{D 0_{U}}$ 。式中S—力系数参考面积;CD0u—测量的最小阻力系数。

(2) 下洗干扰计算

三维试验 $Δ a = 57.3 (1 + τ_{2}) δ \frac{S}{A} C_{L_{u}} (1 - 2 ε)$ ;

$Δ C_{D} = (1 + τ_{2}) δ \frac{S}{A} C_{L}^{2} (1 - 2 ε)^{2}$ ;

ΔCDW=-εSCD0u;

$Δ C_{L} = - 57.3 τ_{2} δ \frac{S}{A} C_{L_{u}}^{a} C_{L_{u}} (1 - 2 ε)$ ;

ΔCm=ΔCmW+ΔCmT;

$Δ C_{m W} = 57.3 \frac{l_{Κ}}{b} τ_{2} δ \frac{S}{A} C_{L_{u}}^{a} C_{L_{u}} (1 - 2 ε)$ ;

$Δ C_{m Τ} = \frac{\partial C_{m}}{\partial a_{t}} Δ a_{t}$ ;

$Δ a_{t} = 57.3 τ_{2} (l t) \cos a δ \frac{S}{A} C_{L_{u}}$ 。

δ—升力干扰因子, 可由文献图中查得;b—机翼平均气动弦长;τ2—小展弦比的机翼的上洗增量因子, 大展弦比机翼取τ2=0, 小展弦比机翼由文献图中查得;lK—力矩参考点至机翼1/2弦长处的距离;τ2 (lt) —可由文献图中查得;一般估算时可取

$\frac{\partial C_{m}}{\partial a_{t}} = - 0.8 \frac{S_{t} l_{t}}{S b} \frac{0.1 λ_{t}}{(2 + λ_{t})}$ 。

2.3.2 阻力特性分析

图3是“皇冠”牌男子标枪的在风速为30 m/s不同攻角下 (α) 的阻力 (Fd) 变化曲线图。可以看出阻力值均为正值。在攻角为-10°—10°的范围内阻力值的增加较为缓慢, 而其他范围却随攻角角值的增大而迅速增大。图形趋势与王肇明[8]、蔡国华[1]、周继和[3]等人的规律变化是一致的。

2.3.3 升力特性分析

图4是该标枪的升力 (FL) 与攻角 (α) 的关系曲线。在-30°—30°的范围内升力呈现反S形。随着攻角的增加升力值逐渐增大。在攻角0o时为升力值趋于零。结合阻力曲线, 在同一攻角下, 标枪所受的升力和阻力的绝对值是不一样的。阻力肯定缩短标枪的水平飞行距离, 那么升力就需要具体分析, 升力太小, 飞行不高, 达不到延长飞行距离的目的;升力太大, 飞的比较高, 然而阻力增大, 同样对延长飞行距离不利。因此, 下面引出升阻比特性的分析。

2.3.4 升阻比特性

图5是该标枪的升阻比 (Q=FL/Fd) 与攻角 (α) 的关系曲线。当攻角为零度, 阻力最小, 此状态的升力也很小, 趋于零。小阻力对飞行标枪减缓速度下降是有利的, 但无延长飞行时间的作用, 同轴投掷不是一种最好的投掷方法。因此, 我们不能孤立地去看一支标枪的升力或阻力, 而是应当用升力和阻力之比值, 也就是空气动力效率有称滑翔数来做为评定标枪气动性能的指标。从实验结果看, 标枪的升阻比在10°—20°间达到了升阻比在2以上。而以15°为最佳。在这一点上周继和测出25°为男子标枪的最佳升阻比, 而王倩推算的出手攻角为20°, 王肇明为-10° 。本文与王倩所测的“皇冠”牌男子标枪的升阻比最值的最佳攻角值一致。综上, 本人认为各品牌标枪有其各自的气动力特性, 不能一概而论男子标枪或女子标枪的最佳攻角值。而应相应的对各标枪进行真正的认识。

2.3.5 力矩特性分析

图6是力矩 (M) 与攻角 (α) 的关系曲线。可以看到在-10° —10°之间, 力矩值随攻角的增大而增大。其他区域则迅速向相反方向增加。而不是运动生物力学课本[9]的线性关系。结合升阻比考虑标枪出手时的最佳攻角, 本人认为是15°—20°。在此范围内标枪被作用了产生低头的力矩, 而升阻比的比值在10°—15°还是较高, 因此可以达到延长标枪在空中的飞行时间的目的。

3 小结

(1) 该“皇冠”牌男子标枪阻力与攻角的关系曲线近似对称, 攻角为0°时阻力最小为0.24 N。升力随攻角呈现反S形, 攻角为0°时升力趋于0 N。结合升阻比, 力矩与攻角的关系曲线图分析得知:标枪出手时的最佳攻角为15°—20°。

(2) 只要通过常规方法测定标枪的外形与重心, 那么就可以通过模拟仿真了解该标枪的气动力特性, 避免对标枪气动力特性的不可逆的测量, 而且可重复性强。

(3) 标枪在飞行过程中是弹性体, 而模拟过程中视标枪为刚体, 虽然对标枪的飞行时真实性还需继续努力, 但对运动员投掷时初始攻角的选择具有指导意义。

(4) 要发挥低重心标枪的气动特性, 提高运动成绩, 必须有适宜的投掷出手条件与之匹配, 建议作六自由度 (6DOF) 模拟投掷, 给运动训练提供理论依据, 再在运动实践中检验。

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

关键词：标枪,气动力,数值模拟

参考文献

[1]蔡国华.标枪的气动力特性和几何、物理参数测量.空气动力学学报, 2004;22 (4) :452—456

[2]王倩.关于女子标枪运动员的实际最佳出手角度和初始攻角的研究.北京体育学院学, 1990;2:30—35

[3]周继和, 石玉琴.低重心标枪的空气动力学研究.中国体育科技, 1990;06:61—66

[4]张绪树.标枪出手初始参数与投掷成绩的关系研究.太原理工大学硕士学位论文, 2003

[5]廖红, 张绪树, 陈维毅, 等.考虑空气阻力时标枪最佳出手角度研究.中国体育科技, 2007;43 (1) :98—101

[6]王倩.新规则男子标枪气动力风洞实验研究.北京体育学院学报, 1989;4:78—82

[7]张厚梅.风洞实验干扰与修正, 北京:国防工业出版社, 2003;77—155

[8]王肇明.标枪动力学的优化计算.力学与实践, 1983;4:40—42, 54

开放式气动力数值模拟系统研究篇2

开放式气动力数值模拟系统研究

气动力数值模拟系统是CFD流场解算技术、网格生成技术、数据可视化技术和网络技术相结合的产物.描述了一种运行于计算机网格环境下多任务、多用户、交互式的数值模拟软件环境-气动力数值模拟系统,该环境能计算模拟亚、跨、超音速绕流的.气动力,也能方便地集成模拟不同复杂气动布局飞机流场的解算程序.介绍了其网格生成、数据可视化、系统接口等核心部分的结构和功能.

作者：任碧宁白文朱培烨周天孝 REN Bi-ning BAI Wen ZHU Pei-ye ZHOU Tian-xiao 作者单位：中国航空计算技术研究所,陕西,西安,710068刊名：航空学报 ISTIC EI PKU英文刊名：ACTA AERONAUTICA ET ASTRONAUTICA SINICA年，卷(期)：20(6)分类号：V211.3关键词：CFD 网格生成接口技术数据场可视化

模拟气动力篇3

摘要：基于显式动力非线性有限元分析方法，利用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左右土体发生隆起的现象.

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.

风力机叶片气动特性数值模拟篇4

风能是一种无污染的可再生能源, 与传统能源相比, 风力发电具有发电成本稳定, 无碳排放以及可利用风能分布广泛等优点。近年来, 世界各国对风能利用开发越来越重视, 这也促进了风力发电技术的迅速发展。风力机叶片是风力发电机组的核心部件之一, 其设计制造质量直接影响着机组内部的其它部件, 而叶片翼型的空气动力学特性又是叶片设计的关键[1,2]。在叶片设计的过程中, 由于实验条件的限制, 要想获得翼型的气动数据比较困难, 这时就需要对翼型进行数值计算和仿真模拟。本文采用FLUENT软件研究NACA4412翼型在不同攻角下的气动性能。

1 湍流模型的选择

湍流模型的数值模拟有很多方法, 而在航空航天领域中, S-A模型因其主要着眼于壁面边界流动的研究而迅速发展。该模型是一种有效的低雷诺数湍流模型, 并且能够对边界层受黏性影响的区域进行恰当求解。对动态流场进行计算时, 由于流动复杂并且计算量大, 若采用双方程模型, 计算复杂且不易收敛, 若采用零方程模型, 计算精度又会受到影响, 而S-A单方程模型相较于双方程模型和零方程模型具有易收敛、求解快等优点, 所以本文选择S-A单方程模型作为湍流模型。

为求解湍流黏度S-A模型给出如下形式的输运方程[3]:

2 计算模型的建立

本文选择相对厚度较小、气动性能较好的翼型NACA4412作为研究对象。在GAMBIT中对翼型周围流场进行网格划分, 网格质量对计算结果及计算效率均有影响, 因此应合理布置网格点。由于翼型周围空气流动变化及前后缘处梯度变化均较大, 网格应密集, 而远场边界处流场梯度接近于0, 所以当逐渐接近远场边界时网格应稀疏。综合考虑计算精度及计算时间, 在翼型上共布置150个节点, 其中尾流方向布置60个节点, 翼型表面法向布置90个节点, 最终得到整个计算物理域长度为32.5倍翼型弦长, 宽度为25倍翼型弦长。

选择S-A湍流模型, 关闭能量方程, 空气密度设为1.225 kg/m3, 边界条件设置上、下及左半圆边界为速度进口, 右边界为压力出口, 并且根据前方来流速度及攻角设置x方向和y方向的速度。

3 翼型二维仿真计算结果及分析

翼型确定以后, 要选择合理的设计点, 使其在该点附近有较大升力及较小阻力。因此, 翼型数据应该选升阻比最大的点所对应的攻角、升力系数、阻力系数等参数值。

设定雷诺数Re=106, 来流风速为12 m/s, 模拟NACA4412翼型攻角为0°~20°下的二维流场, 得到以同攻角的表面压力分布, 速度分布, 升力、阻力系数, 升阻比。

图1是攻角为0°、5.5°、10°时叶片附近压力云图。图2是攻角为0°、5.5°、10°时叶片附近速度矢量图。通过叶片附近压力云图可以看出叶片下表面压力高于上表面压力, 两者的压力差为叶片提供了升力。将压力云图和速度云图进行对比可以看出, 流场压力高的地方气流速度较小, 流场压力低的地方气流速度较大。并且从图中可以看出攻角为5.5°时叶片上下表面的压力差最大, 即升阻比最高。

图3是升力系数、阻力系数随攻角变化曲线。FLUENT软件采用数值模拟的方法对翼型气动特性进行仿真, 而Profili软件根据儒可夫斯基升力定理对其进行解析计算。经验证, 两种计算方法得到的升、阻力系数随攻角变化趋势基本一致。

通过升、阻力系数可以分别计算出各对应攻角下的升阻比。升阻比随攻角变化曲线如图4所示。从图中可以看出当攻角为5.5°时对应的升阻比最大, 所以最佳攻角为5.5°。

4 结语

1) 通过FLUENT软件在定雷诺数和定风速的情况下, 对NACA4412翼型在不同攻角下进行气动特性数值模拟, 得到攻角为0°、5.5°和10°的叶片附近压力云图和速度矢量图, 为后续叶片的气动设计提供参考。

2) 在定雷诺数和定风速的情况下, 得到不同攻角下的升、阻力系数。通过升、阻力系数计算升阻比, 并根据最大升阻比确定了最佳攻角。

参考文献

[1]牛山泉.风能技术[M].北京:科学出版社, 2009.q

[2]赵振宙, 郑源, 高玉琴.风力机原理与应用[M].北京:中国水利水电出版社, 2011.

气动塞式喷管实验与数值模拟研究篇5

气动塞式喷管实验与数值模拟研究

提出了塞锥型面设计的特征线法及示例结果,得出了侧喷管粘性跨音速和超音速内流场数值计算结果,及塞锥外流和尾流的`Euler方程解,总结了塞式喷管工作特性的固体推进剂燃气模拟实验系统、实验结果及主要结论.塞式喷管火箭发动机是一个复杂而有挑战性的研究方向,对于发展新一代先进天地往返运输系统有重要作用.

作者：刘宇张正科程显辰张振鹏 Liu Yu Zhang Zhengke Cheng Xianchen Zhang Zhenpeng 作者单位：北京航空航天大学,宇航学院刊名：北京航空航天大学学报 ISTIC EI PKU英文刊名：JOURNAL OF BEIJING UNIVERSITY OF AERONAUTICS AND ASTRONAUTICS年，卷(期)：25(6)分类号：V430 V433.92关键词：实验流场模拟特征线法塞式喷管

模拟气动力篇6

关键词：流体动力学仿真原理；有限元分析技术；铝合金；动静盘；重力铸造；应用

中图分类号：TP391.9 文献标识码：A 文章编号：1006-8937（2016）27-0041-02

1 概述

在涡旋式汽车空调压缩机铝合金动静盘毛坯件重力铸造领域，大多数企业都是一模出一件铸造件，生产工艺极其落后，效率低，浪费大。为了提高效率，降低劳动力成本，许多企业在探索一模出四件或多件铸件，并进行模具的创新设计和开发。然而，“一模四腔”新模具，设计复杂，造价昂贵，且无前置经验可借鉴，亦无其他设计依据，设计失败率很高；“模组设计→加工制造→试模→修改”，整个过程无科学依据与分析，只能靠不断的资金投入，反复的试验、修改，不断的“尝试”。目的是为了成功，但资金投入大，制造、试验成本高，成功率低。

流体动力学仿真&有限元分析技术，在其他行业发展迅速，如果将其与模具开发结合，可以大大提高设计研发效率和试验、开发的成功率，降低开发成本。

2 涡旋式压缩机动静盘生产工艺分析

2.1 涡旋式压缩工作原理

涡旋式压缩机是由两个涡旋体相互成180o的角，装配在一起。其中一个被固定，即静盘；另一个旋转，即为动盘。这一对渐开线形的涡旋体，组成3对同时工作的压缩腔，动盘一方面沿着很小的偏心距（曲轴回转半径）轨道移动（即摆动），另一方面与静盘接触做相对转动，与静盘形成3个变容积的密封腔，每1转中，第1个腔在吸气，第2个腔在压缩气体，第3个腔在排出气体，也就是说，在每个转角，都在持续循环地进行吸气、压缩、排气，因此没有负荷的起伏变化，所以涡旋式压缩机运转非常平稳，这种特性对发动机非常有利。

2.2 一模四腔的浇注装置及模具存在的问题

目前，涡旋式汽车空调压缩机铝合金动静盘生产，主要的几种制造工艺分别为：重力铸造、低压铸造、液态模锻、热模锻和背压成型工艺，其中重力铸造在规模化生产中应用较广，但是，铝合金动静盘重力浇注，普遍采用一模一腔即：一个模具生产一个动盘或者静盘，效率低，浪费大，成本高；一模四腔的浇注装置及模具，尚在研发试用阶段，存在很多问题和技术难题，如：“一模四件铸件”铸造模具，设计复杂，造价昂贵，且无前置经验可借鉴，亦无其他设计依据，设计失败率很高；浇冒口、留道的尺寸、流量、流速的确定；冒口铝水提前结晶、远端反流补缩效果差；试验过程中存在充型不满、有热结、热裂纹、缩孔、砂眼、夹渣等缺陷。

3 一模四腔模具设计完成后的动力学仿真模拟&有限元分析

3.1 “动盘&静盘铸造方案”流体力学仿真分析

通过对一模四腔模具设计进行动力学仿真模拟根据仿真结果，如图1所示。我们总结出了如下改善方案以更好地改善一摸四腔模具的产出质量。

3.2 铸造方案改善

毛胚顶部表面疏松&缩孔：①增加浇口与冒口高度（浇铸量）；②顶部加强筋疏松部位增加铸造余量，并增加排气功能；③请供应商严格管控设计精准度，保持毛胚件水平，并将浇道详细设计方案提供给研发小组评估；④铸件缺少编号，请供应商纠正设计。

3.3 相关配套改善

避免冒口补浇，降低人力作业强度；模具加热与冷却必须自动化；模具配套需要符合环保要求，增加安全光栅，急停开关等设备，避免发生事故。

4 动力学仿真模拟&有限元分析在“一模四件”铝合金动静盘重力浇注中应用

“一模四腔”重力铸造铝合金动盘模具，设计、仿真模拟、修正后，随即在工厂进行了涡旋式动盘重力浇注试验，经过2次试验，成型毛坯件（一模出4件铸件）终于成功铸造出来。但毛坯件呈现以下缺陷：

①动盘涡旋根部、轴承孔底部出现了不同程度的热裂纹；

②基圆底平面有凹坑；

③帽口较大、较厚不利于后步工序敲掉帽口，甚至会在敲帽口时将动盘基体敲坏（俗称“掉肉”）。

④浇注流道与动盘体联接处凸出，有较大的台阶，车削加工时不便于装夹，且易打坏刀具。

⑤有热结、热缩现象。

针对上述存在的问题和出现的缺陷，我们应用动力学仿真模拟&有限元分析原理，进行缺陷原因分析，并提出改进措施。

5 技术改进效果

①完成&完善。毛胚件（出模整件）之3D图纸，如图2所示。

②利用流体动力学仿真软件对毛胚件做有限元分析。如图3所示。

③按照需求（如表面品质，重力分析，热节分析等），按照环境要求（如材料特性，环境温度，充型速率等）进行仿真，制作成铸造仿真分析动画。如图4所示。

④ 根据仿真分析结果，修改毛胚件（出模整件）设计，优化模具设计，如图5所示。

由图5可以看出样品正面和背面分别出现冷隔和裂纹现象，通过仿真分析我们得知出现冷隔的原因主要有：模具温度和铝液温度较低；排气不畅；冒口的浇口小；浇铸速度太慢或浇铸中断所致。

建议方案：提高铝液与模具温度；适当提高浇铸速度；适当增大冒口浇口。

出现裂纹的主要原因有：铸件有尖角或厚薄过渡不圆滑；局部温度过高；未凝固即开模；缩补不足。

建议方案：检查毛胚图纸，倒圆角，调整厚薄；注意涂料厚度；增大冒口补缩；考虑是否放置冷铁。

⑤“设计-仿真-优化设计”过程中，仿真为主导，优化设计，提升成本优势与设计成功率.

6 结语

综上所述，通过在涡旋式汽车空调压缩机铝合金动静盘毛坯件重力铸造模具设计中，应用流体动力学仿真&有限元分析技术，可以生动形象地观察出实验过程中存在的主要问题，并找出相应的解决方案，大大提高设计研发效率和试验、开发的成功率，降低开发成本。

参考文献：

一维纳米摩擦分子动力学模拟篇7

世界上使用的能源大约有1/3～1/2消耗于摩擦, 如果能够尽力减少无用的摩擦消耗, 便可大量节省能源。另外, 机械产品的易损零件大部分是由于磨损超过限度而报废和更换的, 如果能控制和减少磨损, 则既减少设备维修次数和费用, 又能节省制造零件及其所需材料的费用[1,2]。作为研究相对运动的作用表面间的摩擦、润滑和磨损, 以及三者间相互关系的理论与应用的一门边缘学科的摩擦学自古至今都得到了相应的重视。摩擦学的研究得到了相应的摩擦学库仑定律和Amoton规律, 这些经验规律能够合理地解释生活生产中的摩擦现象, 对于指导人类日常生活和社会生产的进行起到了重要的推动作用;但是本世纪新兴的纳米技术和纳米材料研究的不断深入, 使得人们有能力制造出尺度越来越小, 而性能越来越完善的微型装置, 尤其在备受关注的生物、环境控制、医疗器械、航空航天、数字通讯、传感技术等前沿领域快速发展, 纳米机械的发展使得摩擦学问题再次成为科学研究热点。因为对于纳米尺度的物质系统, 基于连续介质假设的传统理论可能不再有效, 纳米尺度下的摩擦学系统由于纳米材料自身的特性将表现出不同的摩擦学规律, 而“纳米摩擦学”即是为解决微纳米系统中摩擦的规律和控制而逐步形成的。

近年来纳米摩擦学的发展较为迅速, 主要基于两方面原因:一方面是纳米摩擦测试技术的进步使得对纳米摩擦现象的实验测试成为可能。在实验室中可以通过原子力显微镜[3], 摩擦力显微镜[4], 石英晶体微天平[5,6]等对纳米摩擦过程进行观测, 从而能够更好地理解摩擦的起源问题;另一方面是计算机计算能力的快速增长和计算软件的不断更新。计算能力的快速增长使得复杂的摩擦现象可以通过分子动力学进行相对真实的模拟再现, 从而为深层次地理解摩擦现象提供了新的研究手段[7,8,9]。

本文将经典FK模型和能量耗散机制相结合开展类Xe原子在类Ag表面滑动摩擦的分子动力学模拟[10]。为了得到在一维纳米摩擦的一般规律, 而不强调求解某种材料的摩擦特性, 所以在模型中设定相应的参数调节来得到不同的摩擦环境。将Xe与Ag的相互作用势能改变为一种组合周期势能进行描述。模拟过程中原子在受到不同大小附加拉力的作用下, 在不同温度和不同匹配度的衬底表面进行滑移。模拟过程除了设置了界面势能大小, 匹配度等参数 (图1) , 还设置复合周期势能中正弦势能振幅比例UA, 势能周期比例UP两个重要的参数。振幅比例UA, 势能周期比例UP的调节可以很好地体现出一般衬底晶格中含有的两种原子之间的相对距离和两原子与滑移原子作用相对强弱对纳米摩擦过程的调制作用。

本文的第二部分描述分子动力学的模型和模拟过程中的参数设置。模拟的相关结果和讨论在文章的第三部分给出。第四部分对文章做出了总结。

2、模拟模型

本文模拟无限长的一维原子链在调节后的无缺陷原子表面的滑动过程来获得纳米摩擦学匹配度对摩擦的调节作用。滑移原子之间的相互作用采用Lennard-Jones势[11,12,13], 如下式:

其中, (i, j为标记不同的原子) ;0为T=0时的原子间的平衡距离;ε为势能阱深度。计算过程中截断到第三近邻原子, 晶格常数取α=0.99720[14]。

附加表面原子与衬底原子之间的相互作用以U (图1) 来表示, 由两个正弦势能复合而成如下式:

U (x) 是一个关于坐标的周期势函数, 周期为b, b=b1+b2;u1和u2分别为正弦周期势能的振幅的一半;b1和b2分别为两正弦周期势能的周期;定义势能的周期比参数, 即:b1= (1-UP) b, b2=UPb;势能振幅比为, 即:u1=u0, u2=UAu0 (其中u0为势能设置参数) 。

附加原子还受到一个复合高斯分布的随机力i, 随机力通过涨落-耗散机制与电子摩擦系数ηe相关[15], 且满足下式:

其中, 尖括号表示对时间求平均值;σ (t) 表示delta函数;k为Boltzmann常数;T为系统温度。根据涨落-耗散方程选取随机力i可以使系统温度保持在设定的温度。

在模拟过程中, 对所有的附加原子附加以恒定的拉力F, 经长时间驰豫系统进入稳定状态后求得质心速度, 由公式确定摩擦系数[14]。

附加的原子应满足如下的Lagevin方程[16,17,18]:

模拟采用对分子运动的数值积分得到不同状态下的摩擦过程。在模拟计算中, 为减少计算量, 提高效率, 对计算公式进行无量纲化处理。选用0, ε和m作为基本单位, 将其他单位进行换算, 例如时间单位变为。另外, 为了减少模拟尺寸的影响, 在模拟过程中采用了周期性边界条件, 所以在有限长度内附加原子数Nα和衬底原子数Nb之间地有如下关系:

其中, 比例常数C称为附加原子与衬底原子的匹配度。

为了得到一个稳定状态, 首先在附加拉力F=0的情况下, 时间步长△t=0.01t0, 对初始条件下的Lagevin方程进行模拟, 模拟过程进行30万步。模拟结果系统达到一个没有人为干扰的自然稳定的初始状态。然后对附加原子链中所有的原子附加拉力F进行模拟, 直到系统进入另一个稳定状态。在平衡状态下求解质心速度, 从而求出系统的摩擦系数。整个模拟过程设定在一定温度T下, 通过调节两势能周期比UP和振幅比UA, 附加拉力F等参数得到不同的模拟环境, 了解不同情况下的匹配度C对纳米摩擦的影响。

3、结果与讨论

首先, 给出在UA=UP=1极限情况下的摩擦模拟结果, 进行程序的正确性验证;然后计算不同势能周期比例和振幅比例随着匹配度变化纳米摩擦的变化。

3.1 一种极限情况UA=UP=1 (普通的正弦势能形式)

由前文可知, 当UA=UP=1复合准周期势能形式将变为一种简单的正弦形式, 也就是经典纳米摩擦分子动力学模型中所用的势能形式。

图2给出了附加拉力与质心速度关系模拟结果与文献[14]所得到的结果 (c, d为本文得到的模拟结果) 。通过对比发现本文的计算方法得到的摩擦系数会比文献[14]的要稍微大一些 (最大误差小于10%) , 但两种计算方法对摩擦系数的预测的趋势都是一致的。之所以出现误差的主要原因是在求解Lagevin方程时采用了不同的计算方法。但是因为本文的重点在于找出在复合势能形势下滑动摩擦所表现的规律, 所以两种计算结果的不同不会影响到摩擦规律的发现。

3.2 一般情况

为了得到纳米摩擦在一般情况下的规律性, 模拟过程中设置了相关可调参数, 如表1所示。

在KBT=0.1ε, 势能u0 (0.01ε～0.04ε) , UP (0.03～0.09) , UA (0.03～0.09) 情况下, 模拟C在 (0.03～0.09) 之间对摩擦系数的影响。图3给出了不同势能u0大小, 不同UP, UA情况下得到的摩擦系数的倒数随C的变化曲线。

由图可以看出, 最显著的特点是在不同摩擦条件下的摩擦系数倒数曲线随C的变化在C=0.5处出现了汇聚且出现最小值, 而在两侧则离散增加。当C=0.5时出现摩擦系数倒数最小值, 也就是出现了摩擦系数最大值, 即在完全匹配的情况下出现摩擦最大值。完全匹配处之所以出现摩擦最大值的原因在于此时所有的附加原子步伐统一, 能够同时受到衬底原子的拉力而加速移动, 又能够同时受到阻力减速, 相互之间无法将能量转化为势能保存, 而只能以热量的形式进行耗散。曲线的聚合同时说明对于纳米摩擦无论势能形式如何组合变化, 势能大小的如何变化, 在对摩擦的调制作用中界面晶格匹配度C对整个摩擦过程的调制作用最为明显, 在完全匹配的情况下摩擦最大。

由图可以发现另一个特点, 摩擦系数在C为0.5两侧随C的增大或者减少而迅速减小, 但在C为0.8时出现了转折点。此转折点在不同的条件下, 对摩擦的调制作用表现出不同程度的影响, 说明此转折点为复合势能情况下产生的影响, 而这种调节作用与势能周期比例UP和势能振幅比例UA之间存在非线性关系。转折点的出现原因为附加原子在不同的衬底复合势能阱中进行运动时受到衬底原子的拉扯作用产生了强烈的声子振动所带来的影响, 这种拉扯作用主要受到势能振幅比例与势能周期之间的比例关系的影响。当处于衬底势能最高点时, 附加原子将出现失稳并自动跳跃到下一势能最低点, 然后在此平衡位置剧烈震荡和激发声子, 从而使能量不可逆地以声子的形式耗散掉。当附加原子在势能阱较深, 而周期较小的势能阱中运动时受到的衬底原子拉力由最大值迅速变成了反向阻力最大值的过程中, 附加原子无法将拉力所带来的能量完全转化为原子之间的势能保存, 而只能以晶格振荡的方式转化为热的形式放出, 从而导致系统拉力作功成为摩擦产生热量过程, 进而实现能量以晶格振动的形式耗散, 摩擦系数增大。

模拟结果表明, 无论如何改变势能周期和势能大小比例组合势能, 均可以发现摩擦系数受到了匹配度的调制作用最大, 说明摩擦的主要原因来源于衬底与附加原子之间的晶格不匹配造成滑动过程中能量的耗散。

4、总结

本文在不同势能组合形式下, 不同晶格匹配度情况开展了纳米摩擦分子动力学模拟, 得到了一维纳米摩擦的一般性规律。模拟过程中主要设定了势能周期的分配, 势能大小的分配等几个重要参数, 通过固定不同参数模拟了匹配度对摩擦过程的调制作用。通过模拟发现在不同的摩擦环境中, 匹配度对纳米摩擦的影响出现了两个极值点, 一个为匹配度0.5, 此时纳米摩擦的主要调制作用由匹配度来决定, 其他因素对摩擦的调制作用不明显;另一个为匹配度0.8, 此时纳米摩擦的主要调节作用为势能振幅比例和势能周期比例调制, 这种调制作用实现对摩擦的控制作用, 为多原子的表面膜的摩擦系数控制提供了参考依据。通过控制多原子薄膜中不同原子种类和它们之间的晶格常数来控制薄膜的摩擦系数, 从而实现对摩擦磨损的合理控制。

摘要：本文在建立一维纳米摩擦模型的基础上, 对纳米摩擦现象进行了分子动力学模拟研究。研究结果表明, 表面接触界面之间的晶格匹配度对纳米量级的摩擦起到较大的调制作用, 而这种调制作用随着晶格匹配度的变化而产生调制作用的大小也不同。分析表明在完全匹配的情况下摩擦最大, 而其他因素例如界面接触势垒的大小、形式等与匹配度相比对摩擦的影响较小。

分子动力学模拟及应用简介篇8

1.2 分子动力学模拟的基本过程。

分子动力学模拟的基本过程就是体系能量最小化和能量平衡的过程。能量最小化的过程实际上就是求势能函数极小值的过程, 利用力场函数计算分子内各种可能结构的势能, 搜索体系的势能最低点。能量最小化是对所有模拟体系探讨的第一步, 即用系统的方法改变体系的原子坐标, 通过分子内部运动的势能下降的办法, 使分子结构逼近稳定结构, 避免不合理的结构在模拟中出现错误。

只有当体系完成能量最小化后才能进行真正的分子动力学模拟, 即能量平衡过程, 能量最小化的过程是求势能函数极小值的过程, 能量平衡的过程实际上就是解这个多粒子体系的牛顿运动方程来描绘体系中所有粒子随时间演化的规律的过程。这时就可以根据粒子的位置和动量计算任何可观测物理量的稳定值。在分子动力学模拟分析中, 通常通过计算粒子运动的方均根偏差 (root mean square deviation, RMSD) 用来衡量体系相对于初始时刻是否能达到平衡, 同时也可以反应原子偏离比对位置的程度。RMSD值越大, 原子运动的空间范围越大, 原子的空间位阻也就越小, 原子α在时间Ni内RMSD的求算公式:

2 分子动力学模拟常用软件

分子动力学模拟常用的软件有:NAMD、GROMACS、CHARMM、AM-BER等, 这些分子动力学软件应用范围都比较广泛, 但是每个软件自身有具备各自的相对优势[2]。接下来重点介绍常用分子动力学分析软件, NAMD和GROMACS.NAMD是由美国伊利诺斯大学理论与计算生物物理研究组共同开发的分子动力学模拟软件, 适用于计算生物大分子体系和化学材料体系。NAMD并行计算效率非常高, 可以有效并行上千个处理器, 兼容使用Amber、X-PLOR、GROMACS等多种力场以及多种输入和输出文件格式。NAMD几乎支持所有的操作系统, 同时还可以与分子可视化与辅助分析软件VMD配合使用, 用来分析体系原子的运动轨迹。NAMD的最新版本可以从http://www.ks.uiuc.edu/Research/namd/免费下载。GROMACS是由荷兰的格罗宁根大学的生物物理化学系研发的分子动力学的模拟软件, 最初主要目的是为了研究脂类、核酸和蛋白质等生物分子的折叠以及生物分子内部的相互作用[3]。后来Gromacs本身进行了大量的算法的优化, 计算功能更加强大, 其在模拟几百到几百万粒子体系的牛顿运动方面具有极大优势, 适用的研究范围也扩展到了液晶、晶体、聚合物等领域。GROMACS的模拟程序包主要使用GROMACS力场, 结合分子动力学、随机动力学或者路径积分方法模拟体系中的任意分子的运动, 分析体系内部相互作用以及构象的变化。此外, GGROMACS提供轨迹的可视程序, 具有友好的用户界面, 方便用户的操作使用。同时还为轨迹分析提供了大量的插件工具, 用户不必再为常规分析编写任何脚本以及程序。GGROMACS还可以在Windows、Unix和Linux等多种操作系统上使用。Gromacs的最新版本可以从http://www.gromacs.org/免费下载。

3 分子动力学模拟的应用

分子动力学模拟经过几十年的发展, 适用的范围越来越广的同时, 计算的精度也得到极大的提高, 目前已成为生物物理、生物化学、凝聚态、材料、药物设计等方面不可缺少的研究工具。接下来简单介绍分子动力学模拟在研究蛋白质构象变化中的应用。氨基酸的种类, 排列顺序以及多肽链的空间折叠方式决定了蛋白质的特异性, 蛋白质结构基础决定了蛋白质的生物学基础。目前PDB (Protein Data Bank) 库已经明确了大量蛋白的结晶结构, 这为我们利用分子动力学的方法研究蛋白的空间结构、构象变化、相互作用以及运动趋势提够了结构基础。为什么要采用分子动力学模拟的方法研究蛋白质的空间结构、构象变化以及运动趋势等?目前, 蛋白质在行使其生物学功能的过程是一个连续的动态过程, 这是传统的电生理实验无法实现的, 同时分子动力学模拟的方法可以在分子原子水平来观察蛋白氨基酸之间的相互作用, 蛋白的构象变化以及运动趋势等。

我们以钙激活大电导钾通道 (BK channel) 为例, 来介绍分子动力学模拟在蛋白模拟方面的具体应用。首先从PDB库中获取BK通道的关闭态的晶体结构 (PDB:3NAF) , 因为3NAF中的部分氨基酸缺失, 首先利用同源模建的方法, 构建了完整的蛋白结构。接下来, 在完整的晶体结构的基础上构建物理模型, 以模拟蛋白在生物体内的真实生理环境, 具体步骤为:对晶体结构加水盒子, 模拟生物体内的液体环境;在加水体系的基础上添加一定物质的量的K+和CL-来中和体系内多余的电荷, 以模拟生物体内中性环境。经以上两步骤后, 就完成了物理模型的构建。在进行分子动力学模拟之前, 我们需要构建用来执行动力学模拟的配置文件, 配置文件中包含了一系列的配置参数, 配置文件的作用是设定NAMD分子动力学模拟的过程中的各种参数, 比如输出文件的频率、体系的温度、计算时长、PDB文件和PSF文件的存储路径、输入和输出文件的存储路径等。配置文件是分子动力学模拟计算的关键文件, 配置文件构建成功与否直接关系到模拟能否正常运行。接下来运行NAMD:在集群中的目标路径下运行配置文件, 根据构建的蛋白体系的大小, 模拟相同的步数, 其模拟计算时间也会有所差异, 模拟时间的长短受体系内原子数多少的影响。在分子动力学模拟结束后, 对模拟过程中产生的输出文件进行构象变化、蛋白稳定性和氨基酸相互作用的分析[4]。

目前, 我们已经利用分子动力学模拟完成了BK通道关闭态构象的模拟工作, 接下来主要是结合轨迹dcd文件, 分析模拟过程中, BK通道的构象变化、结构的柔性变化、结合位点位置氨基酸相互作用的变化、以及在模拟过程中BK通道的运动趋势等。

参考文献

[1]杨坤.一类非平衡分子动力学模拟中优化方法的研究[D].大连:大连理工大学, 2009.

[2]庞春丽.钙离子与跨膜蛋白TMEM16A通道静电结合机制的研究[D].天津:河北工业大学, 2014.

[3]Hess B, Kutzner C, Spoel D V D, et al.GROMACS 4:algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation[J].Journal of Chemical Theory&Computation, 2008, 4 (3) :435-447.

电网络模型的动力系统模拟应用研究篇9

当前, 随着工业及其自动化的迅猛发展, 模拟方法已经成为各种复杂系统研制工作的必不可少的手段。建立模型是系统模拟分析的基础, 其本质是依据系统之间的相似性原理。在研究不同类型的系统时, 一种系统可能比另一种系统更容易通过试验进行研究, 所以可以在一对系统之间建立一种对应关系, 从而可以利用模型对实际系统进行研究。这样不但可以使实验系统简化, 而且使实验成本和复杂及危险程度大大降低。

随着工业中各学科的不断发展, 学科的交叉和渗透也是当前科技辉煌发展的重要方向之一, 电路及点网络理论已经在机械力学、生物医学、工程热力学, 甚至在经济学等领域上已经初步得到了应用[1,2,3], 尤其在前两个学科上的应用已经日趋成熟[4,5,6,7], 如扩散-对流-化学反应的耦合和相应传输求解问题、网络的拓扑处理、三角化学反应问题、微观对称网络与混合网络、非平衡态热力学与运动学系统的时间特性与进化、分段线性网络、稳定与摇摆等。利用电路系统来模拟实现能源动力系统的分析和优化, 将大大提高对实际系统的分析效率和处理精度。

1电网络系统的物理构成法则

在电学中, 其基本计算处理元件包括电阻、电容及电感。它们的u-i或u-q关系式如下:

$\begin{array}{l} u = R i (1) \\ q = c u (2) \\ u = L d i / d t (3) \\ q = \int i d t (4) \\ L i = \int u d t (5) \end{array}$

从物理意义上讲, 电阻与力和流有关, 电容与电荷和力有关, 电感与冲量和流有关。对于理想化的电阻, 它的作用可以理解为通过网络元件的所表现出来的耗散。在网络的Langrangian公式中, 电阻代表了系统的非守恒思想。电容是一种没有耗散的能量储存形式, 表现为趋于平衡的孤立系统的理想状态。由于在电容的真实过程中总是存在耗散的, 因此现实电路中要求有一个电阻与电容串联以保证在数学上对系统的忠实描述。在大多数系统中, 电容还是系统提供动力的元件。电感是一个理想的惯性元件, 其模拟应用的主要领域是机械学。

2动力系统中的模拟应用

电网络在动力系统中模型的建立过程, 其实就是等效电系统的过程, 其对应原则遵循电路、电阻和电容的广义思想策略, 即在一切含有物流与能流的动态系统中, 不论该系统中是否有电的存在都可找到电系统中电阻、电容和电感等概念的对应物, 电阻可以对应各种阻力和对能量的消耗与节流, 电容可以对应各种形式的储存能量, 而电感则可对应各种惯性存储能量。这样所建立的网络模型很自然要服从Kirchhoff的电流电压定律, 而这两条定律是正交的, 正适合各种强度量与广延量之间的关系, 因此从理论上讲, 电网络模型可以模拟表达任意的力 (强度量) 与流 (广延量) 的耦合参数与信息状况。

2.1 弹簧机械系统模拟描述

机械学中存在的许多物理现象, 都可用对应的电方法描述, 进而为相应问题的求解和分析提供了一种可行的简洁途径。图1为一个由质量、弹簧、阻尼器组成的机械系统 (阻尼器的质量可以忽略不计) , 在对应的图2中, 描述为一个对应的由电阻、电感、电容组成的模拟电路系统[8]。

上述系统描述后所建立的模型分别为:

$\begin{array}{l} m \frac{d^{2} x}{d t^{2}} + f \frac{d x}{d t} + k x = Ρ (6) \\ L \frac{d^{2} q}{d t^{2}} + R \frac{d q}{d t} + \frac{1}{c} q = u (7) \end{array}$

可看出, 这两个方程具有相似的数学描述, 并且在参数上存在一一对应的关系, 相应地两者的响应也具有相似的振荡特性。设图1中, 木块质量m=1 kg, 阻尼系数f=0.2 N/ (m·s) , 弹性系数k=40 N/m, 则对应电路中, L=1 H, R=0.2 Ω, c=0.025 F。用电路仿真软件PSpice仿真, t=0 s时, 向电压源施加激励1 V时, 回路电流随时间的变化曲线中, 开始振幅很大, 随后沿中心线逐渐衰减, 最后回落至零。而用Matlab仿真图1机械系统时, 在施加力P为1 N时, 对应弹簧的振动速度dx/dt随时间的变化曲线也类似于上面电气系统进行模拟仿真的变化规律。这时, 如果参数选择适当, 其结果在数值上, 同利用机械系统进行试验的结果将完全相同, 因此可以将机械系统通过电气系统的相似模型来处理[9], 其相似对应关系如表1所示。

2.2 在生化的动力系统中的应用

在生化的动力系统分析中, 各种化学反应也是微观粒子在浓度差的驱使下发生变化的过程, 这也符合电网络理论体系中力-流的关系[10]。例如建立血管系统的电路模型, 为简化计算过程, 可作如下假设:血液为不可压缩的牛顿液体;血管为薄壁直圆管;除主动脉上升支和主动脉弓外, 其余部位血管中的流动均为轴对称的层流, 其流速呈抛物线分布;只考虑血压沿其轴向的变化;忽略血管壁的运动, 则可由Navier-Stokes方程得到各段血管压力及流量间的关系, 具体的关系式如下:

$\begin{array}{l} q_{n} = \frac{1}{L_{n}} \int (p_{n - 1} - p_{n} + p_{G n} - R_{n} q_{n}) d t (8) \\ p_{n} = \frac{1}{C_{n}} \int (q_{n} - q_{n + 1}) d t + \frac{0.002 (q_{n} - q_{n + 1})}{C_{n}} (9) \end{array}$

即每段血管可由四元件等效电路来模拟, 其中:pn, qn为第n段血管的压力和流量;pn-1为前一段血管的压力;qn+1为下一段血管的流量;式 (9) 中的最后一项为考虑血管的粘滞性而加入的修正项。

图3为所研究对应血管段的等效电网络图, PE表示各血管段所受外加压力;Rn, Ln及Cn分别表示第n段血管的等效流阻、等效流电感及血管等效电容, 其值可由式 (10) ～式 (12) 求出:

$\begin{array}{l} R_{n} = {\begin{cases} 8 π μ_{n}^{3} / V_{n}^{2}, V_{n} > V_{n 0} \\ 8 π μ_{n}^{3} V_{n 0} / V_{n}^{3}, V_{n} < V_{n 0} \end{cases} (10) \\ L_{n} = 2 ρ l_{n}^{2} / V_{n} (11) \\ C_{n} = {\begin{cases} 3 π r_{n}^{3} l_{n} / (2 E h_{n}), V_{n} > V_{n 0} \\ 60 π r_{n}^{3} l_{n} / (2 E h_{n}), V_{n} < V_{n 0} \end{cases} (12) \end{array}$

式中:ρ为血液密度;Rn为血管半径;ln为血管长度;hn为管壁厚度;μ是粘滞系数;E是管壁的弹性模量;Vn为血管容积;Vn0是血管内外跨壁压为0时该血管段的容积。

2.3 热力学的电网络模型

这里打算从一个热力系统开始, 以便构建每个能量节点的端口网络[11]。例如, 存在一个含有强度量温度T和广延量熵S的热力端口。假定流入端口的流为正, 热力学绝对温度的定义为:

$Τ \equiv (\partial U / \partial S) (13)$

设系统中其他所有广延量为常数, 则这个端口能量U的变化dU为:

$d U = (\partial U / \partial S) d S = Τ d S (14)$

把熵dS看作是流, 温度看成力, 这个端口就类似于在非平衡网络中的端口。同样, 各种形式功的端口也可以类似地表示为:

机械功:

$d U = Ρ d V (15)$

化学功:

$d U = A d ξ (16)$

物质输运功:

$d U = μ_{i} d Ν_{i} (17)$

式中:P, A, μi, dV, dξ, dNi分别表示压力、化学亲和力、组份i的化学势、体积变化、反应程度和组份i的摩尔数变化。

还可考虑一个含有两个自由度的均匀流体系统。该系统可以用一个二端口网络表示, 其中端口变量dS, dT为力, dP, dV为流。为便于对称处理, 将dV定义为系统外部的体积变化, 而不是内部体积变化。

该端口网络所表示的系统如图4所示, 对其可列出如下处理模型:

$\begin{array}{l} E_{1} = d Τ = (\partial Τ / \partial S)_{V} d S + (\partial Τ / \partial V)_{S} d V \\ = R_{11} F_{1} + R_{12} F_{2} (18) \\ E_{2} = d Ρ = (\partial Ρ / \partial S)_{V} d S + (\partial Ρ / \partial V)_{S} d V \\ = R_{21} F_{1} + R_{22} F_{2} (19) \end{array}$

对于非平衡态的系统而言, 当且仅当互易条件成立时, 该端口网络才是拓扑连通的, 因此对上述模型, 根据热力学麦克斯韦互易有:

$(\partial Τ / \partial V)_{S} = (\partial Ρ / \partial S)_{V} (20)$

对于构成n端口z参数电网络模型的情况, 同样可取T, P, S及V的任意2个端口增量为自变量[12], 根据麦克斯韦互易关系, 可以分别得到计算用矩阵U, H, A及G (见表2) 。

对于表2中的结论, 可通过运用电网络理论方法以理想气体为例进行分析, 以说明上述结果的合理性。根据理想气体的状态方程:

$Ρ V = n R Τ (21)$

可以得到G参数方程的系数:

$\begin{array}{l} g_{21} = (\frac{\partial V}{\partial Τ})_{Ρ} = \frac{n R}{Ρ} (22) \\ g_{22} = (\frac{\partial V}{\partial Ρ})_{Τ} = - \frac{n R Τ}{Ρ^{2}} (23) \\ g_{12} = - g_{21} (24) \end{array}$

根据微分原理, 设G=G (T, P) 为连续函数, 全微分的充分必要条件为:

$(\frac{\partial g_{11}}{\partial Ρ})_{Τ} = (\frac{\partial g_{12}}{\partial Τ})_{Ρ} = - (\frac{\partial g_{21}}{\partial Τ})_{Ρ} = (\frac{\partial^{2} V^{}}{\partial Τ^{2}})_{Ρ} = 0$ (25)

可看出, g11仅为P的函数, 即g11=g (p) , 因而G参数方程可以表示为:

$[\begin{array}{l} d s \\ d v \end{array}] = [\begin{matrix} g (p) & - \frac{n R}{Ρ} \\ \frac{n R}{Ρ} & - \frac{n R Τ}{Ρ^{2}} \end{matrix}] [\begin{array}{l} d Τ \\ d Ρ \end{array}] (26)$

由以上的分析可看出, 热力学的电网络形式在方法上找到了与热力系统参数一一对应的拓扑矩阵, 这为热力系统的电方法分析提供了严谨的理论基础, 其结果形式简单, 对为进一步开拓能量系统的电网络分析和应用提供了强大的理论及技术支撑。

3结语

阐述了以电路方法模拟动力系统的基本原则和法则, 提出主要以电阻、电容和电感基本元件为基础的模拟计算和分析, 以便于实现动力系统的电方法分析研究。在弹簧机械系统的分析中, 建立了相应的电网络拓扑结构, 实现了该系统的电网络建模, 并进行了仿真计算, 得出的结果说明了建模的合理性。在生态学及生物学的应用中, 围绕其动力学方面的问题, 应用电网络方法进行了分析处理, 通过对某段选定的血管拓扑结构建立实现了其电网络建模, 为进一步的生物学分析研究提供了相关的流体动力学基础。在工程热力学的基础理论方面, 提供了电网络方法和工程热力学在某些本质方面的一致性, 描述了以热力学参数为描述对象的端口网络和描述模型, 建立了工程热力学问题研究的电网络分析方法基础。

参考文献

[1]程伟良, 李艳秋, 夏国栋, 等.网络热力学——热力学发展的一个新逻辑阶段[J].西安交通大学学报, 2004, 38 (5) :548-550.

[2]卢虹冰, 白净, 张立藩.多元非线性人体循环呼吸系统模型及其应用[J].第四军医大学学报, 1999, 20 (3) :190-194.

[3]王向庭, 王芳林, 刘朋.基于PSpice的动力学微分方程的模拟解[J].机械, 2008, 35 (3) :18-20.

[4]陈熙, 李青, 李正宇, 等.热声热机网络模型演化与机理初探[J].低温工程, 2004 (1) :27-31.

[5]包家立, 杨媛媛, 王红.药物经皮转运通道的网络热力学模型[J].生物物理学报, 2004, 20 (5) :413-418.

[6]包家立, 王红, 葛霁光.药物经离体角质层渗透的非平衡态特性[J].药学学报, 2004, 39 (4) :296-300.

[7]包家立.脉冲电场诱致物质经角质层转运的唯象关系[D].杭州:浙江大学, 2003.

[8]陈建业.电力电子电路的计算机仿真[M].北京:清华大学出版社, 2003.

[9]KAY J J, GRAHARN L A, ULANOWICZ R E.A detailedguide to network thermodynamics analysis[J].Coastal andEstuarine Studies, 1989, 32:15-61.

[10]LEWIS Edwin R.Network thermodynamics revisited[J].Bio.Systems, 1995 (34) :47-63.

[11]MIKULECKY D C.Network thermodynamics and com-plexity:a transition to relational systems theory[J].Com-puters and Chemistry, 2001 (25) :369-391.

岩质边坡岩体动力破坏数值模拟分析篇10

一般认为, 岩质边坡岩体在地震作用下的破坏主要是动力引起的剪切破坏。然而在汶川地震后, 大量的滑坡表现上部拉张破坏, 下部剪切破坏。通过理论分析和数值模拟对这种现象进行了分析。

1 地震荷载

在边坡地震稳定性分析中, 通常是将地震作用简化成水平不变的加速度的作用, 以此加速度作用于不稳定体的质心的惯性力[2], 将地震力转化成应力的形式表示。常用的地震波都是以加速度的形式表示, 通过傅立叶转换得到速度时程, 然后根据波动理论公式 (1) 、 (2) [3], 动应力表达式如下:

σn=ρCnvn (1)

σs=ρCsvs (2)

式中:σn 为地震作用产生的垂直应力;undefined为岩土体中纵波传播的波速;ρ为岩土体密度;vn为岩土体中质点的竖向振动速度 (关于时间t的变量) ;σs为地震作用产生的水平应力;undefined为岩土体中横波的波速;vs为质点的水平振动速度 (关于时间t的变量) 。

汶川地震造成的边坡破坏现象表明, 岩体的破坏既有拉张破坏也有剪切破坏。岩体拉张破裂判据如下[4]:

1) 岩石是低抗拉介质。

2) 节理、断层、层理的抗拉强度更低。

3) 岩土体如果同时承受拉应力、剪应力, 首先发生拉张破坏, 然后才可能发生剪切破坏。

4) 当某一节点的最大拉应力大于岩石的抗拉强度时, 此节点就破裂, 即:

σ≥[σ] (3)

式中:σ为节点的最大拉应力。在动力作用下, 总应力为静力、动力同时作用下的应力值之和:σ=σz+σs。

式中:σz, σs分别为静力、动力作用下的应力值。

5) 当主压应力矢量方向与最大拉应力方向垂直时, 岩石破裂。

6) 在众多承受拉应力的应力场中, 承受最大拉应力的点优先破裂。

2 数值模拟分析

某均质岩质边坡, 模型尺寸如图1所示, 材料参数见表1。

输入的加速度地震波持时20s, 加速度峰值为6.4m/s2 (放大5倍) , 时步0.02s, 如图2所示。采用FLAC2D6.0的动力学模块进行数值模拟, 网格最大尺寸划分为2m。采用平面应变假设, 强度准则为摩尔—库仑强度准测。坡顶为自由表面;左侧和右侧边界施加水平方向的位移约束, 同时施加粘滞边界条件和自由场边界条件;模型底部设置水平和垂直方向位移约束, 施加粘滞边界条件, 地震波从模型底部施加。阻尼采用局部阻尼, 阻尼系数取0.157。

图3是岩质边坡分别在静力和动力作用下剪应力云图。从图3 (a) (静力作用下) 中可以看到, 剪应力从坡顶向坡脚逐渐增大, 处于同一量级, 无明显的剪应力集中区, 这与边坡中剪应力分布规律相符。从图3 (b) (动力作用下) 中可以看到, 坡顶附近的剪应力与静力作用下相比大一个量级, 坡腰向下延伸至坡脚部位的剪应力比静力作用下大两个量级, 在坡腰及坡脚附近形成了明显的剪应力集中区。由此可以发现, 在动力作用下, 边坡坡腰及其下部的岩体所受的剪应力明显大于坡腰及其上部岩体所承受的剪应力。同时, 图4分别是边坡在静力和动力作用下的水平向主应力云图。从图4 (a) 中可以看到, 在静力作用下, 边坡坡顶中没有出现拉应力 (正值为拉应力, 负值为压应力) ;而在动力荷载作用下, 如图4 (b) 所示, 在坡顶附近产生了拉应力 (图中的深蓝色部分) , 显然岩体这种低抗拉介质很容易在拉应力作用下发生拉破坏。综合对比分析边坡岩体在静力和动力作用下的应力云图可以看出, 岩质边坡在动力荷载的作用下, 坡顶附近岩体会产生相对较大的拉应力和相对较小的剪应力, 而坡腰向下一直延伸到坡脚的岩体中会形成明显的相对较高的剪应力。根据上述岩体抗拉判据可知, 岩体在同时承受拉应力、剪应力时, 首先发生拉张破坏。所以可以得出这样一个结论, 在地震动荷载的作用下, 坡顶附近岩体发生拉张破坏, 下部岩体主要发生剪切破坏。

3 结构面对岩体破坏影响的分析

岩体中充满了大量的软弱结构面, 其结构面是岩体中力学相对薄弱的部位, 它导致了岩体力学性能的不连续性、不均一性和各向异性。岩体中的软弱结构面通常是岩体稳定性的控制面。对于受结构面控制的边坡工程地质模型, 其变形破坏形式决定于结构面的形态及组合[5], 岩质边坡的稳定性主要由岩体结构面控制。因此, 对于岩质边坡而言, 其内部软弱结构面在动力作用下的强度特征及其动力响应将是决定边坡动力稳定性的重要因素。

根据孙进忠[6]的研究成果表明:在地震作用下, 结构面离坡顶越近, 结构面上部岩土体振动越强烈;结构面倾角从反倾到顺倾的变化过程中, 结构面上部岩土体的振动加速度越来越大。结构面的存在显著增强了结构面上部的动力响应, 并且造成了结构面附近动力响应急剧增大。结构面下部动力反应较无结构面时也略有增大, 幅度明显小于上部。由此可见, 岩体在地震作用下动力响应最强烈的地方是软弱结构面附近及其上部坡肩部位。边坡顶部动力响应幅值较之边坡底部存在明显的放大和集中现象, 尤其是在坡顶、坡肩部位。坡顶以下部位动力响应幅值减小较快。所以, 地震作用造成动应力在坡顶、坡肩部位相对较大, 这些部位在动应力和静应力的耦合作用下容易产生拉应力。所以对于岩质边坡, 坡顶和坡肩部位的破坏应该是以沿结构面、不连续面的拉张破坏为主。随着离坡顶距离的增大, 边坡动力响应逐渐减弱, 岩土体的自重应力增大, 岩土体在动应力和静应力的作用下产生的拉应力逐渐减小, 直至到不产生拉应力, 这时将发生剪切破坏。

4 结论

1) 结构面的存在显著增强了岩质边坡的动力响应, 并且造成了结构面附近动力响应急剧增大, 形成了动应力集中区。

2) 对于岩质边坡, 坡顶附近发生以结构面的拉张破坏为主的破坏, 下部发生剪切破坏。

摘要：汶川地震触发的大量岩质滑坡表明, 坡顶岩体发生拉裂破坏, 下部岩体发生剪切破坏。利用FLAC2D6.0的动力学模块对某岩质边坡在动荷载下的破坏进行了数值模拟分析, 结果表明, 对于岩质边坡、坡顶附近岩体发生的是以结构面的拉张破坏为主的破坏, 坡体下部发生剪切破坏。

关键词：岩质边坡,动荷载,拉破坏,剪切破坏,数值模拟

参考文献

[1]李忠生.国内外地震滑坡灾害研究综述[J].灾害学, 2003, 18 (4) :64-70.

[2]刘红帅, 薄景山, 刘德东.岩土边坡地震稳定性评价方法研究进展[J].防灾科技学院学报, 2007, 9 (3) :20-27.

[3]Itasca Consulting Group Inc.FLAC2D users manual DY-NAMIC ANALYSIS1-27[R].Minneapolis:Itasca Con-sulting Group Inc, 2008.

[4]王来贵, 赵娜, 李天斌.强震诱发单一弱面斜坡塌滑有限元模拟[J].岩石力学与工程学报, 2009, 28 (增1) :3163-3167.

[5]祁生文, 伍法权, 刘春玲, 等.地震边坡稳定性的工程地质分析[J].岩石力学与工程学报, 2004, 23 (16) :2792-2797.

本文来自 360文秘网(www.360wenmi.com)，转载请保留网址和出处