^`
UDEC 实例翻译与命令解析
翻译：珠穆朗玛
1 地震诱发地层坍塌 Seismic-Induced Groundfall
1.1 问题描述
本例展示使用 UDEC 模拟分析地震诱发地层坍塌的一类的问题,模型见图 1.1,该模型基 于加拿大安大略省萨德伯里市鹰桥公司弗雷则矿 34-1-554 切割断面的一个剖面图的结构和 尺寸. 用二维平面应变模型代表垂直于超采轴向方向的平面效应,超采面高 5m,宽 10m.
假定两个连续节理交叉平面分析:一个角度为 45 度,另一个为-9 度,两者节理间距均为 5m,为了演示的目的,一个近似垂直的“虚拟节理”也被添加到块体内开挖面顶部以增强不稳 定性。围岩参数来自试验室平均测试数值，假定岩石块体参数如下：
假定块体仅具有弹性行为，节理假定符合库伦滑动准则，选择典型的教课书数值作为节 理参数，如下：
初始应力状态按各向同性估计为24Mpa（假定垂直荷载由覆盖深度大约800m 的岩层产生）。
1.2 UDEC 分析
UDEC 模拟顺序分三个阶段,首先,模型在初始应力状态下进行无超采固结.其次,进行开挖并且模型循环至平衡状态.本阶段超采面周围的应力分布见图1.2.超采正上方和下方的块体滑动后稳定.在第三阶段.估计了两个不同的峰值速度的地震事件.对所有地震模拟,在问题域的外周边界引入粘滞边界用以消除波的反射.从而模拟有限的岩体,地震事件用施加到模型顶部y 方向的正弦应力波表现.应力波被叠加到已存在的初始地应力上.在第一个模拟中,施加1.25Mpa 的峰值应力,应当注意的是,由于粘滞边界条件实际是在模型顶部, 施加的有效影响应力应该是1.25 MPa/2, or 0.625 MPa.0.02 秒后的开挖面拱顶的应力分布见图1.3,两点的位移被监测,1 点位于开挖面的左角,点2 位于拱顶块体的右角, 图1.4 的位移时间曲线显示两点本质上是弹性反应
.
本例关心的问题是在模型顶部施加的速度和计算速度的对比,下面的公式可以用以估计施加的波速.
使用这个方程,施加的最大波速大概是0.04m/sec,图1.5 显示的峰值波速小于0.06m/sec. 估计的波速和监测波速的不同在于使用的围岩模量.而是没有考虑节理变形的相等变形模量.
在第二个案例中,施加应力波峰值12.5 Mpa(有效应力6.25Mpa).0.02 秒后的开挖拱顶应力分布见图1.6.该图显示出拱顶岩体不受力,表面该块体已经松散并正在下落.对于关心的问题,后来三个时间的几何体和应力分布见图1.8 至图1.10.
在问题的顶部预测的波速(从上面的方程)是0.4m/sec.从模型中计算的波速见图1.11,再次,由于使用的是原岩弹性模量而不是岩体的变形模量导致预测和监测的波速之间的差异.
1.3 节包含了该模型的数据列表,该列表包含了一个FISH 函数(show)被用来创建坍塌的动画文件,每隔0.02 秒俘获一个显示的图片.通过改变FISH 参数time_int 可以改变动画帧的
间隔.视图的总数也可以通过改变snap_shot 的数值进行改变.为了显示80 帧的显示图片而创建的该电影文件需要大概13MB 的硬盘空间.
1.3 数据文件列表
Example 1.1 SEISMIC.DAT
title
SEISMIC INDUCED ROOF COLLAPSE 地震诱发拱顶坍塌
;
round 0.01
; define original boundary of modeled region 定义模型区域的原始边界
block -25,-20 -25,20 25,20 25,-20
; generate joint pattern over entire original region 在整个原始区域生成节理形态
jregion id 1 -25,-25 -25,25 25,25 25,-25
jset 45,0 200,0 0,0 5.0,0 (0,0) range jreg 1
jset -9,0 200,0 0,0 5.0,0 (0,0) range jreg 1
; put in joints needed for the later excavation 为了后面开挖而设置的节理
crack -5.01,-2.51 5.01,-2.51
crack -5.01, 2.51 5.01, 2.51
crack -5,-2.5 -5,2.5
crack 5,-2.5 5,2.5
crack 2.25,2.5 1.93,5.0
; generate fdef zones and assign joint properties (mat=1 & jmat=1;default) 生成单元和设置节理参数generate edge 9.0 range -30,30 -30,30
prop mat=1 d=0.00300 k=39060 g=31780
prop jmat=1 jkn=20000 jks=20000
prop jmat=1 jf=30.0
; apply boundary conditions and initial conditions to 在地应力下施加边界条件和初始条件
; consolidate model under field stresses
bound stress=-24.0, 0.0, -24.0 ygrad=-0.3 0 -0.3
insitu stress=-24.0, 0.0, -24.0 ygrad=-0.3 0 -0.3
bound yvel 0.0 range -26,26 -21,-19
grav 0.0 -10.0
; track the x-displacement, and y-displacement over time 追踪位移
hist solve
hist xdis=0,7 ydis=0,7 type 1
solve rat 1e-5
; save consolidated state
save seismic1.sav
; make excavation
delete range -5,5 -2.5,2.5
solve rat 1e-5
; save excavated state
save seismic2.sav
;
rest seismic2.sav
; apply seismic load from top (peak velocity=0.04 m/sec)
;
; set up nonreflecting boundary
bound mat=1
bound xvisc range -26 -23 -21 21
bound xvisc range 23 26 -21 21
bound xvisc yvisc range -26 26 -21 -19
bound xvisc yvisc range -26 26 19 21
; apply sinusoidal stress wave
bound stress 0 0 -1.25 yhist=cos(100.0,0.0195) range -26 26 19 21 ;
reset time hist disp rot
hist ydis (-4.48,2.57)
hist ydis (0,2.57) yvel (0,2.57) yvel (4,2.57) yvel(-4.48,2.57)
hist yvel (0,20) yvel (25,10) yvel (25,-10) yvel (0,-20)
hist yvel (-25,-10) yvel (-25,10)
hist sxx (25,10) sxx (25,-10) sxx (-25,-10) sxx (-25,10)
hist syy (0,20)
;
damp 0.1 1.0 mass
; 0.02 sec.
cyc time 0.02
save seismic3.sav
;
rest seismic2.sav
; apply seismic load from top (peak velocity=0.4 m/sec)
; set up nonreflecting boundary
bound mat=1
bound xvisc range -26 -23 -21 21
bound xvisc range 23 26 -21 21
bound xvisc yvisc range -26 26 -21 -19
bound xvisc yvisc range -26 26 19 21
; apply sinusoidal stress wave
bound stress 0 0 -12.5 yhist=cos(100.0,0.0195) range -26 26 19 21 reset time hist disp
hist ydis (-4.48,2.57)
hist ydis (0,2.57) yvel (0,2.57) yvel (4,2.57) yvel(-4.48,2.57)
hist yvel (0,20) yvel (25,10) yvel (25,-10) yvel (0,-20)
hist yvel (-25,-10) yvel (-25,10)
hist sxx (25,10) sxx (25,-10) sxx (-25,-10) sxx (-25,10)
hist syy (0,20)
;
damp 0.1 1.0 mass
save seismov.sav
;
; 0.02 sec.
cyc time 0.02
save seismic4.sav
; 0.25 sec.
cyc time 0.23
save seismic5.sav
; 0.50 sec.
cyc time 0.25
save seismic6.sav
; 0.75 sec
cyc time 0.25
save seismic7.sav
;
rest seismov.sav
; make a movie of the groundfall
;
wind -12 12 -12 12
set ovtol 0.05
plot block vel max 2.0 blue stress max 50
movie on
movie file = seismic.dcx
movie step 1000
step 40000
3 隧道支护荷载Tunnel Support Loading
3.1 问题陈述
本例模拟展示了UDEC 在检查衬砌隧道方面的应用,着重强调了荷载在混凝土衬砌中的发展,本例也解释了模拟连续建造操作中独立阶段的模拟程序.
隧道系统的理想几何体见图3.1.系统包含在海床下大约70m(中线)深度,中线间距12m 的两个隧道, 初始水位在隧道中线上方110m 处.服务隧道直径5.24m,衬砌厚度37cm.主隧道直径8.22m,衬砌厚度46cm.服务隧道先于主隧道开挖和衬砌.随后设置主隧道衬砌,水位上升增加到100m.
施工顺序是:
(1)开挖服务隧道excavation of the service tunnel;
(2)衬砌服务隧道lining of the service tunnel; (3)开
挖主隧道excavation of the main tunnel; (4)衬砌主
隧道lining of the main tunnel; and (5)升高水位
raising of the water level.
分析的目的是评价每个施工阶段服务隧道和主隧道支护状况.
本例的材料参数见下:
岩体——开挖隧道的围岩参数为：
弹性模量elastic modulus 0.89 GPa
泊松比Poisson’s ratio 0.35
单轴抗压强度uniaxial compressive strength 3.5 MPa
粘聚力cohesion 1 MPa
密度density 1340 kg/m3
混凝土衬砌——弹性模量为24 GPa ，泊松比为0.19. 假定衬砌为线弹性材料。
3.2 UDEC 分析
本问题创建的UDEC 模型见图3.2，隧道中心线在Y=-70 的位置，注意，模型边界距隧
道开挖面很近，本模型为了演示的目的只是提供一个快速计算，实际计算时应采用大的模型。
底边和侧边界采用滚轴固定，海底之上海水的重量按30m 的水头施加在模型表面上，假定隧道为线性防水衬砌，因此不需要进行瞬态地下水流分析。忽略孔隙水压力，而将岩体单位重量设置为浮容重。垂直对水平应力比率假定为静水压力0.3。
对本例而言，上面列的5 个施工顺序按三个阶段进行模拟，主隧道开挖和服务隧道衬砌按一个瞬时的活动进行模拟。主隧道衬砌和水位上升也假定瞬时产生，这些活动可以通过线性（梯度）降低隧道周围反力或者设定产生一些松弛后安装支护而被分开模拟。
在模拟的第一阶段，在重力应力被初始化到块体中后，服务隧道被开挖，UDEC 循环直到达到平衡状态，弹性位移结果见图3.3.
在第二阶段, 衬砌服务隧道和开挖主隧道.采用16 个梁单元来模拟服务隧道混凝土衬砌, 图3.4 和图3.5 显示了主隧道开挖后的位移和主应力分布状况.注意,整个服务隧道向主隧道转移.
在第三个阶段,一个额外的荷载被施加在模型的顶部用来模拟增加的100m 的水的重量, 使用STRUCT apply pressure 命令来施加静水压力荷载(施加到隧道衬砌上),用以代表隧道中心线上210m 的水位.图3.6 和图3.7 显示了施加额外荷载后相应的轴力和弯矩分布,服务隧道最大轴力为5.59×106N,主隧道最大轴力为9.06×106N.
3.3 数据文件列表
; tunnel support loading
;
block 0 -30 60 -30 60 -90 0 -90 round 0.1
crack 0 -70 60 -70
crack 30 0 30 -90
crack 42 0 42 -90
tun 30 -70 4.11 12
tun 42 -70 2.62 8
tun 30 -70 5.5 12
tun 42 -70 5.5 12
gen edge 2.0Example 3.1 TUNNEL.DAT
海床面
0,-30
30,0
海水面
42,0
60,-30
; initial stress state
bound stress 1.02e5 0 1.02e5 ygrad 1.34e4 0 1.34e4 ;bound stress sxxo sxyo syyo ygrad sxxy sxyy syyy ; sxx = sxxo + (sxxx·x) + (sxxy·y)0,-70
0,-9030,-90 42,-90
60,-70
90
insit stres 1.02e5 0 1.02e5 szz 1.02e5 ygra 1.34e4 0 1.34e4 zgra 0 1.34e4 gravity 0 -10
; rock properties
prop mat=1 d=1340
zone model mohr
zone shear=.33e9 bulk=.99e9 coh=1e6 fric=30.0
;prop mat=1 d=1340 g=.33e9 k=.99e9 coh=1e6 fric=30.0
;change con=3
; elastic joint properties
prop jmat=1 jkn=1e9 jks=1e9 jcoh=1e10 jtens=1e10
; cycle to initial equilibrium
hist solve type 1
solve rat 1e-5
save tun1.sav
; excavate service tunnel
del 40 44 -72 -68
bound -1 1 -91 0 xvel=0.0
bound -1 90 -91 -89 yvel=0.0
bound 59 61 -91 0 xvel=0.0
; histories around tunnel 1
hist ydis 42 -67 sxx 42 -67
hist ydis 42 -73 sxx 42 -73
hist xdis 39 -70 syy 39 -70
hist xdis 45 -70 syy 39 -70
; histories around tunnel 2
hist ydis 30 -65.0 sxx 30 -65.0
hist ydis 30 -75.0 sxx 30 -75.0
hist xdis 25.0 -70 syy 25.0 -70
hist xdis 35.0 -70 syy 35.0 -70
reset disp jdisp
solve rat 1e-5
save tun2.sav
; line service tunnel
struct gen xc=42 yc=-70 npoint=16 mat=5 thick=0.37 fang=-11.25 theta 360 prop mat=5 st_d=2400 st_ymod=24.0e9 st_prat=0.20 st_yield=1e10
prop mat=5 if_kn=1e8 if_ks=1e7 if_coh=1.0e10
; excavate main tunnel
del 28 32 -72 -68
reset disp
solve rat 1e-5
save tun3.sav
; line main tunnel
struc gen xc=30 yc=-70 npoint=8 mat=5 thick=0.46 fang=22.5 theta -360
; add additional load representing raised water level
bound stress 0.0 0.0 -1.0e6 range -1 91 -31 -29
; add hydrostatic loads to tunnel liners
struct apply press 0.0 2.06e6
reset disp
solve rat 1e-5
save tun4.sav
ret
4 流体和动态荷载下的重力坝分析
4.1 问题描述
本例展示了一个坐落在节理岩体基础上的100m 高混凝土重力坝建造问题,节理平均间距是50m,节理方向为20 度和-70 度, 主要研究了两个荷载条件,首先是研究包含流体在岩体裂隙内渗流的水库蓄水效应分析.其次,为了研究潜在地震类型的荷载,一个动态波被施加在模型的基础上进行分析.
4.2 UDEC 分析
本问题的UDEC 模型见图4.1。模型理想化的展示了一个坐落在节理岩体基础上的重力坝，目的是展示这种类型问题的推荐求解过程。数据文件见4.3 节，按下面的顺序进行分析。
第一步：重力荷载——空的水库
假定初始应力状态下的有效应力比率为σH/σV = 0.69，水位假定在y = 0，由大坝重量产生的初始应力状态和未蓄水的水库见图4.1。注意，使用INSITU 命令设定的应力是总应力，被加载在块体上。对节理裂隙，UDEC 计算有效应力并将域压力设置到静水压力。
第二步：水库蓄水本阶段，假定水位上升到大坝的顶部，在大坝上游侧和岩石基础处施
加静水压力作用，
在侧边界由于施加在大坝上的荷载产生的水平反作用力被认为是滚轴支撑。
对于流体渗流假定下面的条件：
1. 假定沿着模型底部和侧边的节理裂隙连接渗透性为零。
2. 大坝上游的岩体面，通过使用BOUND pp 命令将水头固定在100m（0.98Mpa）。下游
面，水头设置为零。
3. 大坝和岩体基础交界面假定有较低的渗透性。
4. 使用稳定流体算法（SET flow steady）
选定的第二步结果见图4.2 至图4.6，蓄水后的位移结果见图4.2。大坝顶部的X 方向和Y 方向位移历史见图4.3。后来的图形显示蓄水后模型处于平衡状态。图4.4 的流体速率显示了大多数流体被集中在大坝基础下面的裂隙内。图4.5 显示了沿着节理（at x = -33.42, y = -30.37）的法向和剪切位移。正法向位移表明在本阶段节理张开。在位置x = -22.1, y = -26.3 处沿着裂隙的流体压力历史见图4.6。
第一部分：数值模拟技术研究文献综述 浅析数值模拟技术 1．引言 近年来，随着我国大规模地进行“西部大开发”和“南水北调”等巨型工程，越来越多的岩土工程难题摆在我们面前，单纯依靠经验、解析法显然已不能有效指导工程问题的解决，迫切需要更强有力的分析手段来进行这些问题的研究和分析。自R.W. Clough 上世纪60年代末首次将有限元引入某土石坝的稳定性分析以来，数值模拟技术在岩土工程领域取得了巨大的进步，并成功解决了许多重大工程问题。特别是个人电脑的普及及计算性能的不断提高，使得分析人员在室内进行岩土工程数值模拟成为可能。在这样的背景下，数值模拟特别是三维数值模拟技术逐渐成为当前中国岩土工程研究和设计的主流方法之一，也使得岩土工程数值模拟技术成为当今高校和科研院所岩土工程专业学生学习的一个热点。 采用大型通用软件对岩土工程进行数值模拟计算，在目前已成为项目科研、工程设计、风险评估等岩土类项目的必须，学习和掌握Ansys、FLAC3D、UDEC 等数值计算软件已成为学校、科研院所对工程从业人员的基本要求。 数值模拟方法主要有限元法、边界元法、加权余量法、半解析元法、刚体元法、非连续变形分析法、离散元法、无界元法和流形元法等，各种方法都有其对应的软件。 2.数值模拟的发展趋势 可以说, 继理论分析和科学试验之后, 数值模拟已成为科学技术发展的主要手段之一。随着软件技术和计算机技术的发展, 目前国际上数值模拟软件发展呈现出以下一些趋势： （1）. 由二维扩展为三维。早期计算机的能力十分有限，受计算费用和计算机储存能力的限制，数值模拟程序大多是一维或二维的，只能计算垂直碰撞或球形爆炸等特定问题。随着第三代、第四代计算机的出现, 才开始研制和发展更多的三维计算程序。现在，计算程序一般都由二维扩展到了三维，如LS-DYNA2D 和LS - DYNA3D、AUTODYN2D 和AUTO-DYN3D。 （2）.从单纯的结构力学计算发展到求解许多物理场问题。数值模拟分析方法最早是从结构化矩阵分析发展而来，逐步推广到板、壳和实体等连续体固体力学分析，实践证明这是一种非常有效的数值模拟方法。近年来数值模拟方法已发展到流体力学、温度场、电传导、磁场、渗流等求解计算，最近又发展到求解几个交叉学科的问题。例如内爆炸时，空气冲击波使墙、板、柱产生变形，而墙、板、柱的变形又反过来影响到空气冲击波的传播，这就需要用固体力学和流体动力学的数值模拟结果交叉迭代求解。 （3）.由求解线性问题进展到分析非线性问题。随着科学技术的发展，线性理论已经远远不能满足设计的要求。诸如岩石、土壤、混凝土等，仅靠线性计算理论就不足以解决遇到的问题，只有采用非线性数值算法才能解决。众所周知，非线性的数值模拟是很复杂的，它涉及到很多专门的数学问题和运算技巧，很难为一般工程技术人员所掌握。为此，近年来国外一些公司花费了大量的人力和资金，开发了诸如LS- DYNA3D、ABAQUS和AU-TODYN等专长求解非线性问题的有限元分析软件，并广泛应用于工程实践。这些软件的共同特点是具有高效
Project #1 硅的晶格常数和体弹模量的计算 一、平衡晶格常数和内聚能 自然条件下硅为金刚石结构(dc )。计算模拟时，我们可以假定它为各种结构，f cc, bcc, sc, dc. 可以预测，模拟的dc 结构的硅的体系能量最低，也即最稳定。下面我们将运用LAMMPS 来对硅的各种结构进行模拟。 定义晶格能量为Φ， 数密度为 ρ： pot E N Φ= N V ρ= 其中E pot 为势能， N 为体系总原子数，V 为体系的体积。选取 Stillinger-Weber (SW)，以下面命令执行 lammps 运算： 其中，lmp_serial 为 lammps 命令；”<” 符号为读取符；in.Silicon 为输入文件，里面包含运算所需要的各种数据和命令；-log 指定输出文件的名称。 可以看到屏幕上显示出lammps 运行的信息。这个计算量很小，所以很快就结束。接下来以如下命令来查看计算得到的数据： grep 是linux 中一个很重要的命令，用来搜索文本，读取匹配的行并打印出来。这里是搜索 dc.log 文件，将 @ 开头的行打印出来。如下： 晶格参数为5.4305埃,数密度为0.0499540303，每个原子的能量为-4.336599609eV.
下面具体来看刚才给的输入文件，in.Silicon . dc.log 文件中有原子总数的信息， 每个金刚石晶胞中有8个原子，383216?=，所以是216个原子。如下给出各种结构下的体系的原子数：
晶体结构类型 晶胞中的原子数 总原子数 简单立方SC 1 27 体心立方BCC 2 54 面心立方FCC 4 108 金刚石DC 8 216 表1. 不同晶体结构中的原子数 下图是计算模拟得出的各种结构下的数密度与每个原子能量的关系图。 横坐标为数密度， 以金刚石为例，ρ= 8/5.4315^3=0.049926，也即我们直接通过 grep 命令得到的第二项值；纵坐标为每个原子的能量，为第三项值。 金刚石之外，还需计算其他结构。只需对 in.Silicon 做稍微改动： 首先，将in.Silicon 复制成in.fcc : 然后编辑 in. fcc 改动如下几项： 然后如下命令执行： 相应的，如下命令查看log 文件中的数据：
华中科技大学体育馆数值模拟分析 6.1分析模型的建立 采用有限元软件ANSYS建立该网壳结构有限元分析模型。整体屋盖结构共计1481个节点，4430个单元，16种截面类型。建模时，网壳结构主体结构部分 （包括主桁架、次桁架、水平支撑和檩条）采用ANSYS的LINK8杆单元建模，两侧翼的主梁、次梁和支承钢管柱均采用BEAM4梁单元，网壳结构屋面下部混凝土支承结构亦采用BEAM4梁单元。分析时，屋面板、设备管线等荷载等效为节点荷载，施加在结构节点上。 在网壳结构有限元分析中，对于杆件采用的LINK8 3-D Spar单元为三维单元，假设材料为均质等直杆，且在轴向上施加载荷，可以承受单向的拉伸或者压缩，每个节点上具有三个自由度，即沿X、丫和Z坐标轴方向。该单元具有塑性、蠕变、应力硬化和大变形等功能，能较好的模拟三维空间桁架单元。 对于两侧翼结构和下部支撑体系的柱、梁等结构采用的BEAM4单元是一个轴向拉压、扭转和弯曲单元，每个节点有三个平动自由度和三个转动自由度，具有应力刚化和大变形功能。 施工过程模拟分析时考虑时，同时考虑温度效应影响，计算时材料假定为理想弹塑性材料。
6.2分析工况选取 按照实际施工顺序，将网壳结构屋盖施工过程划分为 5个工况进行施工数值 模拟，计算温度取为该阶段施工完成时的环境温度。 工况1: 7榀拱形主桁架安装完毕，但临时支撑未撤除，计算温度为温度15C ； (a )短轴立面 1' 裁灌 1::“ 麻 MM ■■号 -? ■?■ 水 * S-ES-B- Kir *￡靈曲 r-…；： 1 ;u * 图6-1有限元分析模型
图6-4工况2中屋盖结构平面图图6-5工况2中屋盖结构立面图工况3：次桁架、水平支撑及楼梯安装完毕，临时支撑拆除，计算温度为29C; 图6-6工况2中屋盖结构平面图 (b )长轴立面 图6-7工况2中屋盖结构立面图工况4：檩条及设备管线安装完毕，计算温度为41C; 完成后拆除其临时支撑，计算温度为8C； (a)短轴立面 (b)长轴立面 (a)短轴立面
Project #2 金属中的点缺陷：空位和间隙原子 一、空位 从晶体中移去一个原子，即可形成空位。本例将运用 LAMMPS 计算空位形成能, E v. LAMMPS 输入文件为in.vacancy 1) 在 fcc 结构的完整Cu晶体中引入一个空位 沿<100>方向构造一个 4 ×N×N×N 的晶体。N为input 文件中lattice命令指定的个方向上的晶胞重复单元数。 2) 弛豫 当一个原子从晶体中移走之后，周围的原子将相应地调整位置以降低体系势能。为得到稳定的构型，需要对体系进行弛豫，relaxation. LAMMPS提供两种能量最小化方式，cg 和 sd。本例中选用 sd 方式进行能量最小化。 如下是输入文件，in.vacancy:
3) 运行lammps 4) 计算空位形成能 空位浓度由下式给出： [n ] = exp( ? F v / k B T ). 其中 F v = E v ? TS v 为形成一个空位所需要的Helmholtz 自由能. 忽略熵S v , 空位浓度公式简化为 [n ] = exp( ?E v / k B T ). 设 E 1 为完整晶体能量，含N 个原子；E 2 为弛豫后的晶体能量，含N – 1个原子。空位形成能 E v 为： 211v N E E E N -≡- 或 ()21v coh E E N E ≡--, 其中 E coh = E 1 / N , 为完整晶体的内聚能。 本例中以EAM 模型计算4×(20×20×20)=32000个原子的体系，得到空位形成能E v ~ 1.26 eV ，文献中的实验值为~1.28 eV ，符合较好。 另由上式计算得到，300K 温度下的空位浓度为~ 7.59×10-22 ，1350 K (T m ) 时的空位浓度~ 2.2×10-5（文献中的实验值为~2×10-4 ）。换算时注意(1 eV/k B = 1.1604×10?4 K) 图1. 空位处于4×(6×6×6) 的 FCC 晶体中心，106c a =,206c a =,306c a =. 颜色依据原子势能标注。
Project #4 表面与界面能 铜的表面能 当物体形成表面时，表面上的原子键发生断裂，接近表面的几层原子不再如之前处于平衡状态，从而导致能量的升高，升高的温度便是物体的表面能。 利用LAMMPS 做出 20*20*40 fcc 的盒子，删去边缘的原子制造出一段真空层；算出此时体系的总能量0E ，然后从中间把盒子切成两半并移至足够远的距离，此时的体系总能量为E final ， 从而表面能： 02final surface E E A γ=? A 为表面的面积 (100) 面与 (111) 面 如下是输入文件in.surface_Cu_100 # LAMMPS Cu _Surface_100 units metal boundary p p p atom_style atomic lattice fcc 3.61 region box block 0 20 0 20 0 40 create_box 1 box create_atoms 1 box timestep 0.005 thermo 5 pair_style eam/alloy pair_coeff * * jin_copper_lammps.setfl Cu region boundary1 block INF INF INF INF 29.9 INF region boundary2 block INF INF INF INF INF 9.9 group boundary1 region boundary1 group boundary2 region boundary2 group boundary union boundary1 boundary2
数值模拟步骤
数值模拟 1、CFD方法简介 利用CFD方法，采用流体力学分析软件Fluent对三相分离器的流场进行了研究和分析，为实验研究提供理论支持。 CFD是英文Computational Fluid Dynamics(计算流体动力学) 的缩写,是一门用数值计算方法求解流动主控方程以发现各种流动现象规律的学科]。用CFD 技术进行数值求解的基本思想是: 把原来在空间与时间坐标中连续的物理量的场, 用一系列有限个离散点上的值的集合来代替, 通过一定的原则来建立离散点上变量值之间关系的代数方程, 求解代数方程以获得所求解变量的近似值。其主要用途是对流态进行数值仿真模拟计算，因此，CFD技术的用途十分广泛，可用于传质、传热、动量传递及燃烧等方面的研究。 流体机械的研究中多用CFD方法对分离器进行仿真模拟，其基本应用步骤如下： 1) 利用Gimbit进行前处理 a. 根据分离的形状、结构及尺寸建立几何模型； b. 对所建立的几何模型进行网格划分； 2) 利用Fluent进行求解 a. 确定计算模型及材料属性； b. 对研究模型设置边界条件； c. 对前期设置进行初始化，选择监视器，进行迭代计算； 3)利用Fluent进行后续处理，实现计算结果可视化及动画处理。
上述迭代求解后的结果是离散后的各网格节点上的数值,这样的结果不直观。因此需要将求解结果的速度场、温度场或浓度场等用计算机表示出来，这也是CFD 技术应用的必要组成部分。 利用CFD方法进行仿真模拟可以对分离器的结构设计及参数选择作出指导，保证设计的准确度，也可以为分离器样机的试验提供理论参考。由于CFD仿真模拟的广泛使用及其重要性，国内外很多学者，如Mark D Turrell、M.Narasimha、师奇威等都对其进行了研究，尤其是A.F. Nowakowski及Daniel J.SUASNABAR等人]对CFD 技术在旋流器模拟方面的应用做了详细的介绍，这些工作对CFD技术的发展起到了积极的促进作用。 2、控制方程 流体流动要受物理定律的支配，基本的守恒定律包括：质量守恒定律，动量守恒定律、能量守恒定律。如果流动处于湍流状态，系统还应要遵守附加的湍流输运方程。 1、基本假设 （1）鉴于我国各主力油田采出液含水已达到80％以上，故以水代替采出液进行分析计算；天然气的主要成分是是甲烷，故采用甲烷的替代天然气的性质； （2）在工作状态下，流动不随时间变化，流动为稳态； （3）水在管内的流动可以简化成二维流动； （4）不考虑温度的影响，服从绝热流动基本方程。 2、基本控制方程 （1）连续性方程 在直角坐标系下的质量守恒方程又称连续性方程：
第40卷第2期 当 代 化 工 Vol.40，No. 2 2011年2月 Contemporary Chemical Industry February，2011 收稿日期： 2010-08-17 作者简介： 魏显达（1983-），男，硕士，黑龙江北安人，2007年毕业于大庆石油学院电子信息工程，研究方向：塔顶流出系统的腐蚀与防 基于 Fluent 的三通管数值模拟及分析 魏显达，王为民， 徐建普 （辽宁石油化工大学石油天然气工程学院， 辽宁 抚顺 113001） 摘 要：Fluent 软件作为流体力学中通用性较强的一种商业CFD 软件应用范围很广。通过利用Fluent 计算流体动力学（CFD）的软件，对石油工业系统中常见的三通管内部流体进行了模拟分析，得到了三通管内在流体流动时的速度、压力和温度场分布图，为石油管道中的流体输送提供了理论依据。 关 键 词：Fluent；三通管；模拟分析；分布图 中图分类号： TQ 018 文献标识码： A 文章编号： 1671-0460（2011）02-0165-03 Numerical Simulation and Analysis of Fluid in Three-way Connection Pipe Based on Fluent Software WEI Xian-da ，WANG Wei-min ，XU Jian-pu (Institute of Petroleum and gas engineering , Liaoning Shihua University, Liaoning Fushun 113001，China ） Abstract : As a commercial CFD software with good universality, the Fluent software has been used extensively. In this paper, Simulation analysis on fluid in the three-way connection pipe of the oil industry was carried out by the software of fluid mechanics computation .Then distribution graphs of velocity , pressure and temperature of fluid in the three-way pipe were gained ，which can offer theoretical basis on fluid transportation in the petroleum pipeline. Key words : Fluent three-way ；Connection pipe ；Simulation analysis ；Distribution graphs Fluent 是目前国际上比较流行的商用CFD 软件包，在美国的市场占有率为60%，广泛应用于流体、热传热和各种化学反应等有关工业。软件包括前处理器（利用Gambit 进行物理建模、网格划分和划定边界层条件）、求解器（根据专业条件不同，采用不同的求解器，并规定物性、外部工作环境和进行数值迭代）和后处理器（把一些数据可视化，满足用户的特定要求）。 三通管在石油工业中应用广泛，采用传统的设计开发方法，存在经济成本高，研发周期长等缺陷，耗费大量的人力、物力 [1-2] 。应用CFD 软件，能够在 相对较短的设计周期内，较低的成本运行下，准确模拟流动具体过程，如速度场、压力场和温度场等的时变特性等。CFD 技术已经成为不可缺少的设计手段。 本文利用Fluent 的超强数值计算和分析能力对三通管道内原油流动时的速度、压强和温度场进行了数值模拟和分析，为石油管道中的流体输送提供了可靠的理论依据。 1 数学模型的建立和分析 输油管道管中，原油在三通管内的流动属于湍流，简化方程管道内的流体流动满足质量守恒、动量守恒、能量守恒、状态方程等。 连续性方程（连续性方程式质量守恒定律在流体力学中的表现形式）在直角坐标系下表示为（（1）方程） [3-5] ： 0)()()(=??+??+??+??z y x t z y x νννρρρρ （1） 式中：V x ，V y ，V z 是速度矢量ν在x 、y 和z 轴方向的分量，t 是时间，ρ是密度。 最常用的湍流求解模型是标准k -ε湍流模型。它需要求解湍动能k ((2)方程)和耗散率ε((3)方程)，具体如下所示： Y G G x x M b k i t i k t k ?+++??+??=ρεσμρ μ)[(d d （2） K K k t C G C G C x x b K i t i εμρεσμερεεε2 231)(])[(d d ?++??+??= （3）
LAMMPS手册-中文解析一、简介 本部分大至介绍了LAMMPS的一些功能和缺陷。 1．什么是LAMMPS? LAMMPS是一个经典的分子动力学代码，他可以模拟液体中的粒子，固体和汽体的系综。他可以采用不同的力场和边界条件来模拟全原子，聚合物，生物，金属，粒状和粗料化体系。LAMMPS可以计算的体系小至几个粒子，大到上百万甚至是上亿个粒子。 LAMMPS可以在单个处理器的台式机和笔记本本上运行且有较高的计算效率，但是它是专门为并行计算机设计的。他可以在任何一个按装了C++编译器和MPI 的平台上运算，这其中当然包括分布式和共享式并行机和Beowulf型的集群机。LAMMPS是一可以修改和扩展的计算程序，比如，可以加上一些新的力场，原子模型，边界条件和诊断功能等。 通常意义上来讲，LAMMPS是根据不同的边界条件和初始条件对通过短程和长程力相互作用的分子，原子和宏观粒子集合对它们的牛顿运动方程进行积分。高效率计算的LAMMPS通过采用相邻清单来跟踪他们邻近的粒子。这些清单是根据粒子间的短程互拆力的大小进行优化过的，目的是防止局部粒子密度过高。在并行机上，LAMMPS采用的是空间分解技术来分配模拟的区域，把整个模拟空间分成较小的三维小空间，其中每一个小空间可以分配在一个处理器上。各个处理器之间相互通信并且存储每一个小空间边界上的”ghost”原子的信息。LAMMPS(并行情况)在模拟3维矩行盒子并且具有近均一密度的体系时效率最高。 2．LAMMPS的功能 总体功能： 可以串行和并行计算 分布式MPI策略 模拟空间的分解并行机制 开源 高移植性C++语言编写 MPI和单处理器串行FFT的可选性（自定义） 可以方便的为之扩展上新特征和功能 只需一个输入脚本就可运行 有定义和使用变量和方程完备语法规则 在运行过程中循环的控制都有严格的规则 只要一个输入脚本试就可以同时实现一个或多个模拟任务 粒子和模拟的类型： （atom style命令）原子粗粒化粒子DNA 全原子聚合物，有机分子，蛋白质，联合原子聚合物或有机分子金属粒子材料粗粒化介观模型延伸球形与椭圆形粒子点偶极粒子刚性粒子所有上面的杂化类型力场：）（命令：pair style, bond style, angle style, dihedral style, improper style, kspace style, tabulated.
一、简介 1.SiC热分解制备石墨烯 自2004年Novoselov、Geim和合作者们从石墨上剥离出世界上第一种二维材料——单层石墨：石墨烯（Graphene）以来，石墨烯就受到了科技界的广泛重视[1]。Novoselov 和Geim两人因此在2010年获得了诺贝尔物理学奖。因为石墨烯的独特特性，在许多技术领域例如光电子学上它都被寄予厚望。研究石墨烯这种材料相关的物理化学特性和发展大面积、高质量生长石墨烯的技术，同时将其与器件物理学联系起来是我们研究和应用石墨烯的必由途径。 石墨烯是由碳元素组成的二维六边形材料，其在光学、电学、热学、力学等性质十分优异。它有可能在后摩尔定律时代成为硅（Silicon）的继任者，在单分子气体传感器[2]、自旋电子学[3]、量子计算[4]、太赫兹振荡器[5]等等领域发挥重要作用。如今，从石墨上剥离出石墨烯仍然是一种重要的石墨烯制备方方法。然而，这种方法产生的石墨烯大小通常不超过1000 μm2，只适合实验室研究，尚不能在工业上大规模应用。科学家发展了其他的石墨烯制备方法，包括将石墨烯视作一种薄膜来生长的化学气相沉积（Chemical Vapor Deposition, CVD）法、热分解碳化硅法（SiC thermal decomposition）、氧化石墨烯还原法（Graphene oxide reduction）等。 CVD法通过使含碳气源在有催化作用的金属表面分解或者使溶入到这些有催化作用的金属中的碳（C）发生表面偏析，使得在金属表面生成石墨烯或者多层石墨烯（Few-Layer Graphene, FLG）。能否直接在半导体/绝缘体上生长石墨烯呢？碳化硅热分解成功的解决了这一问题。最早试图使六方晶系的SiC晶体石墨化的研究报告见于1961年，Badami在高温和真空环境下得到了发生了一定石墨化的SiC[6]。在一定的退火条件下，SiC晶体表面发生热分解，Si原子发生解吸附，而C原子留下来重新排列和组合可以生长成外延型的石墨烯层[7]。更细致的研究发现用热退火的方法在六方SiC的Si面上生长的石墨烯比C面有更好的可控性，例如：可以更好的控制石墨烯的层数。Si面上生长的石墨烯生长方向与基底晶体结构有密切关系，这样提供了在基底上均匀覆盖和特定方向生长石墨烯的可能性。特别地，石墨烯直接生长在半导体SiC上使得我们无
CFD数值模拟—某厂房车间案例 ——索斯系统送风设计的经典案例 以往对中央空调系统内气流组织的理解和室内空间速度场、温度场的设计，只停留在经验值和独立单元射流计算的高度，而现在有了CFD（Computational fluid Dynamics）－计算流体力学软件的简称，是专门用来进行流场分析、流场计算、流场预测的软件。 杜肯索斯与知名高校合作，建立了大空间气流组织实验室，并使用当前流行的CFD计算程序（GAMBIT EXCEED FLUENT），将先前的经验与现代计算机强大的迭代能力相结合（采用1.83GHZ双核处理器的计算机进行模拟，迭代次数3000－5000次，数据收敛一般需要数十小时）。通过此软件的应用，可以显示并分析封闭空间中的流场；在较短的时间内能预测效果，并通过改变各种参数，使送风效果设计达到最佳。同时将此设计在大空间气流组织实验室中模拟验证。CFD的数值模拟能使我们更加深刻地理解问题产生的机理，为实验提供指导，节省实验所需的人力、物力和时间，并对实验结果的整理和规律的得出起到很好的指导作用。 在此，我们对某车间厂房内空气流场的CFD数值模拟应用作一案例介绍，具体如下： ●案例环境 该工程车间厂房生产精密器件，对工作区域的风速、温度有相当高的要求，同时还要求工作区外测能形成风幕，有效隔离工作台面两侧设备产生的气味和油污。 ●CFD技术应用: 传统的点式送风模式是根本不能实现这种在整个长达几十米的矩形工作台面上形成均匀的速度、温度场，并形成如此大的风幕。这时，我们选择采用具有线式及立体面送风特点的纤维织物空气分布系统，同时采用CFD技术对该设计进行指导。该设计关键在于： 一、纤维材料渗透率的设计 二、空气分布系统射流孔孔径和孔间距的设计 三、整个系统风量、压力的匹配性设计 对该厂房采用CFD模拟主要目的是验证在初步设计的基础上，密闭空间的气流组织情况，速度场、温度场、密度场、压力场的分布情况。该项目设计由durkeesox技术中心完成，分别在5：30、6：30、8：30、3：30四个方向开0.5”孔（详细参照CAD
一年一度的毕业设计就要到来了，CAE软件依然是流体专业众多学子毕设的拦路虎，为了使各高校流体同学顺利完成毕业设计以及有志于在流体行业有一番作为的青年才俊迅速掌握一门技能，从而更好地适应职场需求，北京经纬云图仿真科技有限公司感谢各位同行的支持和厚爱，特组织各方面CAE软件专家长期进行免费在线网络培训，诚邀您的参与！特此声明：本公司的所有培训将主要以工程实例为基础进行，让您真正的学到知识，懂得原理，而不仅仅是简单的软件操作。最重要的一点是为了回馈广大同仁们，本公司的培训完全免费，机不可失失不再来啊！ 主办单位：北京经纬云图仿真科技有限公司 培训时间：每周五晚8：00---9：30 培训方式：在线免费网络培训 培训2000人QQ群：281194860（参与培训请加入，注明：FLUENT培训） 培训内容：见附录 附录：1，基础流动计算 以AICD装置为例，说明利用fluent进行基础的流动模拟的步骤，包括计算设置和简单的后处理 2，两相流Mix模型应用 以某烟雾报警器为例，利用mix两相流模型，预测烟雾报警器内部的烟雾浓度分布 3，空化问题 以某拉瓦尔喷管型的空化装置为例，利用mix两相流模型和空化模型，计算装置内的水的空化情况 4，水的蒸发和凝结 以某水蒸气动力装置为例，利用mix两相流模型和蒸发凝结模型（fluent自带蒸发凝结模型和udf编程），并利用瞬态计算的方法，得到装置的瞬时流场分布和水凝结情况 5，湿空气的露点问题 以某实验装置为例，说明露点问题与蒸发凝结问题的本质区别，利用mix两相流模型和udf 露点模型，计算过冷空气中的水析出的问题，并与理论结果进行对比 6，萃取问题（溶液间传质问题） 以某反应塔的一层为例，利用mix两相流模型、组分输运模型和udf传质模型，计算甲苯的萃取过程，以及利用瞬态计算的方法，得到脉动进口条件对于萃取过程的影响 7，水中气泡上升 以某鱼缸中通过气泡的模型为例，利用vof两相流模型和udf气泡源，利用瞬态计算的方法，得到水中的连续气泡上升的流动动画 8，融化、蒸发联合作用问题 以低压力环境中的固态铝加热为例，利用vof两相流模型、融化凝固模型、蒸发凝结模型和udf加热源，利用瞬态计算的方法，计算铝融化后的液面形状 9，粒子冲蚀问题
例一 units metal # 单位为lammps 中的metel 类型 boundary p p p # 周期性边界条件 atom_style atomic # 原子模式 lattice fcc 3.61 # Cu 的晶格常数3.61 region box block 0 4 0 4 0 4 # x,y,z 各方向上的晶胞重复单元数，也即区域大小 create_box 1 box # 将上述区域指定为模拟的盒子 create_atoms 1 box # 将原子按晶格填满盒子 pair_style eam # 选取Cu 的EAM 势作为模型 pair_coeff * * Cu_u3.eam # EAM 势文件名称 run 0 # 运行0 步，仅为启动lammps 的热力学数据计算 variable E equal pe # 定义变量E 为系统总势能 variable N equal atoms # 定义变量N 为系统总原子数 print "the number of atoms & system energy now are $N $E" # 打印信息 create_atoms 1 single 2.45 2.05 2.05 # 在该位置插入一个原子 min_style sd # 能量最小化模式，sd minimize 1.0e-12 1.0e-12 1000 1000 # 能量最小化参数，指数越大最小化程度越深 print "interstitial introduced, minimized: $N atoms, energy is $E" fix 1 all nvt 100 100 100 drag 0.2 # nvt 系综，原子数、体积和温度保持不变;T=100K timestep 0.005 # 步长0.005fs run 1000 # 运行1000 步 print "nvt performed, temperature up: $N atoms, total energy is $E" fix 1 all nvt 100 0.0001 100 drag 0.2 # nvt 系综，温度由100K 到0.0001K run 1000 # 运行1000 步 print "nvt performed, temperature down: $N atoms, total energy is $E" compute 3 all pe/atom # 计算每个原子的势能 compute 4 all ke/atom # 计算每个原子的动能 compute 5 all coord/atom 3.0 # 计算每个原子的近邻原子数 dump 1 all custom 1 dump.atom id xs ys zs c_3 c_4 c_5 # 将信息写入dump.atom min_style sd minimize 1.0e-12 1.0e-12 10000 10000 # 再次能量最小化 print "the final state: $N atoms, total energy is $E" # 打印信息 ~ 例二 boundary p s s # 边界条件，拉伸方向是周期性，其余是自由边界；如果是薄膜拉伸则是两个周期性，块体则是三个周期性 units metal #单位制定义为metal atom_style atomic #原子类型自动 neighbor 2.0 bin #截断半径相关的东西 neigh_modifydelay 1 check yes # 邻近原子列表更新速度
# Big colloid particles and small LJ particles units lj This command sets the style of units used for a simulation. It determines the units of all quantities specified in the input script and data file, as well as quantities output to the screen, log file, and dump files. Typically, this command is used at the very beginning of an input script.这个指令设定模拟的格式，它决定了在输入脚本文件和数据文件以及在屏幕上显示的输出物，日志文件和垃圾文件中所有的单元格式。典型的是，它经常用于输入脚本的开头 For style lj, all quantities are unitless. Without loss of generality, LAMMPS sets the fundamental quantities mass, sigma, epsilon, and the Boltzmann constant = 1. The masses, distances, energies you specify are multiples of these fundamental values. The formulas relating the reduced or unitless quantity (with an asterisk) to the same quantity with units is also given. Thus you can use the mass & sigma & epsilon values for a specific material and convert the results from a unitless LJ simulation into physical quantities. 对于lj格式，所有数量都是无量纲的。不失一般性的，lammps设定了基本量质量sigma epsilon 和波尔兹曼常数=1。你定义的质量，距离，能量是这些基本值的倍数。公式把衰减量或者是无量纲量（用*）和单元给出的相同的量关联起来。因此你可以用mass，sigma，epsilon值表示一个具体材料并且改变一个无量纲lj模拟的物理量的结果。 mass = mass or m distance = sigma, where x* = x / sigma time = tau, where tau = t* = t (epsilon / m / sigma^2)^1/2 energy = epsilon, where E* = E / epsilon velocity = sigma/tau, where v* = v tau / sigma force = epsilon/sigma, where f* = f sigma / epsilon torque = epsilon, where t* = t / epsilon temperature = reduced LJ temperature, where T* = T Kb / epsilon pressure = reduced LJ pressure, where P* = P sigma^3 / epsilon dynamic viscosity = reduced LJ viscosity, where eta* = eta sigma^3 / epsilon / tau charge = reduced LJ charge, where q* = q / (4 pi perm0 sigma epsilon)^1/2 dipole = reduced LJ dipole, moment where *mu = mu / (4 pi perm0 sigma^3 epsilon)^1/2 electric field = force/charge, where E* = E (4 pi perm0 sigma epsilon)^1/2 sigma / epsilon density = mass/volume, where rho* = rho sigma^dim atom_style atomic 原子类型atomic（原子的) Define what style of atoms to use in a simulation. This determines what attributes are associated with the atoms. This command must be used before a simulation is setup via a read_data, read_restart, or create_box command. 定义在模拟中用到的原子类型，这个决定了与原子有关的属性，这个命令必须用在模拟被设定之前通过read_data, read_restart, or create_box Atomic only the default values coarse-grain liquids, solids, metals Atomic 附加属性是只是默认值用于粗粒液体，固体和金属的模拟
相变分析数值模拟实例详解 实例——水结冰过程分析 1、问题描述 有一圆柱体水缸，缸内盛放着高度为100mm的水，如图7.16所示，缸内水的初始温度为0℃，周围空气温度为-10℃，对流系数为12.5W/(m2?℃)，水的热性能参数见表7.3（水缸材料对水温的影响忽略不计）。试求： 在时间t=30分钟时，水或冰的温度场分布； 在时间t=120分钟时，水或冰的温度场分布，并绘制X轴和Y轴上各点温度随距离的变化关系曲线；绘制A、B、C、D各点温度随时间的变化关系曲线。 表7.3 水热性能参数 温度℃密度 Kg/m3 导热系数 W/(m?℃) 焓 J/m3 -10 1000 0.6 0 -1 1000 0.6 3.78e7 0 1000 0.6 7.98e7 10 1000 0.6 1.22e8 图8.16 水缸纵截面示意图 2、三维建模 应用Pro-E软件对流体计算域进行三维建模，实体如图7.17所示：
图7.17 水缸三维实体图 3、网格划分 采用流动传热软件CFX的前处理模块ICEM对计算域进行网格划分，得到如图7.18所示的六面体网格单元。流场的网格单元数为1920，节点数为2511。 图7.18 水缸网格图 4、模拟计算结果及分析 采用流动传热软件CFX非稳态计算，定义圆柱水缸内水的热传导系数为0.6 W/(m?℃)，水的初始温度为0℃，周围空气对流传热系数为12.5W/(m2?℃)，空气温度为-10℃。求解时选取Thermal Energy传热模型。水缸表面边界条件为-10℃温度载荷。求解方法采用高精度求解，时间步长为60s，总的时间为7200s。计算收敛残差为10-4。 图7.19为t=30分钟时，温度场分布等值线图 图7.20为t=120分钟时，温度场分布等值线图 图7.21为X轴上各点温度随距离的变化关系曲线；
Project #5 熔化与凝固：氩，铜，铝 铜和铝的熔化转变: 对于铜和铝，LAMMPS建立8×8×5的FCC晶格体系;充分弛豫后利用Nose-Hover方法，保持压强为零，使体系从T=2.5K开始加热，直至发生熔化转变。 下面是铜熔化的输入文件： # LAMMPS Melt_Cu or Al units metal # 单位，指定为lammps里的金属类的单位，长度为?，能量为eV。boundary p p p # 周期性边界条件 atom_style atomic # 原子模式 variable x equal 2.5 # 定义变量 x为初始温度 lattice fcc 3.61 # Cu 的晶格常数3.61 #lattice fcc 4.05 region box block 0 8 0 8 0 5 # x,y,z各方向上的晶胞重复单元数，也即区域大小 create_box 1 box # 将上述区域指定为模拟的盒子 create_atoms 1 box # 将原子按晶格填满盒子 timestep 0.01 # 步长 0.005fs thermo 1000 #每隔1000步输出热力学结果 pair_style eam/alloy # 选取 Cu 的EAM势作为模型 pair_coeff * * jin_copper_lammps.setfl Cu # EAM 势文件名称 #pair_style eam/fs #pair_coeff * * Al_FM.eam.fs Al neighbor 0.5 bin neigh_modify every 5 delay 0 check yes #velocity all create $x 825577 dist gaussian fix 1 all nvt $x $x 1.0 drag 0.2 #保持初始温度，在NVT下弛豫 #compute 3 all pe/atom
热应力数值模拟分析实例详解 实例1——圆筒热应力分析 1、问题描述 有一短圆筒，其横截面结构如图7.24所示，筒内避温度为200℃，外壁温度为20℃，圆筒材料参数如表7.4所示，求圆筒内的温度场、应力场分布。 表7.4 材料性能参数 弹性模量E GPa 泊松比ν 线膨胀系数α ℃-1 导热系数K W/(m?℃) 220 0.28 1.3e-6 70 图8.24 圆筒横截面结果示意图 2、三维建模 应用Pro-E软件对固体计算域进行三维建模，实体如图7.25所示： 图7.25 短圆筒三维实体图 3、网格划分 采用采用ANSYS有限元分析软件对计算域进行网格划分，得到如图7.26所示的六面体网格单元。流场的网格单元数为5760，节点数为7392。
图7.26 短圆筒网格图 4、模拟计算结果及分析 采用ANSYS有限元分析软件稳态计算，设置短圆筒导热系数为70W/(m?℃)，弹性模量为220Gpa，泊松比为0.28ν，线膨胀系数为1.3e-6℃-1。筒内壁加载温度载荷为200K，筒外壁加载温度载荷为20K。求解时选取Thermal Energy传热模型。求解方法采用高精度求解，计算收敛残差为10-4。 图7.27为圆筒内的温度场分布等值线图； 图7.28为圆筒轴截面上的温度场分布等值线图； 图7.29为圆筒轴截面上的径向应力场分布等值线图； 图7.30为圆筒轴截面上的轴向应力场分布等值线图； 图7.31为圆筒轴截面上的周向应力场分布等值线图； 图7.32为圆筒轴截面上的等效应力场分布等值线图。数据文件及结果文件在heat stress文件夹内。 图7.27 圆筒内的温度场分布等值线图