当前位置:首页 >> 材料科学 >> 计算材料学

计算材料学


计算材料学
杨振华

第一性原理计算方法
?
第一性原理方法是一种理想的研究方法,物理学家常 称第一性原理方法,化学家常称为“从头算”,但是 本质都是一样的。就是从材料的电子结构出发,应用 量子力学理论,只借助于普朗克常数h、电子的静止 质量m0、电子电量e、光速c和波尔兹曼常数k这五个 基本的物理常量,以及某些合理的近似而进行计算。 这种计算不需要任何其他可调的(经验的或拟合的)参 数就可以如实地求解材料的一些基本物理性能参数。 通过求解多粒子系统总能量的办法来分析体系的电子 结构和原子核构型的关系,从而确定系统的性质 。

绝热近似
? 波恩(Born M)和奥本海默(Oppenheimer J.E) 提出了绝热近似,根据这种近似,可 以将原子核运动和电子的运动分开。通过 绝热近似,可以获得多电子的薛定谔方程

H? (r, R) ? E ? (r, R)
H

H ? H e ? H N ? H e? N
电子作 用项 原子核 作用项 电子和原子核 相互作用项

? 波恩(Born M)和奥本海默(Oppenheimer J.E) 提出了绝热近似 单粒子算 双粒子算符


? ? ? 1 ' 1 ? 2 ? ? ? V ( r ) ? ? ? H ? H ? ? ri ? i ? ?? i ? ii' ?? ? E? ? 2 i ,i ' | ri ? ri ' | ? ? i ,i ' ? i ? ? ?

多电子的薛定谔方程
,成功地分开了电子的运动与 原子核的运动

? 哈特利方程
' 2 ? ? | ? ( r ) | 2 ' i' ?? ? ? V (r ) ? ? ? dr ?? i (r ) ? Ei? i (r ) ' |r ?r| ? i ' ( ?i ) ? ? ?

此方程以位于r处的单个电子为研究对象, 描述其在晶格势和其他所有电子的平均 势中的运动规律

? 将多电子问题变为了单电子问题,但是没 有考虑电子的交换反对称性 。为了研究电 子的交换反对称性的影响,采用Slater行 列式来求能量,经过合适的变换,得到了 如式所示方程:
' 2 * ' ' ? | ? ? ' (r ) | ? ' ( r )? i ( r ) 2 ' ' i i ? ( r ) ? dr ? ? Ei? i (r ) ?? ? ? V (r ) ? ? ? dr ? ? i ' ' ? |r ?r| ? |r ?r| i ' ( ?i ) i ' ( ? i ),|| ? ? ?

单电子的哈特利-福克方程, 比哈特利方程多了交换相互作用项。

? 多电子的薛定谔方程可通过哈利特-福克近 似简化为单电子有效方程,如式所示。 ?

?? ?

2

? Veff (r ) ?i (r ) ? Ei?i (r )

?

包含了电子与电子的交换相互作用,但自旋 反平行电子间的排斥相互作用没有被考虑, 即还需考虑电子关联相互作用。

? 为了更加准确地描述多电子系统, Hohenberg P和Kohn W提出了两个基本 的定理: ? (1) 定理1:不计自旋的全同费密子系统的 基态能量是粒子数密度函数的唯一泛函; ? (2) 定理2:能量泛函在粒子数不变条件下 对正确的粒子数密度函数取极小值,并等 于基态能量。

? 定理1的主旨思想是粒子数密度函数是一个 决定系统基态物理性质的基本变量;定理2 的要点是在粒子数不变条件下能量泛函对 密度函数的变分就得到系统基态的能量。 密度泛函理论的理论基础是这两条基本定 理,其基本的思想是原子、分子和固体的 基态物理性质可以用粒子密度函数来表示。

? Hohenberg-Kohn定理说明了粒子数密度 是确定多粒子系统基态物理性质的基本变 量以及能量泛函对粒子数密度函数的变分 是确定系统基态的途径。但是仍然存在三 个问题未解决: ? (1) 如何确定粒子数密度函数; ? (2) 如何确定动能泛函; ? (3) 如何确定交换关联能泛函。

? 为了解决这三个问题,Kohn W与 Sham L.J共同合作,提出了Kohn- Sham方程 。

?? ?

2

? VKS [ ? (r )] ?i (r ) ? Ei?i (r )
N 2 i ?1

?

? (r ) ? ? | ?i (r ) |
N i ?1

Ts [ ? ] ? ? dr?i* (r )(?? 2 )?i (r )

VKS [ ? (r )] ? v(r ) ? Vcoul [ ? (r )] ? Vxc [ ? (r )]

? Kohn W和Sham L.J成功地提出了KohnSham方程,用无相互作用的粒子模型代替有相 互作用粒子哈密顿量中的相应项,将有相互作用 粒子的全部复杂性归入交换关联作用泛函。将多 粒子系统的基态求解转化为单粒子系统的等效求 解,解决第一和第二个问题,对于第三个问题, 需要采用局域密度近似来解决。为了求解KohnSham方程,必须构造合适的交换关联能。目前 比较常用的交换关联能主要有以下两种形式:局 域密度近似(LDA)和广义梯度近似(GGA)。

局域密度近似
? 局域密度近似最早是由Kohn W和Sham L.J提出来的,这是一种既简单可行而又很 有效的近似,其基本思想是在局域密度近 似中,利用均匀电子气密度函数来获得非 均匀电子气的交换关联泛函。 ? 交换关联能可以写为式

E xc [ ? ] ? ? dr ? (r )? xc [ ? ( r )]

? Kohn- Sham方程中的交换关联势近似为 式
密度为 ? (r )

?E xc [ ? ] d Vxc [ ? (r )] ? ? ( ? (r )? xc [ ? (r )]) ?? d? (r )
? xc [ ? (r ) :均匀无相互作用电子气的交换-关联密

度,在实际的计算过程中,通常把交换关联密度分成两部分:交换项和关联项。

? xc [ ? (r )] ? ? x [ ? (r )] ? ? c [ ? (r )]
交换能 关联能

LSDA E xc [ ? ? ? ? ? ] ? ? d 3 r? (r )? xc ( ? ? (r ), ? ? (r ))

? ? d 3 r? (r )[? x ( ? ? (r ), ? ? (r )) ? ? c ( ? ? (r ), ? ? (r ))]
考虑了自旋

? Local Density Methods

假设局域电子密度可以被认为是均匀电子气,或等效地说,电子密度是随空间缓慢 变化的函数。
交换项

Local Density Approximation (LDA)
4/3 LDA 1/ 3 ExLDA ? ? ? ? r ?? ? ? ?Cx ? ? ? r ? dr?????????? x ? ? ? ? r ?? ? ? ?Cx ?

Local Spin Density Approximation (LSDA)
1 4/3 4/3 4/3 4/3 LSDA 1/ 3 ? ? ? ExLSDA ? ? ? ? ?21/ 3 Cx ? ? ? ? ? d r ????????? ? ? ? ? C ? 1 ? ? ? 1 ? ? ? ? ? ? ? ? ? ? x x ? ? ? ? 2

关联项 Vosko,Wilk,and Nusair (VWN)

?

VWN c

? f ?? ? ? 4 4 ? 1 ? ? ? ? r ,1 ? ? r , 0 f ? ? ? ? ? rS , ? ? ? ? c ? rS , 0 ? ? ? a ? rS ? ? '' ? ? ? ? ? ? ? ? c S c S ? ? ? ? ? f ? 0? ?

2 ? bx0 ? ? x ? x0 ? 2 ? b ? 2 x0 ? ?1 ? Q ? ? ? x2 2b ?1 ? Q ? ? ? ? c / a ? x ? ? A ?ln ? tan ? ? ln ? tan ? ? ? ? ? ? X x Q 2 x ? b X x X x Q 2 x ? b ? ? ? ? ? ? ? ? ? ?? 0 ? ? ? ?? ? ?

?1 ? ? ? f ?? ? ?

4/3

2?2

? ?1 ? ? ?
1/ 3

4/3

? 1?

?2

x ? rS ????? X ? x ? ? x 2 ? bx ? c?????Q ? 4c ? b 2

GGA (见下)中的 PW91 修改了 VWN 的泛函形式:

?

PW 91 c/a

? ? 1 ? 1? ? x ? ? ?2a? ?1 ? ? x ? ln ? ? 2a ? ?1 x ? ? 2 x 2 ? ?3 x3 ? ? 4 x 4 ? ? ? ?
2

? Gradient Corrected Methods

Gradient Corrected or Generalized Gradient Approximation (GGA): 泛函不仅决定 于电子密度,还决定于电子密度的梯度。
交换项 Perdew and Wang (PW86): 修正 LSDA 的泛函形式:加入高阶项。

?

PW 86 x

??

LDA x

?1 ? ax

2

? bx ? cx
4

6 1/15

?

?????x ?

??

? 4/ 3

Becke (B or B88): 正确的能量密度渐进行为。

? xB88 ? ? xLDA ? ?? xB88 ??????? xB88

2 x ? ??? 1/ 3 1 ? 6? x sinh ?1 x

Becke and Roussel (BR): 加入轨道波函数的导数项。

?

BR x

2 ? 2 e? ab ? ab e? ab ?? ?????a3 e? ab ? 8?? 4b ?
2

a ? ab ? 2 ? ? b

?? ? ? ? ? 2D 2 ?????D ? ? ??i ? ? 4?
N i

2

Perdew and Wang (PW91)

? ?? xa sinh ?1 ? xa ? ? a ? a e?bx2 x 2 ? 1 2 3 4 ??????x ? ?? ? xPW 91 ? ? xLDA ? ? ? 1 ? xa1 sinh ?1 ? xa2 ? ? a5 x 2 ? 4/ 3 ? ? ? ?
关联项

?

?

Lee, Yang, and Parr (LYP)
?
LYP c

? ? e? c? 8/ 3 8/ 3 ? ? ?18 ? 22 / 3 ? CF ? ?? ? ?a ? ab ? ?? ? 18? tW ? ?? ? 2tW ? ? 2 ?? ? ? ? ? ? 2tW ? ? 2 ? ? ?? ? ?1/ 3 ?1/ 3 8/ 3 ? ? ?1 ? d ? ? 9 ?1 ? d ? ? ?
2 2 ? ?? ? ? ?? 1 ? ??? ? ? ? 2 ?1 ? ? ????????????tW ? ? 2 ? 8 ? ?? ? ? ? ? 2

?1/ 3

? ? ? ?? ? ? ?
2

Perdew:修正 LSDA 的梯度项。

? cP86 ? ? cLDA ? ?? cP86 ??????? cP86
f ?? ? ? 2
1/ 3

e? C ? ? ? ?? C ? ? ? ?? ? ?????? ? a 7/3 f ?? ? ? C ? ? ? ? 7/6
2
5/ 3

? 1? ? ? ? ? ? 2 ?

5/ 3

? 1? ? ? ?? ? ? 2 ?

b2 ? b3rS ? b4 rS2 ?????C ? ? ? ? b1 ? 1 ? b5 rS ? b6 rS2 ? b7 rS3

Perdew and Wang(PW91 or P91):。

?? cPW 91 ? ? ? H0 ?t, rS , ? ? ? H1 ?t, rS , ? ?? ? t 2 ? At 4 ? 3 ?1 H 0 ? t , rS , ? ? ? b f ?? ? ln ?1 ? a 1 ? At 2 ? A2t 4 ? ? ? 3 2 ? dx 2 / f ?? ?2 ? 16 ? 2 1/ 3 H1 ? t , rS , ? ? ? ? ? ? 3? ? ? C ? ? ? ? c? f ?? ? t e ? ? ?? ? 1/ 6 ?? ?? 3 ? 192 ? t ?? 2 ? ??????? ? a ?exp ?b? c ? rS , ? ? / f ?? ? ? 1? ?????? 7/6 ? ? ? ? ? 2 f ?? ? ?

?

?

1 2/3 2/3 f ?? ? ? ?1 ? ? ? ? ?1 ? ? ? 2

?

?

其中

? c ? rS , ? ? 在 LSDA 部分已经给出。

Becke(B95):更好地满足一些基本的物理约束。
?? ?? ? cB95 ? ? c ? ?c ? ? c??
??
2 2 ?1

PW 91,?? ? ?c ? ? 1 ? a x ? x ? ? ? ? ?? c ? ?? 2 ?2 D? PW 91,?? ? ?c ? ? 1 ? bx ? ?? LDA c ?

D?

LDA 5/3 D? ? 25/3 CF ??

? 混合方法 混合 HF 和 DFT 给出的能量项。

Becke 3 parameter functional (B3)
B3 GGA Exc ? ?1 ? a ? ExLSDA ? aExHF ? b?ExB88 ? EcLSDA ? c?Ec

广义梯度近似
? 为了对局域密度近似进行提高和改善,引入 了电荷密度梯度,即粒子密度的空间分布不 仅仅与局域密度有关系,而且与对应点附近 的密度有关系。其中最为常用的是广义梯度 近似(GGA),在GGA近似下,在交换相关 能泛函中引入电子密度的梯度来完成。

E xc [ ? ] ? ? d r ? (r )? xc [ ? (r ), | ?? (r ) |]
3
考虑了电子密度的非局域性,改善了 LDA的计算结果。一般GGA的计算结 果与实验结果较为吻合。

DFT+U方法简介
? 基于密度泛函理论(DFT)的第一性原理计算 方法已在材料的晶体结构、磁结构、电子 结构以及材料的力学性能计算等方面取得 了巨大的成功,但是对于Mott绝缘体(如过 渡金属氧化物和稀土氧化物),由于其d电 子或f电子的强关联作用,传统的第一性原 理方法已不能很好地描述其基本性质。

? 在Mott绝缘体中,当电子从一个一个原子 位置跳跃到另外一个原子位置时,如果那 个原子位置已经拥有一个电子,电子之间 就会产生库伦排斥力作用,这种跳跃需要 一定的能量以致能克服这种库伦斥力作用, 如果这个能量大于能带带隙,即使能带没 有全部占满,电子也很难自由输运,从而 使材料体现绝缘体的特征。

? 当采用传统的第一性原理计算Mott绝缘体时,只 考虑了交换参数J,没有考虑Hubbard参数U,而 在Mott绝缘体中,其决定性的参数是Hubbard参 数U值,因此采用传统的计算方法往往会导致失 败。为了解决计算Mott绝缘体的问题, Anisimov等提出了Anisimov 模型,在该模型中, 将所研究的电子分为两个部分:(1) 传统的DFT 算法,在此过程中没有考虑Hubbard参数U;(2) 对于d轨道电子或f轨道电子,能带模型为 Hubbard模型,考虑了d轨道或f轨道电子的强关 联作用。

? LDA+U方法为例,电子的总能量计算可 考虑了d轨道或f轨道电子的强关 以通过下式进行表述: 联作用,并采用Hartree表达式
所计算的能量

LDA E[?, ?ni ?] ? ELDA [? ] ? EH [?ni ?] ? Edd [nd ]

局域态的轨 道占据数

总的局域 电子数

原来传统LDA 计算过程所包 含的关联能, 采用LDA+U 方法后,此项 应该减去

nd ? ? ni
1 E H [?ni ?] ? U ? ni n j 2 i? j
E
LDA dd

1 [nd ] ? Un d (nd ? 1) 2

U 为Hubbard参数

E[ ? , ?ni ?]

对轨道占据数进行微分
当轨道占据数分别为1和0时, 相应的值表示将采用传统LDA 计算所得的轨道能量分别偏移

?E[ ? , ?ni ? 1 ?i ? ? ? LDA ? U ( ? ni ) ?ni 2
电子轨道 势

1 Vi (r ) ? VLDA (r ) ? U ( ? ni ) 2

VASP计算软件包简介
? VASP,其全称是Vienna Ab-initio Simulation Package。它基于1989年的CASTEP(1989版), 最早是由Gerorgo Kresse 和Jürgen Furthmüller合作,共同开发出来的,1995年被 正式命名为VASP,随后被开发者不断完善。 VASP是一种使用赝势和平面波基组进行从头量子 力学分子动力学计算和第一性原理计算的软件包, 主要用于具有周期性的晶体或表面的计算,可以 采用大单胞,也可以用于处理小的分子体系。

? 与同类的软件相比,它比较早地实现了超 软赝势,计算量相对于一般的模守恒赝势 方法大为减少。其对计算领域最大贡献无 疑是在Bl?chl的基础上发展的投影缀加平面 波(PAW)方法,这是最重要的。这使得 VASP不仅计算速度快,而且精度是abinit 和pwscf没法比的。VASP的精度,比如磁 性计算,很多可以跟FLAPW相比,并且计 算速度比FLAPW快很多。

? 在实空间计算势的非局域部分并保持正交 化的数目减少,使得计算时间小于N3; VASP在电子自洽迭代计算中,采用了 RMM-DISS和blocked Davidson等非常 有效算法并能自动确定体系的对称性;此 外,VASP的代码使用FORTRAN语言编写, 可读性好,几乎支持所有的计算机平台, 已广泛应用于材料科学领域。

?VASP基本原理简介

?基本知识
?常用关键词使用说明

?计算结果处理

VASP程序基本原理
VASP是基于赝势平面波基组的密度泛函程序,其前身
是CASTEP 1989版本,其基本原理如下: 根据Bloch定理,对于周期体系,其电子波函数可以写 为单胞部分和类波部分的乘积:
?? ? ? ik ?r ? i (r ) ? e f i (r )

其中,单胞部分的波函数可以用一组在倒易空间的平面
波来表示:
?? ? i G ? e ?r f i (r ) ? ? c i ,G ? G

这样,电子波函数可以写为平面波的加和:
? ? ? ? i ( k ? ? e ?G )?r ? i (r ) ? ? c i , k ?G ? G

根据密度泛函理论,波函数通过求解Kohn—Sham方程来确定:

?2 2 ? ? ? ? ? [? ? ? Vion (r ) ? VH (r ) ? V XC (r )] ? i (r ) ? ? i? i (r ) 2m ?i:Kohn—Sham本征值
Vion:电子与核之间的作用势 VH和VXC:电子的Hartree势和交换—相关势

? n ( r ') 3 ? ? 2 VH (r ) ? e ? ? ? d r ' | r ? r '|

? ? ?E XC [n(r )] V XC (r ) ? ? ?n(r )

基于平面波表示的Kohn—Sham方程:
? ? ? ? ? ? ? ? ?2 ? ? 2 ? ? ? V (G ? G ' ) ? V (G ? G ' ) ? V ? ? ? ? | k ? G | ?G ? ion H XC (G ? G ' ) ? ci , k ? G ? ? i ci , k ? G G' ? ? G ' ? 2m ?

上式中动能项是对角化的,通过求解上式方括号中的哈密顿矩

阵来求解KS方程,该矩阵的大小由截至能(cutoff energy)来决定。

尝试电子密度和尝试波函数

程序流程:
写出交换相关势表达式

构造哈密顿量

子空间对角化,优化迭代

自由能的表达式E

新电子密度,与尝试电子密度比较 是



输出结果,写波函数

与原子轨道基组相比,平面波基组有如下优点:

1) 无需考虑BSSE校正;
2) 平面波基函数的具体形式不依赖于核的坐标,这样,一 方 面 , 价 电 子 对 离 子 的 作 用 力 可 以 直 接 用 H el lm a n -

Feymann 定理得到解析的表达式,计算显得非常方便,
另一方面也使能量的计算在不同的原子构象下具有基本 相同的精度;

3) 很方便地采用快速傅立叶变换 (FFT) 技术,使能量、力
等的计算在实空间和倒易空间快速转换,这样计算尽可 能在方便的空间中进行; 4) 计算的收敛性和精确性比较容易控制,因为通过截断能 的选择可以方便控制平面波基组的大小。

平面波基组方法的不足之处: 1) 所求得的波函数很难寻找出一个直观的物理或化学图象与

化学家习惯的原子轨道的概念相联系,即其结果与化学家
所感兴趣的成键和轨道作用图象很难联系出来,这就为我 们计算结果的分析带来了困难; 2) 考察某些物理量时,例如原子电荷,涉及到积分范围的选 取,这造成所得物理量的绝对值意义不大;

3) 有些方法,例如杂化密度泛函方法不易于采用平面波基组
方法实现。

VASP程序基本知识
1. VASP程序主要功能:
1) 能量计算

J. Phys. Chem. C, 2008, 112, 191

2) 电子结构(能带结构、DOS、电荷密度分布)

能带结构

DOS

电荷密度分布

J. Phys. Chem. B, 2005, 109, 19270

3) 构型优化(含过渡态)和反应途径

J. Phys. Chem. B, 2006, 110, 15454

4) 频率计算和HREELS能谱模拟

J. Phys. Chem. C, 2007, 111, 7437

5) STM图像模拟

Surf. Sci., 2007, 601, 3488

6) UPS能谱图像模拟

Surf. Sci., 2007, 601, 3488

7) 材料光学性质计算

8) 其它性质计算,包括功函、力学性质等

2. 重复平板模型(或层晶模型):
VASP程序采用重复平板模型来模拟零维至三维体系

零维分子体系

Dv: Vacuum thickness (~10 A)
二维固体表面

说明: 重复平板模型中的平移矢量长度必须合理选择,以保证: 1) 对于分子体系,必须保证相邻重复单元中最近邻原子之 间的距离必须至少7~10埃以上; 2) 对于一维体系,相邻两条链最近邻原子之间的距离必须 至少7~10埃以上; 3) 对二维体系,上下两个平板最近邻原子之间的距离必须 至少7~10埃以上;

4) 严格意义上,通过考察体系总能量/能量差值对真空 区大小的收敛情况来确定合理的平移矢量长度。

tal energy

3. K网格大小的选择:
对于一维至三维体系的计算,需涉及k点数目的选择,对

于K点的确定,它与布里渊区的形状以及对称性有关。VASP的
K点输入方法有多种,其中最常用的是直接给定K-mesh的大小, 然后程序根据布里渊区的形状以及对称性自动生成各K点的坐 标和权重。 对于K-mesh的确定方法,通常通过考察总能量/能量差的收敛

程度来确定,能量的收敛标准是1meV/atom。
多数情况下,对半导体或绝缘体较小的K-mesh能量就可以 收敛,对于导体,一般需要较大的K-mesh。

-10.2

-10.3

-10.4

Total energy(eV)

-10.5

-10.6

-10.7

-10.8

-10.9 2 4 6 8 10 12

Size of k-mesh 硅体相总能量随K-mesh 大小的变化情况

4. Cutoff energy大小的选择:
截至能的大小直接影响到计算结果的精度和计算速度,

因此,它是平面波计算方法的一个重要参数。
理论上截断能越大计算结果也可靠,但截至能大小决定 了计算中平面波的数目,平面波数目越多计算时间约长、内 存开销越大。 一般根据所求物理量来确定截至能,例如计算体模量以

及弹性系数时,需要较高的截至能,而通常的构型优化只要
中等大小的截至能即可,另外动力学模拟时,可选取低的截 至能。

不同元素在构造其赝势时,有各自的截至能,对于VASP, 在缺省情况下,选取的是中等大小的截至能,这对于求解多

数物理量是足够的。严格意义上,截至能的确定与K-mesh大
小的确定类似,也是通过考察在总能量的收敛情况来确定(即 保证总能量收敛至1meV/atom)。

-10.55

-10.60

Total energy (eV)

-10.65

-10.70

-10.75

-10.80

硅体相总能量随cutoff energy 大小的变化情况 Cutoff energy (eV)

100

150

200

250

300

5. VASP输入和输出文件:
输入文件(文件名必需大写) INCAR : 其内容为关键词,确定了计算参数以及目的; POSCAR : 构型描述文件,主要包括平移矢量、原子类 型和数目、以及各原子坐标; KPOINTS : K点定义文件,可手动定义和自动产生; POTCAR : 各原子的赝势定义文件。

主要输出文件 OUTCAR : 最主要的输出文件,包含了所有重要信息; OSZICAR : 输出计算过程的能量迭代信息; CONTCAR: 内容为最新一轮的构型(分数坐标,可用于续算);

CHGCAR、CHG、PARCHG :用于电荷密度图绘制;
WAVECAR : 波函数文件; EIGENVAL: 记录各K点的能量本征值,用于绘制能带图; XDATCAR: 构型迭代过程中各轮的构型信息(分数坐标,用于 动力学模拟);

DOSCAR : 态密度信息。

POSCAR文件内容说明:
Silicon bulk (Title) 2.9 (Scaling factor or lattice constant) 0.0 1.0 1.0 (第一个平移矢量的方向)

1.0
1.0

0.0
1.0

1.0 (第二个平移矢量的方向)
0.0 (第三个平移矢量的方向)

2(单胞内原子数目以及原子种类) Selective dynamics(表示对构型进行部分优化,如果没这行,则表示全优化) Direct (表示所采用的为分数坐标,如果内容为Car,则坐标单位为埃) 0.125 0.125 0.125 T T T (各原子坐标以及哪个方向坐标放开优化)

-0.125 -0.125 -0.125 T T T

surface of mgo(100) (2*2)Mg
1.00000000000000 5.9459999999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 5.9459999999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 20.0000000000000000 20 20 (体系中有2种元素,各自的原子数目分别为20,20) Selective dynamics Direct 0.0000000000000000 0.0000000000000000 0.0000000000000000 F F 0.5000000000000000 0.0000000000000000 0.0000000000000000 F F 0.5000000000000000 0.5000000000000000 0.0000000000000000 F F 0.0000000000000000 0.5000000000000000 0.0000000000000000 F F …… 0.2500000000000000 0.2500000000000000 0.0000000000000000 F F 0.7500000000000000 0.2500000000000000 0.0000000000000000 F F 0.2500000000000000 0.7500000000000000 0.0000000000000000 F F 0.7500000000000000 0.7500000000000000 0.0000000000000000 F F

F F F F
F F F F

……

POTCAR文件内容说明: VASP程序本身有提供了赝势库,只需将体系各类原子的 赝势合并在一起即可,但需注意到: 1) 赝势类型: LDA US型赝势 GGA PW91 US 型赝势所需截至能 较小,计算速度快, PAW 赝 势 截 至 能 通 常 较大,而且考虑的电子 数多,计算慢,但精确 度高。

PBE
LDA PAW型赝势 GGA PBE PW91

2) POTCAT中各原子赝势定义的顺序必需与POSCAR中相同:
surface of mgo(100) (2*2)Mg 1.00000000000000 5.9459999999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 5.9459999999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 20.0000000000000000 20 20 Selective dynamics Direct ……

3) 对各原子的赝势参数,我们最关心的是截至能以及电子数; 4) POTCAR的泛函类型必需与INCAR中GGA关键词定义的 类型一致; 5) 使用zcat命令产生和合并POTCAR文件。

电子数目和组态

构造该赝势时,所采用的泛函类型, 这里为PW91

对应于中等大小的截至能 (构型优化时采用)

对应于低的截至能 (动力学模拟时采用)

KPOINTS文件内容说明: 一般有两种定义K点的方法:

1) 通过定义K-mesh大小,由程序自动产生各K点: Automatic mesh (title)
0 (为0时,表示自动产生K点) M (表示采用Monkhorst-Pack方法生成K点坐标)

5 5 5(对应于5x5x5网格)
0 0 0(原点平移大小)

2)手动定义各K点的坐标(一般仅在计算能带结构时使用): k-points for MgO(100) (title) 31 (K点数目) Rec (字母R打头表示为倒易空间坐标,否则为实空间的坐标) 0.0 0.0 0.0 1.0 (各K点的坐标以及权重) 0.05 0.0 0.0 1.0 0.1 0.0 0.0 1.0 0.15 0.0 0.0 1.0 0.2 0.0 0.0 1.0 0.25 0.0 0.0 1.0 0.3 0.0 0.0 1.0 0.35 0.0 0.0 1.0 0.4 0.0 0.0 1.0 0.45 0.0 0.0 1.0 0.5 0.0 0.0 1.0 ……

6. VASP安装和运行:
(1) VASP程序安装:

a. 设置编译环境:安装Fortran编译器,常用为IFC
b. 对于并行版本vasp的编译,还需安装MPICH c. 编译vasp自带的库文件 d. 对makefile进行修改,包括BLAS和Lapack库文件所在 目录,一般可采用IFC所带的数学库 e. 运行make命令进行编译 (2) 创建输入文件,包括INCAR,KPOINTS,POSCAR

和POTCAR

(3) 运行vasp:
单机版: ~/bin/vasp.4.5-ifc-mk-sp > vasp.out & 版本号 编译环境 多个K点Single process CPU数目

并行版本:

mpirun –np 4 –machinefile ./hosts ~/bin/vasp.4.5-ifc-mk-mp > & vasp.out & 存放要并行运算的机器名 或者IP

常用关键词使用说明
(部分参考清华大学物理系苏长荣编写的VASP安装和使用说明)

(1)

(2)

(3)

一般单胞尺寸大时,选实空间,小单胞选取倒易空间。 (4)

(5) ENCUT=数值 用户手动定义截至能,如果没有,则由PREC选项确定。 (6) EDIFF=1e-4

EDIFFG=EDIFF× 10 当数值为负数时,表示以力作为收敛标准,多数情况均采用 力作为收敛标准。
(7) ALGO=38|48 该关键词确定能量计算迭代方法 38-采用Davidson优化方法;(可靠,但速度慢) 48-采用RMM-DIIS算法;(常用,速度快) (8) ISYM=0|1|2 该关键词确定能量和构型优化时是否使用对称性(将影响到K 点数目和计算量大小) 0-不使用对称性; 1-采用对称性; 2-用于PAW型赝势;

(9) NELM=整数 该关键词确定能量自洽场最大迭代轮数,缺省为60轮; NELMIN=整数 在构型优化中,计算每个构象能量时最少 迭代轮数,一般为3~4,以保证能量和力的稳定性; (10) 定义DFT泛函类型,注意要与 POTCAR 中的赝势类型一致。

(11) ISPIN=1|2 1-非自旋极化计算(缺省) 2-自旋极化计算, 将给出体系磁矩大小(对含有过渡金属原 子体系,一般均要采用自旋极化方法)。

(12)

(13)

ISMEAR选择: 1) 对半导体或绝缘体选取-5,如果单胞较大时,或者所选取k 点数目少时,用0; 2) 对导体,通常用0; SIGMA取值: SIGMA取值的原则是使得计算得到的TS项(OUTCAR中),

分摊到每个原子上时小于1meV,否则得到的总能量不准确,
对导体尤其要注意该参数的选择。

(14) 以下为构型优化所用关键词: ?NSW=整数 构型优化的最大轮数

?IBRION = -1|0|1|2
构型优化方法: -1-构型不变更; 0-分子动力学模拟; 1-采用准牛顿方法确定新的构型(当初始构型较合理时使用); 2-采用CG方法确定构型(当初始构型离平衡位置较远时使用)。 ?POTIM=数值 控制构型优化步长,缺省为0.5,对动力学模拟则为时

间步长(单位为fs)

(15) 输出控制关键词:
LCHARG = .FALSE. (输出电荷密度?)

LWAVE = .FALSE. (输出波函数?)
LVTOT = .FALSE.(输出静电势,求功函时使用)

(16) 其他关键词: NPAR = 8 (CPU数目,并行计算时使用) LPLANE = .TRUE.(与并行算法有关)

实例: SYSTEM = Silicon bulk NPAR = 2 LPLANE = .TRUE. Elecronic minimisation ISTART = 0 LREAL= .FALSE. PREC = Medium precission: Mediun/High/Low EDIFF = 1e-4 converge criterion: default = 1e-4 EDIFFG = -0.02 converge criterion for relation loop IALGO = 48 algorithm (8-CG, 48-RMM) NELMIN = 3 the minimum number of electronic SC steps ISYM = 1 symmetry (2-PAW on, 1-US-PP's on, 0-off) ISIF = 3 Relax ions ISPIN = 1 ISMEAR = -5 (tetrahedron/gaussian/m-p) SIGMA = 0.1

OUTPUT CONTROL LCHARG = .FALSE. LWAVE = .FALSE. the i/o cost is not worth it. LVTOT = .FALSE. IONIC RELAXATION NBLOCK = 1 steps for inner block NSW = 300 number of steps for IOM IBRION = 1 -1:no update 0-MD 1-quasi-New 2-CG POTIM = 0.50 default 0.5 of IBRION=1-3

练习: 对MgO体相的构型(包括原子位置和单胞外形) 进行优化,并考察 k 网格和动能大小对计算结果 的影响,计算体系的力学性质。

MgO体相的构型参数: 空间群—No. 225 FM-3M 单胞参数—a=b=c=4.2112?

步骤: 1. 确定构型,包括平移矢量以及各个原子的坐标。 可通过其它软件(如Crystal, Materials Studio等)间接获得。

从MS构型库里面读取MgO构型:

转化为原胞构型:

原胞构型

选择Castep计算设置窗口,选择Files选项获取构型输入文件:

到相应的目录下找到 CELL 文件(该文件为隐藏文件 ),从中 找到平移矢量和各原子分数坐标:

2. 构造VASP四个输入文件,包括INCAR、KPOINTS、 POSCAR和POTCAR。 其中POSCAR文件内容由前一步骤获得: MgO bulk 1.0 0.0 2.1056 2.1056 2.1056 0.0 2.1056 2.1056 2.1056 0.0 11 direct 0.0 0.0 0.0 0.5 0.5 0.5 INCAR和KPOINTS文件内容见前.

3. 运行vasp 4. 力学性质中弹性系数计算方法见参考文献:Comput. Phys. Commun., 2010, 181, 671. 对立方晶系,独立的弹性系数有c11, c12和c44:

体模量和弹性系数的关系为:

实验值:c11=297.08 GPa; c12=95.36 GPa; c44=156.13 GPa

赝势· 选择

? POTCAR 赝势文件

? 需要按照POSCAR里的顺序,将各元素的 POTCAR按顺序连接起来就可以了如以下 命令: cat file1 file2 file3 > POTCAR
?

? 软件包自带的绝大多数赝势是超软赝势(US-PP) 了,但不少元素有两个版本,如何选取呢? ? ? 一个简单的办法是看后缀 标准的没有后缀 _h 硬一点 _s 软一点_pv,_sv,_d 就是说semicore的p,s或者d也当做价态处理了如果是数字的 话,表示的可能是不同的半径截距

k点的选择

? 一般来说,k点越密越多,计算精度也越高, 当然计算成本也越高。 嗯,对于k点的需求, 金属>>半导体,绝缘体,不过呢,很多时 候主要还是受硬件限制 简约化可以使k点的 数目大大下降。对于原子数较多的体系的 计算,就需要谨慎的尝试 k点数目,在避免 或者预先评估wrap-around error的前提 下尽量减少k点数目。

? 另一个问题是k空间网格(k-points grid) 的位置和形状, 是否包括Г点(Gamma点, 也可理解为原点)?(一般不包括的话很可能 会带来误差,尤其是使用tetrahedron方法 的时候。暂时还不知道不包括的好处,为 了减少k点?) 方形?线形?还是长方形? 或者奇形怪状?

? Line1: comment line 注释行 no problem Line2: k点总数 或者 '0'自动生 成网格(Automatic k-mesh generation) 如果是前者,给出k点总数,又分两种情况 M.全手动 Entering all k-points explicitly Line3: 输入格式标识。直角坐 标 (Cartesian)或者 倒格坐标 (Reciprocal)

VASP中电子态密度计算的流程
? 主要分成三步:一、结构优化;二、静态 自洽计算;三、非自洽计算

? ? ? ? ? ? ? ? ? ? ? ? ?

第一步 结构优化 输入文件(INCAR, POTCAR, POSCAR, KPOINT) INCAR文件 System=Al ISTART=0 ISMEAR=1 SIGMA=0.2 ISPIN=2 GGA=91; VOSKOWN=1; EDIFF=0.1E-05; EDIFFG=-0.01 IBRION=2 NSW=50 ISIF=2 (OR 3) NPAR=10

? ? ? ? ? ? ? ? ? ? ? ? ?

POTCAR 文件直接在势库中拷贝 POSCAR文件 Al 4.05 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 4 Direct 0.0 0.0 0.0 0.5 0.5 0.0 0.5 0.0 0.5 0.0 0.5 0.5

? KPOINT 文件 ? Automatic generation ?0 ? Mohkorst Pack ? 15 15 15 ? 0.0 0.0 0.0

? 第二步 静态自洽计算 ? INCAR: PREC = Medium,ISTART = 0,ICHARG = 2,ISMEAR = -5 ? 输入文件(INCAR, POTCAR, POSCAR, KPOINT) ? INCAR文件 ? System=Al ? ISTART=0 ? ISMEAR=1 ? SIGMA=0.2 ? ISPIN=2 ? GGA=91; VOSKOWN=1; EDIFF=0.1E-05; EDIFFG=-0.01 ? #IBRION=2 ? #NSW=50 ? #ISIF=2 (OR 3) ? NPAR=10

? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

POTCAR 文件直接在势库中拷贝 POSCAR文件 Al 4.05 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 4 Selective Dynamic Direct 0.0 0.0 0.0 T T T 0.5 0.5 0.0 T T T 0.5 0.0 0.5 T T T 0.0 0.5 0.5 T T T

? KPOINT 文件 ? Automatic generation ?0 ? Mohkorst Pack ? 15 15 15 ? 0.0 0.0 0.0

? ? 第二步计算是在结构优化的结果上进行的, 所以开始第二步的时候,将第一步中的输 入文件INCAR, POTCAR, POSCAR, KPOINT 以及 C* 文件放入静态自洽计算 中去,并且将CONTCAR 拷贝到 POSCAR中,然后运行VASP。计算结果 中的Fermi能是准确的,需要

? 第三步 非自洽计算 ? INCAR: PREC = Medium,ICHARG = 11,ISMEAR = -5,LORBIT = 10 或者 11(这时可不设RWIGS),ISTART = 1

? 在第二步自洽计算的基础上进行,修改输入文件INCAR, POTCAR, POSCAR, KPOINT。 ? INCAR文件 ? System=Al ? ISTART=1 ? ISMEAR=-5 ? SIGMA=0.2 ? ICHARG=11 ? RWIGS=1.402 ? ISPIN=2 ? GGA=91; VOSKOWN=1; EDIFF=0.1E-05; EDIFFG=-0.01 ? #IBRION=2 ? #NSW=50 ? #ISIF=2 (OR 3) ? NPAR=10

? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

POTCAR 文件直接在势库中拷贝 POSCAR文件 Al 4.05 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 4 Selective Dynamic Direct 1.0 0.0 0.0 T T T 0.5 0.5 0.0 T T T 0.5 0.0 0.5 T T T 0.0 0.5 0.5 T T T

? KPOINT 文件 ? Automatic generation ?0 ? Mohkorst Pack ? 21 21 21 ? 0.0 0.0 0.0

? 在进行能带和DOS计算时,ISMEAR 不能 使用N阶MP方法。因为MP方法在空轨道上 有负的占据,所以求得的能带和DOS是不 正确的。但是从其它地方看到---“提示: 在计算能带结构时,采ISMEAR = 0或1对 结果的影响非常小,可以认为是一样的。 但是不能采用ISMEAR = -5 或4。”ISMEAR到底多少?

? 计算能带:ICHARG = 11 导体的话,用ISMEAR=1; 半导体或绝缘体,用ISMEAR=0 。 ? 计算 DOS: ICHARG = 11 ISMEAR = -5 ? ? 计算的时候,金属可用0、1,非金属不要 大过0,体材料可用-4、-5(面的话就用 -1、0

? 设置完成后进行计算,计算完后,得到包含了态密度值的DOSCAR文 件,采用split_dos对态密度文件DOSCAR进行分割,得到总态密度 DOS0,各个原子的分波态密度DOS1,DOS2……。另外在运行 split_dos程序对DOSCAR文件分割时,要保证当前目录下有对应的 OUTCAR和POSCAR文件。 ? 分割后的DOS0,DOS1…等文件的能量值是以费米能级作为能量参 考零点。DOS0的第一列数据是能量值,单位为eV;第二列数据是总 态密度的值,单位 State/eV.unit cell;第三列数据是总态密度的积 分值,也就是电子数,单位为electrons。DOS1是第一个原子的分 波态密度值,其中的第一列数据是能量值,单位为eV;第二、三、四 列数据分别对应于s、p、d态的分波态密度值,单位为 State/eV.atom。其他的DOS文件与DOS1类似。

基本任务

? 计算电子态密度,能带,电荷密度 ? 优化晶体参数

? 内部自由度弛豫
? 结构弛豫

INCAR输入文件: 程序控制参数
System =diamond Si ISTART = 0 ENCUT = 150.0 NELM= 200 EDIFF = 1E-04 EDIFFG = -0.02 GGA=91 NPAR=4 NSW=100 IBRION = 2 ISIF=2 ISYM = 1 LWAVE = F LCHARG = F

POSCAR输入文件: 原胞中的原子位置

Diamond Si 3.9 0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 1 Direct 0.0 0.0 0.0
? ? ? 1 a1 ? a( j ? k ) 2 ? ? ? 1 a2 ? a (i ? k ) 2 ? ? ? 1 a3 ? a (i ? j ) 2

基矢的公因子 基矢a1 基矢a2 基矢a3 原胞中的原子个数 坐标系选为基矢构成的坐标系 基矢坐标系下原子的位置

KPOINTS输入文件: 控制K-点的选取方式

K-Points 0 Monkhorst Pack 11 11 11 0 0 0

POTCAR输入文件: 赝势文件
US Si 4.00000000000000000 parameters from PSCTR are: VRHFIN =Si: s2p2 LEXCH = CA EATOM = 115.7612 eV, 8.5082 Ry GGA = -1.4125 -1.4408 .0293 -.9884 eV TITEL = US Si LULTRA = T use ultrasoft PP ? IUNSCR = 1 unscreen: 0-lin 1-nonlin 2-no RPACOR = 1.580 partial core radius POMASS = 28.085; ZVAL = 4.000 mass and valenz RCORE = 2.480 outmost cutoff radius RWIGS = 2.480; RWIGS = 1.312 wigner-seitz radius (au A) ENMAX = 150.544; ENMIN = 112.908 eV EAUG = 241.945 …………

输出文件
OUTCAR CONTCAR CHGCAR CHG WAVECAR DOSCAR EIGENVAL OSZICAR LOCPOT PROOOUT

示例1: 用VASP求硅的电子态密度和能带
分如下几步:
(1). 生成4个输入文件: POSCAR POTCAR INCAR KPOINTS (2). 优化晶格参数,求出能量最低所对应的晶格参数

(3). 固定晶格参数, 求出能态密度(DOSCAR), 确定费米能量
(4). 修改KPOINTS和INCAR输入文件,固定电荷密度,做非自洽 计算,得到输出文件EIGENVAL (5). 提取数据,画图

(1). 生成4个输入文件: POSCAR POTCAR INCAR KPOINTS
Diamond Si 5.5 0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 2 Direct 0.0 0.0 0.0 0.25 0.25 0.25 System =diamond Si ISTART = 0 ENCUT = 150.0 NELM= 200 EDIFF = 1E-04 EDIFFG = -0.02

K-Points VASP提供 0 各种POTCAR Monkhorst Pack 21 21 21 0 0 0

NPAR=4 NSW=1 IBRION = 2 ISIF=2 ISYM = 1

(2). 优化晶格参数,求出能量最低所对应的晶格参数

运行VASP程序, 查看SUMMARY.fcc输出文件:

(3). 固定晶格参数, 求出能态密度(DOSCAR), 确定费米能量

(i)

找到平衡晶格常数后, 把该值写入到POSCAR文件中,并增加K点数 作一个离子步自洽计算(NSW = 0, IBRION = -1) .

(ii) 从DOSCAR输出文件中读出态密度和费米能级,费米 费米能级也可从OUTCAR中读出.

(4). 做非自洽计算, 求电子结构
? 修改INCAR文件: 将参数ICHARG设为 11 ? 修改KPOINTS输入文件 ? 运行VASP程序,从输出文件EIGENVAL中提出电子结构

画出电荷密度
? VASP输出电荷密度文件CHGCAR ? 采用免费程序LEV00处理数据文件CHGCAR www.cmmp.ucl.ac.uk/lev
1

0

(?)

-1

0 0.07 0.14 0.21 0.28 0.34 0.41 0.48 0.55

-2

-3 -3 -2 -1 0 1 2 3

(?)

示例2: 用VASP求Mg的电子态密度和能带
分如下几步:
(1). 生成4个输入文件: POSCAR POTCAR INCAR KPOINTS (2). 优化晶格参数,求出能量最低所对应的晶格参数

(3). 固定晶格参数, 求出能态密度(DOSCAR), 确定费米能量
(4). 修改KPOINTS和INCAR输入文件,固定电荷密度,做非自洽 计算,得到输出文件EIGENVAL (5). 提取数据,画图

(1). 生成4个输入文件: POSCAR POTCAR INCAR KPOINTS

Hcp-Mg 3.208 0.5 -0.866 0 0.5 0.866 0 0.0 0.0 1.6 2 Direct 0.0 0.0 0.0 0.66667 0.33333 0.5

VASP提供 各种POTCAR

K-Points 0 Monkhorst Pack 21 21 21 0 0 0

System =hcp Mg ISTART = 0 ENCUT = 150.0 NELM= 200 EDIFF = 1E-04 EDIFFG = -0.02 NPAR=4 NSW=1 IBRION = 2 ISIF=2 ISYM = 1

c/a

? 1? 3 ? a1 ? a( i ? j) 2 2 ? 1? 3 ? a2 ? a ( i ? j) 2 2 ? ? a3 ? ck

(2). 优化晶格参数,求出能量最低所对应的晶格参数

hcp结构晶体含有一个内部自由度, 晶格参数优化过程要比立方 结构费时

Mg: a=3.208, c/a=1.6

(3). 固定晶格参数, 求出能态密度(DOSCAR), 确定费米能量

(i)

找到平衡晶格常数后, 把该值写入到POSCAR文件中,并增加K点数 作一个离子步自洽计算(NSW = 0, IBRION = -1) .

(ii) 从DOSCAR输出文件中读出态密度和费米能级,费米 费米能级也可从OUTCAR中读出.

0.6

0.5

DOS

0.4

0.3

0.2

0.1 -6 -4 -2 0 2 4 6 8 10

Energy

(4). 做非自洽计算, 求电子结构

? 修改INCAR文件: 将参数ICHARG设为 11 ? 修改KPOINTS输入文件 ? 运行VASP程序,从输出文件EIGENVAL中提出电子结构

? 1? a1 ? a( i ? 2 ? 1? a2 ? a ( i ? 2 ? ? a3 ? ck

3? j) 2 3? j) 2

? 2? ? b1 ? (i ? a ? 2? ? b2 ? (i ? a ? 2? ? b3 ? k c

3 3 3 3

? j) ? j)

? ? ? ? ? 0b1 ? 0b2 ? 0b3 ? (0,0,0) 1 ? ? 1 1 K ? (b1 ? b2 ) ? ( , ,0) 3 ? 3 3 M ? 0.5b1 ? (0.5,0,0) ? A ? 0.5b3 ? (0,0,0.5) ? ? ? H ? 0.5b1 ? 0.5b2 ? 0.5b3 ? (0.5,0.5,0.5) ? ? ? L ? 0.5b1 ? 0b2 ? 0.5b3 ? (0.5,0,0.5)

KPOINTS文件:

k-points along high symmetry lines 100 ! 100 intersections Line mode rec 0 0 0 ! gama 0.3333 0.3333 0 ! K 0.3333 0.3333 0 !K 0.5 0.0 0.0 ! M 0.5 0.0 0.0 ! M 0 0 0 ! gama 0 0 0 0 0 ! gama 0.5 ! A

0 0 0.5 ! A 0.3333 0.3333 0.5 ! H
0.3333 0.3333 0.5 ! H 0.5 0.0 0.5 ! L 0.5 0.0 0.5 ! L 0 0 0.5 ! A

电荷密度
Be(1120)
Be(0001)

示例3: 用VASP求铅锌矿结构CoO的电子结构
设CoO呈铁磁性,故需做自旋极化计算
(1). 生成4个输入文件: POSCAR POTCAR INCAR KPOINTS (2). 优化晶格参数,求出能量最低所对应的晶格参数

(3). 固定晶格参数, 求出能态密度(DOSCAR), 确定费米能量
(4). 修改KPOINTS和INCAR输入文件,固定电荷密度,做非自洽 计算,得到输出文件EIGENVAL (5). 提取数据,画图

(1). 生成4个输入文件: POSCAR POTCAR INCAR KPOINTS
2.98 0.5 -0.866 0.0 0.5 0.866 0.0 0.0 2 2 Direct 0.0 0.0 0.0 0.66667 0.33333 0.5 0.66667 0.33333 0.1337 0.0 0.0 0.6337
? 1? 3 ? a1 ? a( i ? j) 2 2 ? 1? 3 ? a2 ? a ( i ? j) 2 2 ? ? a3 ? ck

0.0 1.735 VASP提供 各种POTCAR K-Points 0 Monkhorst Pack 21 21 21 0 0 0

System =hcp Mg ISTART = 0 ENCUT = 150.0 NELM= 200 EDIFF = 1E-04 EDIFFG = -0.02 ISPIN = 2 NPAR=4 NSW=1 IBRION = 2 ISIF=2 ISYM = 1

B A B A

(2). 优化晶格参数,求出能量最低所对应的晶格参数

wurtzite晶体含有两个内部自由度, 晶格参数优化过程要比立方 结构费时

CoO: a=2.98, c/a=1.735, u=0.367

8

8

4

4

E (eV)

0

0

-4

-4

-8

-8

-24

A

L

M

?

A

H

K

?

-24

A

L

M

?

A

H

K

?

DENSITY OF STATES

5

0

-5

-20

-10

0 ENERGY (eV)

10

20

固体材料表面的第一原理计算介绍

? VASP程序
? 构造超原胞 ? 计算表面性质 ? 应用

Building surfaces (1)asymmetric setup

(2)symmetric setup

unit cell coordinates are optimized Fixed layers (bulk)

vacuum

示例1: 用VASP求1*1Mg(0001)的表面性质
分如下几步:
(1). 生成4个输入文件: POSCAR POTCAR INCAR KPOINTS (2). 优化晶格参数,求出体Mg的晶格参数

(3). Mg(0001)的原子层数,构造超原胞的POSCAR
(4). 计算表面性质

(5). 提取数据,画图

POSCAR
Mg(0001): 3.208000 0.5 0.5 0.0 6 Direct 0.0 0.6666667 0.0 0.6666667 0.0 0.6666667

-0.8660254 0.8660254 0.0

0.0 0.0 10.22441

0.0 0.3333333 0.0 0.3333333 0.0 0.3333333

0.3048788 0.3829273 0.4609758 0.5390242 0.6170728 0.6951212

表面性质
Surface energy

??1 2 ( Esurf ? Natom ? Ebulk )
Geometry
getting relaxed structure from CONTCAR

di ? didea ?100% Relaxation of surface layers : didea
Heat of formation of overlayers of A on substrate B

H n( A) ?

Eslab( B) ? 2n( A) ? ( Eslab( B) ? 2nEbulk( A) ) 2

(Should use the same energy cutoff for each calculation)

Local Density of states INCAR: RWIGS = γ ? (works only for NPAR = 1 or serial version) LORBIT = 11 (only for PAW) ISMEAR = -5 (use tetrahedron for DOS calculations) NPAR = 1 Output file : DOSCAR (energy, s-dos, p-dos, d-dos for each atom) PROCAR (dos for each band and k-point)

Work function INCAR: LVTOT = .TRUE.
Output file : LOCPOT (same format as CHGCAR)
WRITE(IU,FORM) (((V(NX,NY,NZ),NX=1,NGX),NY=1,NGY),NZ=1,NGZ)

LOCPOT only contain electrostatic part of potential, if exchange correlation potential is to be included, change one line in main.F :
! comment out the following line to add exchange correlation ! INFO%LEXCHG=-1

1. Search ―E-fermi‖ in OUTCAR to get fermi-level 2. Analyze and plot data in LOCPOT

Z-axis(?)

Coulomb potential (eV)

Coulomb potential

(eV)

Fermi energy

Z-axis(?)

W(111)

Coulomb potential

(eV)

Coulomb potential

(eV)

Z-axis(?)

W(100)

Φ

Z-axis(?)

W(110 )

W(211 )

Mg(0001)的表面性质
10 8
Energy (eV)

6 4 2 0 -2 -4 -6
?

M

L

A

?

L

2
Energy (eV)

1 0 -1 -2
?

(a)

(c)

M

L

A

??

L ?

M

L

A

??

L

(b)
DOS

(d)

-6

-4

-2

0

2

4

6

8 -6

-4

-2

Energy (eV)

0 2 4 Energy (eV)

6

8

示例2: 原子氢在Mg(0001)表面的吸附性质
POSCAR
Mg(0001)+H: 3.208000 0.5 0.5 0.0 6 2 Direct 0.0 0.6666667 0.0 0.6666667 0.0 0.6666667 0.3333333

-0.8660254 0.8660254 0.0

0.0 0.0 10.22441

0.0 0.3333333 0.0 0.3333333 0.0 0.3333333 0.6666667

0.3048788 0.3829273 0.4609758 0.5390242 0.6170728 0.6951212 0.721

示例3: 原子氢在

表面的吸附性质

最佳吸附位置

示例4: Mg在Si(111)表面的吸附性质

形成能 (formation energy)

Si(111)/Mg系统的电子结构性质

Clean Si(111)

0.25ML

0.5ML

1ML

3
-0.040 -0.022

2

-0.005 0.013

1

Mg

0.030 0.048 0.065 0.083

(?)

0

0.100

-1

Si

Si

-2

-3 -3 -2 -1 0 1 2 3

3

(?)
-0.040 -0.022

2

-0.005 0.013 0.030

1

Mg

0.048 0.065 0.083

(?)

0

0.100

-1

Si

Si

-2

-3 -3 -2 -1 0 1 2 3

(?)
3
-0.017 -0.010

2

-0.003 0.004 0.011

1

Mg

0.018 0.026 0.033

(?)

0

0.040

Si
-1

Si

-2

-3 -3 -2 -1 0 1 2 3

(?)

示例5: ZnO在MgO(111)表面上的极向生长

Formation energy:

吸附能

O-O 反相层

More and More by VASP

晶格动力学: 声子谱,声子态密度
结合phonon程序, 冻结声子法

hcp-Be

高压相变

弹性系数

熔化曲线,第一原理分子动力学
金属钼的熔化曲线

扩散系数和粘性系数的计算结果
V 15.48 10.98 T(K) 3500 3869 5500 6500 7000 7000 7500 8000 D(109m2/s) 1.460 5.116 0.189 0.192 6.669 0.885 1.324 ?(mPa?s) 16.35 5.158 259.2 222.6 8.028 62.751 35.694

9.84

200GPa附近不同温度的径向分布函数 (钼)

Linux与并行计算

主要内容
? 并行计算简介 ? Linux与HPC ? 搭建Linux Cluster ? 小结

? 并行计算(Parallel Computing) ? 高性能计算(High Performance Computing) ? 并行计算,就是在并行机上所做的计算 ? 理论、实验与计算,被认为是人类认识自然 的三大支柱。
?降低单个问题求解的时间 ?提高问题求解精度 ?增加问题求解规模 ?无法替代的模拟计算

并行计算--高性能计算

为什么要做并行计算? --应用需求

并行计算机
? 并行计算机就是由多个处理单元组成的计 算机系统,这些处理单元相互通信和协作 以快速、高效求解大型复杂问题。
处理单元有多少 处理单元的功能有多强 处理单元之间怎样连接 处理单元的数据如何传递 各处理单元如何相互协作 并行程序如何编写

并行计算机的发展
? (1)并行机的萌芽阶段 (1964-1975) ? (2)向量机/ SMP的发展 和鼎盛阶段(19761990) ? (3)MPP出现和蓬勃发 展阶段(1990-1995) ? (4)Cluster出现,并成 为并行计算机的主流 (1995年后)

并行计算机结构模型

各种架构对比
? 共享存储体系结构
– SMP瓶颈在访存带宽,限制了可扩展性(性能会严 重下降),通过增大Cache来缓解,又有数据一致 性的问题。 – DSM具有较好的可扩展性,且具有与SMP相同的编 程模式。共享存储的机器可以运行各种编程模式的 程序,包括串行程序、共享存储并行程序和消息传 递并行程序。 – 共享存储的机器首先可以看作是一个内存扩充了的 串行机器,可以运行数据量极大的串行程序,如大 数据量的数据库应用。利用多个处理器的并行处理 性能,自动并行(靠编译器)作用有限,需要其他 编程方式的应用,如OpenMP。

各种架构对比(续)
? 分布式存储体系结构
– 无共享的机器,MPP和Cluster,具有最好的可扩展性, 采用消息传递编程,可编程性不如OpenMP,也可以 通过DSM软件的方式来模拟实现共享存储(其实还是 通过消息传递,只不过对上层通明),性能受到影响, 但可以运行OpenMP应用程序。 – MPP/PVP在构造大规模系统,应用饱和性能方面具有 优势,资金充足的依然会选择; – Cluster由于无可比拟的性价比优势占据主流位置。
? ? ? ? 更高的计算能力 系统的高可用性 良好的可扩展性 更高的性价比

主要内容
? 并行计算简介 ? Linux与HPC ? 搭建Linux Cluster ? 小结

www.top500.org

Linpack
? Linpack是国际上流行的用于测试高性能计算机系 统浮点性能的benchmark ? 通过对高性能计算机求解线性代数方程组能力的测 试,评价高性能计算机的浮点性能 ? HPL(High Performance Linpack)
– 针对大规模的并行计算机系统的测试 – HPL版Linpack一般用于并行超级计算机 – 用户可以选择矩阵的大小(问题规模)、使用各种优化方法来执行测 试程序,寻求最佳的测试结果 – 计算量(2/3 * N^3 – 2*N^2) / 计算时间 – http://www.netlib.org/benchmark/hpl/

Linux在HPC领域的优势
硬件要求 低 软件成本 低 维护成本 低 并行代码 多 开源工具 多 可定制性 好

Linux在HPC领域的劣势

发行版本太多 推广普及不够

主要内容
? 并行计算简介 ? Linux与HPC ? 搭建Linux Cluster ? 小结

集群系统体系结构

采用Linux来搭建集群系统
? 基本流程
– 对头结点/管理结点
? 操作系统安装 ? 集群服务配置

– 对每个计算结点
? 操作系统安装 ? 集群服务配置

– 集群系统联调

Linux的七种武器
? Linux的七种武器(RedHat)
– dd:操作系统克隆 – rsh:远程登录、执行 – mpich:并行编程和运行环境 – nfs:网络共享文件系统 – ntp:集群时间同步 – nis:集群帐号同步 – openpbs:作业调度系统

七种武器之一:dd
? 操作系统克隆
– echo scsi add-single-device 0 0 1 0 > /proc/scsi/scsi
? Host: scsi0 Channel: 00 Id: 01 Lun: 00

– dd if=/dev/sda of=/dev/sdb bs=4096
? if:系统源盘 ? of:系统镜像盘 ? bs:块大小

– 73GB SCSI硬盘 17分钟 – 要求源盘与镜像盘型号完全一致

七种武器之二:rsh
? 远程登录、执行
– 实现节点之间的无密码切换 – 打开服务rsh、rexec、rlogin
? chkconfig rsh on ? chkconfig rexec on ? chkconfig rlogin on ? /etc/init.d/xinetd restart

– 将所有节点名加入/etc/hosts.equiv

– 对于root用户
? echo rsh >> /etc/securetty ? echo rexec >> /etc/securetty ? echo rlogin >> /etc/securetty ? cp /etc/hosts.equiv /root/.rhosts

– ssh:更安全的remote shell

七种武器之三:mpich
? 并行编程和运行环境
– 并行编程标准
? 多线程库标准
– Win32 API. – POSIX threads.

? 编译制导标准
– OpenMP – 可移植共享存储并行编程标准.

? 消息传递库标准–最重要的并行程序设计方式
– MPI – PVM

– mpich的安装
? ? ? ? ? tar -xzvf mpich-1.2.6.tar.gz cd mpich-1.2.6 ./configure --prefix=/usr/local/mpich-1.2.6 make make install

– 环境变量设置
? vi /etc/profile 加入下面的行
export MPI_PATH=/usr/local/mpich-1.2.6 export PATH=$MPI_PATH/bin:$PATH export MANPATH=$MPI_PATH/man:$MANPATH

? source /etc/profile

– MPICH的使用
? 最简单MPI程序
#include <stdio.h> #include "mpi.h“ int main( int argc,char *argv[] ) { MPI_Init( &argc, &argv ); printf( "Hello, world!\n" ); MPI_Finalize(); return 0; }

? MPI程序的编译
– mpicc -o hello hello.c

? MPI程序的运行
– mpirun -np 4 hello

– mpirun -np 8 -machinefile ma hello

七种武器之四:nfs
? 网络共享文件系统
– 头节点
? 1.打开nfs服务
– chkconfig nfs on – /etc/init.d/nfs start

? 2.编辑/etc/exports 假设将/public共享给其它节 点 /public *(rw,no_root_squash,async)
– 不限制用户、可读写、信任客户端、先写入内存

? 3.执行exportfs -a,共享/public目录
/public <world> – 执行 exportfs ,可以看到

– 计算节点
? 1.建/public目录
– mkdir /public

? 2.手动挂载nfs
– mount server:/public /public

? 3.看mount 情况

? 4.开机自动挂载nfs
– echo mount node38:/public /public >> /etc/rc.local

七种武器之五:ntp
? 集群时间同步
– 头节点
? 1.设置系统时间
– date -s 20071215 – date -s 15:35 – hwclock --systohc

? 2. 修改/etc/ntp.conf文件
– vi /etc/ntp.conf # -- CLIENT NETWORK ------restrict 192.168.0.0 mask 255.255.255.0 notrust nomodify notrap

? 3. 打开ntp服务
– chkconfig ntpd on – /etc/init.d/ntpd start

– 计算节点
? 1.与头节点手动同步 – ntpdate node1
– hwclock --systohc

? 2. 修改/etc/ntp.conf文件
– vi /etc/ntp.conf # --- OUR TIMESERVERS ----restrict 192.168.0.1 mask 255.255.255.255 nomodify notrap noquery server 192.168.0.1

? 3. 打开ntp服务
– chkconfig ntpd on – /etc/init.d/ntpd start

七种武器之六:nis
? 集群帐号同步
– 头节点
? 1.设定NIS域名
– nisdomainname dawning – echo NISDOMAIN=dawning >> /etc/sysconfig/network

? 2.打开NIS服务
– chkconfig ypserv on – chkconfig yppasswdd on – /etc/init.d/ypserv start – /etc/init.d/yppasswdd start

? 3.初始化NIS数据库
– /var/lib64/yp/ypinit -m – /etc/init.d/ypserv restart

? 4.nfs共享头节点/home目录
– vi /etc/exports
/home *(rw,no_root_squash,async)

– exportfs –a

– 计算节点
? 1.设定NIS域名
– nisdomainname dawning – echo NISDOMAIN=dawning >> /etc/sysconfig/network

? 2.打开NIS服务
– – – – chkconfig ypbind on chkconfig yppasswdd on /etc/init.d/ypbind start /etc/init.d/yppasswdd start

? 3.使用头节点NIS数据库
– vi /etc/nsswitch passwd: shadow: group: files files files passwd: shadow: group: files nis files nis files nis

? 4.挂载头节点/home目录
– mount node1:/home /home – echo mount node1:/home /home >> /etc/rc.local

? 5.测试NIS
– yptest

– 增加/删除/修改用户
? 在头节点上用useradd/userdel/usermod ? 重新生成NIS数据库
– /var/lib64/yp/ypinit -m – /etc/init.d/ypserv restart

– 修改用户密码
? 不能使用passwd,而要使用yppasswd ? 修改 yp 密码会改变系统密码﹐修改系统密码却不 改变 yp 密码

七种武器之七:openpbs
? 作业调度系统
– 系统组成

– 安装
? 解压源文件包
– tar -zxvf openpbs-2.3.16.tar.gz – cd OpenPBS_2.3.16

? 编译设置
– ./configure --disable-gui –enable-docs --set-server_home=/var/spool/pbs --x-libraries=/usr/X11R6/lib64

? 编译安装
– make – make install

– 配置

? 系统配置文件/etc/pbs.conf
头节点 #!/bin/sh pbs_home=/var/spool/pbs pbs_exec=/usr/local start_server=1 start_sched=1 start_mom=1

计算节点 #!/bin/sh pbs_home=/var/spool/pbs pbs_exec=/usr/local start_server=0 start_sched=0 start_mom=1

? 系统启动脚本
– – – –

/etc/init.d/openpbs Server的系统启动脚本 /etc/init.d/pbs_server Scheduler系统启动脚本 /etc/init.d/pbs_sched Mom系统启动脚本 /etc/init.d/pbs_mom

– 头节点
? 初始化server: (第一次运行或者重新配置)
– /usr/local/sbin/pbs_server –t create

? Server配置目录
– /var/spool/pbs/server_priv/

? 节点属性声明: node2 R220A np=2 /var/spool/pbs/server_priv/nodes
node3 R220A np=2 node4 dualcore np=4 node5 dualcore np=4 node6 R4280A np=4 node7 R4280A np=4

? 队列设置
– 导入server配置文件:qmgr < queue.conf – 输出配置文件:qmgr –c “print server” > queue.conf – 配置文件例子:
create queue default set queue default queue_type = execution set queue default max_running = 20 set queue default enabled = True set queue default started = True set server scheduling = True set server max_user_run = 20 set server default_queue = default set server query_other_jobs = True

– 计算节点
? mom配置目录: /var/spool/pbs/mom_priv/ ? mom配置文件: /var/spool/pbs/mom_priv/config # MOM server configuration file
# if more than one value, separate it by comma. ## rule is defined by the name $ideal_load 1.5 $max_load 2 ## host allowed to connect to Mom server on unprivileged port $restricted *. ## log event : # 0x1ff log all events + debug events # 0x0ff just all events $logevent 0x0ff ## host allowed to connect to mom server on privileged port $clienthost node1 ## alarm if the script hang or take very long time to execute $prologalarm 30

– 使用
? qsub job.pbs
#PBS –N vasp.Hg #PBS –l nodes=8:ppn=2 echo "This jobs is "$PBS_JOBID@$PBS_QUEUE cd $PBS_O_WORKDIR mpirun -np 16 -machinefile $PBS_NODEFILE ./vasp

? qstat

NPACI Rocks
? National Partnership for Advanced Computational Infrastructure (NPACI) 的Rocks软件,是和Linux操作系统捆绑的 集成化的中间件系统,安装和管理都十分 方便。http://www.rocksclusters.org/ ? NPACI Rocks软件被广泛地应用于学术和 政府机构,其中包括西北大学、西北太平 洋国家实验室、斯坦福大学等。

Rocks Basic Approach
? Install a frontend
1. Insert Rocks Base CD 2. Insert Roll CDs (optional components) 3. Answer 7 screens of configuration data 4. Drink coffee (takes about 30 minutes to install)

?

Install compute nodes:
1. Login to frontend 2. Execute insert-ethers 3. Boot compute node with Rocks Base CD (or PXE) 4. Insert-ethers discovers nodes 5. Goto step 3 Optional Rolls
– – – – – – – – – Condor Grid (based on NMI R4) Intel (compilers) Java SCE (developed in Thailand) Sun Grid Engine PBS (developed in Norway) Area51 (security monitoring tools) Many Others …

? ?

Add user accounts Start computing

Rocks的七种武器
? Rocks的七种武器
– insert-ethers:操作系统快速部署 – cluster-fork:集群并行命令 – mpich/openmpi:并行编程和运行环境 – pvfs:并行文件系统 – ganglia:集群状态监控 – 411:集群帐号同步 – sge/torque:作业调度系统

主要内容
? 并行计算简介 ? Linux与HPC ? 搭建Linux Cluster ? 小结

小结
? 集群是高性能计算机的发展方向 ? Linux在高性能计算机中大量采用 ? 配置Linux Cluster非常方便 ? 采用NPACI Rocks可以大大减少工作量

13.1 Numerical Methods (数值方法) History of FEM (有限元发展历史)
?In 1943,Courant proposed a approximative method in his mathematic paper which was similar to FEM (1943年Courant在他的应用数学论文中首先提出了 一个与有限单元法类似的近似方法。) ?In 1960,Clough proposed Finite Element Method (FEM) ( 1960年,Clough提出了有限元法)

Clough

13.2 Basics of the Final Element Method (有限元法基础)
Element and Node (单元和节点)
node (节点)

element(单元)

FEM discretizes the domain under study by dividing the region into subdomains called elements. (有限元方法将研究区域离散成许多子域,这些子域叫单元)

识别无效的结果
? ? ? ? 分析的对象的一些行为 计算出的几何项 求解的自由度及应力 反作用力或节点力

1.分析的对象的一些基本的行为:
? ? ? ? ? 重力方向总是竖直向下的 离心力总是沿径向向外的 没有一种材料能抵抗 1,000,000 psi 的应力 轴对称的物体几乎没有为零的 环向应力 弯曲载荷造成的应力使一侧受压,另一侧受拉

如果只有一个载荷施加在结构上,检验结果比较容易. 如果有多个载荷,可单独施加一个或几个载荷分别 检验,然后施加所有载荷检验分析结果.

236

2.计算出的几何项:
在输出窗口中输出的质量特性,可能会揭示在几何 模型、材料属性(密度)或实常数方面存在的错误.

3.检验求解的自由度及应力:
? 确认施加在模型上的载荷环境是合理的. ? 确认模型的运动行为与预期的相符 - 无刚体平动、 无刚体转动、无裂缝等. ? 确认位移和应力的分布与期望的相符,或者利用物 理学或数学可以解释.
237

4.反作用力或节点力
模型所有的反作用力应该与施加的点力、压力和惯性力 平衡.
在所有约束节点的竖 直方向的反作用力...

…必须与施加的竖直方 向的载荷平衡 在所有约束节点水平方向的反 作用力必须与水平方向的载荷 平衡. 所有约束节点的反作用力矩必 须与施加的载荷平衡. 注意包含在约束方程中自由度 的反力,不包括由这个约束方 238 程传递的力.

反作用力和节点力 (续)
在任意选取的单元字集中的节点力,应与作用在结构 此部分的已知载荷向平衡,除非节点的符号约定与自 由体图上所示的相反.

未选择的单元上的竖 直方向的节点总力...

…必须与被选择的 单元上施加的竖直 方向的载荷平衡

注意包含在约束方程中 自由度的反力,不包括 由这个约束方程传递的 力. 239

ANSYS网格划分精度估算
? 网格误差估算 ? 局部细化 ? P方法&举例

ANSYS网格误差估计
ANSYS通用后处理包含网格离散误差估计.
误差估计是依据沿单元内边界的应力或热流的不连续性,是平均 与未平均节点应力间的差值.

?avg = 1100

? = 1000 Elem 1 ? = 1100

? = 1200 Elem 2 ? = 1300

节点的 ?s 是积分点 的外插)
(

?avg = 1200

241

ANSYS网格误差估计
误差估计作用条件:
? 线性静力结构分析及线性稳态热分析 ? 大多数 2-D 或 3-D 实体或壳单元 ? PowerGraphics off

误差信息:
? ? ? ? 能量百分比误差 sepc 单元应力偏差 sdsg 单元能量偏差 serr 应力上、下限 smnb smxb
242

能量百分比误差
能量百分比误差是对所选择的单元 的位移、应力、温度或热流密度的 粗略估计. 它可以用于比较承受相 似载荷的相似结构的相似模型.

这个值的通常应该在10%以下. 如 果不选择其他单元,而只选择在节 点上施加点载荷或应力集中处的单 元,误差值有时会达到50%或以上.
SEPC ~ 2 %

PowerGraphic off Main menu > general postproc > plot results > deformed shape 选 :Def+undefedge
243

应力偏差
要检验某个位置的网格离散应 力误差,可以列出或绘制应力 偏差. 某一个单元的应力偏差是此单 元上全部节点的六个应力分量 值与此节点的平均应力值之差 的最大值. 应力偏差:

节点n的应力矢量:

??? ?? ?? ?? ?? ?
i n a n i n

所关心位置上的应力偏差值~450 psi (30,000 psi 应力的1.5%)

N en 察看应力偏差:Plot Results > Element Solu > Error Estimation > Stress deviation (SDSG)

?? ??
a n

i ? ? n i ?1

n Ne

244

举例

? 平均应力为4421 (nodal solution) ? 应力偏差为689.598 ? 误差=689.598/4421=15.53%(局部细 化)

能量误差
每个单元的另一种误差值是能量误差. 它与单元上节点应力差值 有关的, 用于计算选择的单元的能量百分比误差.
ei ? 1 T ?1 ? ? ??? ?d (vol) ? ? ? ? D ? 2 vol

其中:ei ? 单元i的能量误差

?D? — 单元的应力 ? 应变矩阵 ??? ? — 应力误差矢量
整个模型的能量误差: e ? ? ei Nr单元数
i ?1 nr

vol — 单元体积

察看能量误差:Plot Results > Element Solu > Error Estimation > Energy error (ENER).

246

应力上下限
应力上下限可以确定由于网格离散误 差对模型的应力最大值的影响.

显示或列出的应力上下限包括: ? 估计的上限 - SMXB ? 估计的下限 - SMNB 应力上下限限并不是估计实际的最 高或最小应力。它定义了一个确信 范围。 如果没有其他的确凿的验证 ,就不能认为实际的最大应力低于 SMXB.

?

mnb j

? min( ?

a jm

? ?? n )

X stress SMAX ~ 32,750 psi SMXB ~ 33,200 psi (difference ~ 450 psi ~ 1.5 %) 247

例如:SMX=32750是节点解的实际值 SMXB=33200是估计的上限

a ? mxb ? max( ? j jm ? ?? n )

局部的细化
采用plane42单元网格局部细化与未细化
能量百分比误差 局部细化
Displacement DMX=0.88E-03 SEPC=14.442 应力偏差 Element Solution(SDSG) SMN=63.453 SMX=426.86

未细化
DMX=0.803E-03 SDSG SMN=64.528 SMX=689.589

能量误差估计

Element Solution(SERR) SMN=0.365E-03 SMX=0.600595
Nodal solution(SEQV) SMN=725.21

SERR SMN=0.005173 SMX=0.38503
SEQV SMN=773.769 SMNB=708.94

应力上下限

P方法及p单元的应用
P 单元的位移形函数 u=a1+a2x+a3y+a4x2+a5xy+a6y2

v=a7+a8x+a9y+ a10x2+a11xy+a12y2
P方法的优点:
如果使用 p-方法 进行结构分析,可以依靠p单元自动调整单元多项式阶数(28),达到收敛到设定的精度. 对这种方法的相信程度,与使用经验有关.

P方法应用控制:
P方法用于线弹性结构分析—实体和壳 体。 P单元由以下5种单元: ?2-D Quadrilateral (Plane145) ?2-D Triangle (Plane146) ?3-D Brick (Solid 147) ?3-D Tetrehedron (Solid 148) ?3-D Shell (Solid 150)

规定 0.1% 局部应力差,使用p方法计算的最 249 大X方向应力约为 34,700 psi (比普通h方法高出大约 5% )

P方法进行静力分析的步骤
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1.选择P方法作业 GUI:Main Menu > Preference > P-Method 定义一个P单元,P方法被激活。 2.建模 建模过程与H-单元分析相同,单元类型必须用P单元 (a)指定P单元 水平 定义局部P-水平等级 定义P单元时用KeyOpt选项定义 定义整体p-水平等级 命令: PPRANGE , START, MAX GUI: Main Menu > Solution > P-Method > Set P Range (b)定义几何模型 应用实体建模 (c) 用P单元分网。 自适应网格对P方法是无效的 3.施加载荷、求解 应用实体模型加载,而不是有限元模型 求解:推荐采用条件共轭梯度法(PCG),但PCG对于壳体P单元无效 4.后处理 察看结果

举例: platep.dat
20 in pres=-100lb/in 2

R=5 in

10 in

? ? ? ?

E=30e6 lb/in2 V=0.29 Thick=0.25in 在节点(0,5,0)处的收敛标准设为1%

高级网格划分技术
? 延伸网格划分 ? 映射网格划分 ? 层状网格划分

延伸网格划分 & 举例
? 将一个二维网格延伸生成一个三维网格;三维网格生成后 去掉二维网格 ? 步骤: ? 1.先生成横截面 ? 2.指定网格密度并对面进行网格划分 ? 3.拖拉面网格生成体网格

?
?

指定单元属性

拖拉,完成体网格划分。

? 4.释放已选的平面单元

举例:飞机模型机翼
y

弹性模量
x

Ex=38E03 psi 泊松比:0.3 密度:
10

斜度=0.25

z

D=1.033e-3 slugs/in3

2

机翼沿着长度方向轮廓一致,且它的横截面由直线和样条曲线 定义。机翼的一端固定在机体上,另一端为悬空的自由端。
采样点:A(0,0,0) B(2,0,0) C(2.3,0.2,0) D(1.9,0.45,0) E(1,0.25,0)

延伸网格划分:作业

截面宽度:10mm 手柄长度: 20cm 导角半径: 1cm

截面形状:正六变形 杆长 : 7.5cm

弹性模量: 2.07E11pa

映射网格划分
? 有两种主要的网格划分方法: 自由划分和映射划分. ? 自由划分 – 无单元形状限制. – 网格无固定的模式. – 适用于复杂形状的面和体.
? 映射划分

– 面的单元形状限制为四边形,体的单元限制为六面 体 (方块). – 通常有规则的形式,单元明显成行. – 仅适用于 “规则的” 面和体, 如 矩形和方块.

映射网格划分
网格划分的优缺点:
自由网格 + 易于生成; 不须将复杂形状的 体分解为规则形状的体. – 体单元仅包含四面体网格, 致 使单元数量较多. – 仅高阶 (10-节点) 四面体单元 较满意, 因此DOF(自由度)数 目可能很多.
+ + – –

映射网格 通常包含较少的单元数量. 低阶单元也可能得到满意的结 果,因此DOF(自由度)数目较少. 面和体必须形状 “规则”, 划 分的网格必须满足一定的准则. 难于实现, 尤其是对形状复杂 的体.

...映射网格划分
自由网格
? 自由网格是面和体网格划分时的缺省设置. ? 生成自由网格比较容易: – 导出 MeshTool 工具, 划分方式设为自由 划分. – 推荐使用智能网格划分 进行自由网格划分, 激活它并指定一个尺寸级别. 存储数据库. – 按 Mesh 按钮开始划分网格. ? 按拾取器中 [Pick All] 选择所有实 体 (推荐). – 或使用命令 VMESH,ALL 或 AMESH,ALL.

映射网格划分&举例
映射网格划分
? 由于面和体必须满足一定的要求,生成映射网格不如生成自由网格容 易: – 面必须包含 3 或 4 条线 (三角形或四边形). – 体必须包含4, 5, 或 6 个面 (四面体, 三棱柱, 或六面体). – 对边的单元分割必须匹配. ? 对三角形面或四面体, 单元分割数必须为偶数.

...映射网格划分
? 因此 ,映射网格划分包含以下三个步骤: – 保证 “规则的”形状, 即, 面有 3 或4 条边, 或 体有 4, 5, 或 6 个面. – 指定尺寸和形状控制 – 生成网格0

...映射网格划分
1.保证规则的形状
? 在许多情况下, 模型的几何形状上有多于4条边的面,有多于6个面 的体. 为了将它们转换成规则的形状, 您可能进行如下的一项或 两项操作: – 把面 (或体) 切割成小的, 简单的形状. – 连接两条或多条线 (或面) 以减少总的边数.

...映射网格划分
? 切割 (divide)可以通过布尔减运算实现.
– 您可以使用工作平面, 一个面, 或一条线 作为切割工具. – 有时, 生成一条新的线或面会比移动或定向工作平面到正确的 方向容易得多.

...映射网格划分
? 连接 操作是生成一条新线 (为网格划分) , 它通过连接两条或
多条线以减少构成面的线数. – 使用 LCCAT 命令或 Preprocessor > -MeshingConcatenate > Lines, 然后拾取须连接的线. – 对面进行连接, 使用 ACCAT 命令或Preprocessor > Meshing- Concatenate > Areas – 若两条线或两个面 相切交汇可考虑用加 (布尔) 运算

连接这两条线使其成为一个 由4条边构成的面Concatenate

举例:一个由六条线围成的面
L1 and L2 are added. New Line L# L4 and L5 are concatenated.

产生四条线围成的面,适于网格划分

New concatenated line.

连接边的单元数为8条

原始边的单元数为4条

...映射网格划分
? 您也可以简单地通过一个面上的3个或4个角点 暗示 一个连接. 此时, ANSYS 内在地 生成一个 连接. – 在MeshTool中选择Quad shape 和 Map 网格. – 将 3/4 sided 变为 Pick corners. – 按 Mesh 键, 拾取面, 然后拾取 3 或 4 角点 形成一规则的形状.

...映射网格划分
? 使用连接时注意:
– 它仅仅是一个网格划分操作,因而应为网格划分前的最后一步, 在所有的实体建模之后. 这是因为,经连接操作得到的实体 不能在后续的实体建模操作中使用. – 可以通过删除产生的线或面 “undo(取消)” 一个连接. – 连接面 (为在体上映射网格) 通常比较复杂,因为您也应该连 接一些线. 只有在对相邻的两个4边形面作连接时其中的线会 自动连接.

...映射网格划分
指定尺寸和形状控制 ? 这是映射网格划分3个步骤中的第2步. ? 选择单元形状非常简单. 在 MeshTool中, 对面的网格划分选择 Quad,对体的网格 划分选择 Hex, 点击 Map. ? 其中通常采用的尺寸控制和级别如下:
– 线尺寸 [LESIZE] 级别较高. – 若指定了总体单元尺寸, 它将用于 “未给定尺 寸的” 线. – 缺省的单元尺寸 [DESIZE]仅在未指定 ESIZE时用于 “未给定尺寸的” 线上. – (智能网格划分 无效.)

...映射网格划分
? 若您指定线的分割数, 切记: – 对边的分割数必须匹配, 但您只须指定一边的分割数. 映射网 格划分器 将把分割数自动传送到它的对边. – 如果模型中有连接(Concatenate)线, 只能在原始(输入)线 上指定分割数,而不能在合成线上指定分割数.
每条初始线上指定6份分割.

此线上将自动使用12 份分 割 (合成线的对边).

其它两条线上会采用几 份分 割 呢? (后面的演示将会回 答这一问题.)

...映射网格划分
生成映射网格 ? 只要保证了规则的形状 并指定了合适的份数, 生成网格将非常简 单. 只须按MeshTool中的 Mesh 键, 然后按拾取器中的 [Pick All] 或选择需要的实体即可.

映射网格划分举例:轮
说明 ? 用自由及映射网格对轮模型进行混合的网格划分.


1. 按指定的工作目录,以 “wheelb-3d”为作业名,进入ANSYS. 或清除 ANSYS 数据库,改换作业名为 “wheelb-3d”: – Utility Menu > File > Clear & Start New ... – Utility Menu > File > Change Jobname ... 2. 恢复 “wheelb.db1” 数据库文件: – Utility Menu > File > Resume from … ? 选择 “wheelb.db” 数据库文件, 然后选择 [OK] – 或用命令: RESUME,wheelb,db1 3. 进入前处理器,用工作平面切分体 : – Main Menu > Preprocessor > -Modeling- Operate > Booleans- Divide > Volu by WrkPlane + ? 拾取[Pick All] – Utility Menu > Plot > Volumes – 或用命令: /PREP7 VSBW,1 VPLOT

4. 平移工作平面到19号关键点: – Utility Menu > WorkPlane > Offset WP to > Keypoints + ? 选择如图所示的19号关键点, 然后选择 [OK] – 或用命令: KWPAVE,19

5. 以工作平面切分体: – Main Menu > Preprocessor > -Modeling- Operate > Booleans- Divide > Volu by WrkPlane + ? 拾取[Pick All] – Utility Menu > Plot > Volumes – 或用命令: VSBW,4 VPLOT

6. 关闭工作平面,设置总体单元尺寸为 0.25: – Utility Menu > WorkPlane > Display Working Plane – Main Menu > Preprocessor > MeshTool … ? 设置大小控制为Global,按 [Set] ? 设置SIZE = 0.25 ? 按[OK] – 或用命令: WPSTYLE ESIZE,0.25 7. 用SOLID45单元,对四个外部的体进行映射网格划分 (TYPE 1): – Main Menu > Preprocessor > MeshTool … ? 在Shape下选择 “Hex” 和 “Mapped” : ? 按[Mesh] ? 拾取四个外部的体 (体号 1, 2, 3, 和 5) ? 按[OK] – 或用命令: MSHAPE,0,3D MSHKEY,1 VMESH,1,3,1

VMESH,5

8. 用SOLID95单元,对内部的体进行自由网格划分 (TYPE 2): – Main Menu > Preprocessor > MeshTool … ? 单击单元属性(Element Attributes)下的 [Set] : ? TYPE = “2 SOLID95”, 然后选择 [OK] ? 设置大小控制为Global,按 [Set] ? 设置SIZE = 0.2 ? 按[OK] ? 在Shape 下选择 “Tet”和“Free” : ? 按[Mesh] ? 拾取内部的体 (体号 6) ? 按[OK] – 或用命令: TYPE,2 ESIZE,0.2 MSHAPE,1,3D MSHKEY,0 VMESH,6

9.


– 10. –

– –

将 SOLID95单元转变为 SOLID92单元: Main Menu > Preprocessor > -Meshing - Modify Mesh > Change Tets ... ? 按[OK] 或用命令: TCHG,95,92 选择并画出 SOLID95 四面体单元: Utility Menu > Select > Entities ... ? 选择 “Elements”, “By Attributes”, “Elem type num” ? 设置 Min,Max,Inc = 2 ? 按[OK] Utility Menu > Plot > Elements 或用命令: ESEL,S,TYPE,,2 EPLOT

11. – – –



11. –



选择 “全部实体”并保存数据库 : Utility Menu > Select > Everything Utility Menu > Plot > Elements Utility Menu > File > Save as … ? 输入数据库文件名“wheelb-3d-mesh.db”, 然后选择 [OK] 或用命令: ALLSEL,ALL EPLOT SAVE,wheelb-3d-mesh,db 退出ANSYS: 在工具条中选择 “QUIT” ? 选择 “Quit - No Save!” ? 按[OK] 或用命令: FINISH /EXIT,NOSAVE

层状网格划分
? 适用于2D情况,生成线性过渡的 自由网格 ? 平行于边线方向的单元尺寸相当 ? 垂直于边线方向的单元尺寸和数 目急剧变化 ? 当分析要求边界单元高精度时, 层状网格很有用

层状网格划分
GUI:Main Menu: Preprocessor > MeshTool > Layer > Set button 指定:线上的单元尺寸,线上两端单元的比 率和内部网格层的厚度。 ?线间距比率(space),对层状划分一般取1.0

?内部网格层厚度(layer1)线上单元尺寸系
数: size factor=2 沿线生成两行尺寸均匀的单元

?外部网格层厚度(layer2)这层的单元尺寸
会从layer1缓慢增加到总体单元尺寸,layer2的 厚度可以用一个网格过渡系数 如:Transition factor=2 生成大约等于前面垂 直于线网格2倍尺寸的单元

在完成此章的学习之后,给出一个已经划分好网 格的模型的数据库文件,我们应该能够使用耦合 或约束方程来建立节点自由度之间的联系

第一讲 耦合 定义耦合设置 说明耦合的三种普遍应用. 采用3种不同的方法建立耦合关系.

第二讲 约束方程 定义“约束方程” 说明约束方程的四种普遍应用 采用四种不同的方法生成约束方程.

2015-3-16

281

耦合设置
?

耦合是使一组节点具有相同的自由度值.
– 除了自由度值是由求解器计算而非用户指定外,与约

? ?

束相类似。 – 例如:如果节点1和节点2在UX方向上耦合, 求解器将 计算节点1的UX值并简单地把该值赋值给节点2的UX。 一个耦合设置是一组被约束在一起,有着同一方向的节 点 (即一个自由度)。 一个模型中可以定义多个耦合,但一个耦合中只能包含 一个方向的自由度。

282

耦合的三种一般应用
1. 施加对称性条件: ? 耦合自由度常被用来实施移动或循环对称条件. ? 考虑在均匀轴向压力下的空心长圆柱体,此3-D 结构可用下面右图所示的2-D轴对称模型表示.
y y

11 x 1

12 2

13 3

14 4

15 5 x

由于结构的对称性,上面的一排结点在轴向 上的位移应该相同
283

2. 无摩擦的界面
? 如果满足下列条件,则可用耦合自由度来模拟接触面:
–表面保持接触, –此分析是几何线性的(小变形) –忽略摩擦 –在两个界面上,节点是一一对应的.

?通过仅耦合垂直于接触面的移动来模拟接触. ?优点:
–分析仍然是线性的 –无间隙收敛性问题

284

3.

铰接

考虑一个2D的梁模型,每个节点上有三个 自由度ux、uy和rotz,A点为一铰链连接 。将同一位置节点的自由度ux、uy耦合起 来。

耦合可用来模拟力耦松 弛,例如铰链、无摩擦 滑动器、万向节

1

2

A

节点1和节点2 处于同一位置 ,但为于清楚 起见,在图上 分开显示。.

为了模拟铰接,将同一位置两个节点 的移动自由度耦合起来,而不耦合转 动自由度
285

用耦合施加循环对称性

在循环对称切面上的对应位置实 施自由度耦合。

286

三种建立耦合关系的方法

进入创建耦合关系的菜单路径: Main Menu: Preprocessor > Coupling / Ceqn > Couple DOFs 1. 拾取将要耦合的结点

4. 单击OK. 2. 单击OK

3. 输入耦合设置参考 号,选择自由度卷 标.
287

在零偏移量的一组节点之间生成附加耦合关系:
Main Menu: Preprocessor > Coupling / Ceqn > Gen w/Same Nodes

1. 输入现存耦合 设置的参考号.

2. 对每个设置指定 新的自由度卷标 . 3. 单击OK
288

在同一位置的节点之间自动生成耦合关系:
Main Menu: Preprocessor > Coupling / Ceqn > Coincident Nodes 1. 指定自由度卷标. 2. 指定节点位置的 容差

3. 单击OK

289

练习-用耦合关系来模拟接触
在此练习中,将用耦合/ 约束选项在两部分间产 生耦合DOF设置来模拟 接触问题
1. 2.

恢复数据库cpnorm.db1,并在 图形窗口中画单元. 在重合节点的所有节点对上建 立UY耦合关系
a. 选择耦合重合的结点.
b. 拾取UY

3.

求解并进行后处理

290

约束方程
约束方程定义节点自由度之间的线性关系 约束方程的特点 ? 自由度卷标的任意组合. ? 任意节点号. ? 任意实际的自由度方向――在不同的节点上ux可 能不同. 例 Constant = Coef1 * DOF1 + Coef2 * DOF2 + ...
291

四种约束方程的应用
1. 连接不同的网格:
? 实体与实体的界面

? 2-D或3-D
? 相同或相似的单元类型 ? 单元面在同一表面上,但结点 位置不重合

292

2. 连接不相似的单元类型: 壳与实体 垂直于壳或实体的梁.

? 建立转动自由度和移动自由度之间的关 系
293

3.

建立刚性区
? 在某些特殊情况下,全刚性区给出了 约束方程的另一种应用

? 全刚性区和部分刚性区的约束方 程都可由程序自动生成

4. 过盈装配
? 与接触耦合相似,但在两个界面之间 允许有过盈量或穿透 ? 典型方程: 0.01 = UX (node 51) - UX (node 251)
294

四种建立约束方程的过程
人工建立约束方程的菜单路径: Main Menu: Preprocessor > Coupling / Ceqn > Constraint Eqn

1. 输入常数项, 节点号,自由 度卷标和方程 系数.

2. 单击OK.
295

1. ..... 2. ..... 3. ..... Procedure

以现有的约束方程为基础生成约束方程:
1. 生成第一个约束方程: Main Menu: Preprocessor > Coupling / Ceqn > Constraint Eqn. 2. 生成其余的约束方程: Main Menu: Preprocessor > Coupling/Ceqn > Gen w/Same DOF.
生成的约束方 程数. 现存约束方 程中的节点 增量 生成的约束方 程的起始序号 ,终止序号和 增量

3. 选择 OK

296

1. ..... 2. ..... 3. ..... Procedure

通过“刚性区”来建立约束方程
Main Menu: Preprocessor > Coupling / Ceqn > Rigid Region >

拾取将要连在一起的结点,然后单击OK
1. 选择将要使用的

刚性区的类型( 自由度设置)

2. 单击OK

297

在相邻的区域生成约束方程:
1. 从网格较密的区域中选择节点 2. 从网格较稀的区域中选择单元. Main Menu: Preprocessor > Coupling / Ceqn > Adjacent Regions
3.

指定容差,此容 差作为单元区域 中最小单元长度 的比率.

4.

在约束方程中 将要使用的自 由度
298

5. 单击OK

Exercise

在此练习中,将使用约束方程 将具有不同单元类型和不同网 格的两部分连接起来。这两部 分分别是涡轮叶片段及叶片连 接的基座

Blade

1. 2. 3.

恢复数据库文件(eblade.db1) 并在图形窗口中显示单元.

Base

选择基座上的单元(mat2) 选择叶片底面上的节点
a.首先,unselect附在底座单元

上的节点(接第2步)
b.然后,在 4. 5.

位置Z=0处 reselect节点
6. .

在所选的相邻区域生成约束 方程. 选择everything 求解并进行后处理
299

Introduction to Finite Element Mechod (有限元方法介绍)
? 13.1 Numerical Methods (数值方法) ? 13.2 Basics of the Finite Element Method (有限元法基础) ? 13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数) ? 13.4 Virtual Work Formulation for Plane Elasticity (平面弹性虚功方程) ? 13.5 A Simple Example (简单例子) ? 13.6 FEM Code Application (有限元序程的应 用)

13.1 Numerical Methods (数值方法) Three numerical methods (三种数值方法)
Finite dfifference method (FDM) (有限差分法)

Finite element method (FEM) (有限元法)

Boundary element method (BEM) (边界元法)

13.2 Basics of the Final Element Method (有限元法基础)
Typical basic steps (典型基本步骤)
1. Discretize the body into a finite number of element subdomains (离散) 2. Develop approximate solution over each element in terms of nodal values (写出单元刚度矩阵) 3. Based on system connectivity, assemble elements and apply all continuity and boundary conditions to develop an algebraic system of equations among nodal values (形成总体 刚度矩阵并引入连续条件和边界条件) 4. Solve assembled system for nodal values; post process solution to determine additional variables of interest if necessary (求解节点值和后处理)

13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数)
Displacement interpolation functions (位移插值函数)
we wish to investigate procedures necessary to develop a linear approximation of a scalar variable u(x,y) over an element.(需要选 定一个位移模式,来表示单元中的位移函数) Assume a linear approximation: (假设为线性的)

u(x, y) ? c1 ? c2 x ? c3 y
For every node

u(x1 , y1 ) ? c1 ? c2 x1 ? c3 y1 u(x 2 , y 2 ) ? c1 ? c2 x 2 ? c3 y 2 u(x 3 , y3 ) ? c1 ? c2 x 3 ? c3 y3

13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数)
Displacement interpolation functions (位移插值函数)
the constants ci can be expressed in terms of the nodal values ui (常 数ci可用节点的值ui表示)

1 ??1u1 ? ? 2u2 ? ? 3u3 ? c1 ? 2 Ae 1 ??1u1 ? ? 2u2 ? ?3u3 ? c2 ? 2 Ae 1 ?? 1u1 ? ? 2u2 ? ? 3u3 ? c3 ? 2 Ae

Where :

? i ? x j y k ? xk y j ? i ? y j ? yk ? i ? xk ? x j

13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数)
So we obtain :

u ( x, y ) ?

? ??1u1 ? ? 2u2 ? ? 3u3 ?x

1 [??1u1 ? ? 2u2 ? ? 3u3 ? 2 Ae
3

? ?? 1u1 ? ? 2u2 ? ? 3u3 ? y ] ? ? ui? i ( x, y )
i ?1

Where Ψi are the interpolation functions for the triangular element given by : (Ψi是三角单元插值函数)

1 ?? i ? ?i x ? ? i y ? ? i ( x, y) ? 2 Ae

Lagrange interpolation functions (拉格朗日插值函数)

13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数)

1 ?? i ? ?i x ? ? i y ? ? i ( x, y) ? 2 Ae
It is noted that the form of the interpolation functions depends on the initial approximation assumption and on the shape of the element. (插值函数的形式取决于初始近似假设和 单元形状)

13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数)
description of element displacement filed (单元位移场 的表达) ? u1 ? ?v ? ? 1? ?u ? ?? 1 0 ? 2 0 ? 3 0 ? ?u2 ? ? ??? ? ? ? ?? ? ??? ? ? v ? ? 0 ? 1 0 ? 2 0 ? 3 ? ? v2 ? ?u 3 ? ? ? ? ? v3 ? ?
Where:

? i ( x, y) ?

1 ?? i ? ?i x ? ? i y ? 2 Ae

13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数)
description of element strain filed (单元应变场的表达)
?? ? ? ? x ? ? ?x ???0 ? ( x, y ) ? ? ? y ? ? ? ? ?? z ? ? ?? ? ? ? ?x ? 0? ? ? ? ?u ( x, y )? ? [? ]u ? [? ][? ]? ? [ B ][? ] ? ? ? ?y ? v( x, y ) ? ?? ? ?y ? ?
?? 3 ?x 0 ?? 3 ?y ? 0 ? ? ?1 ? ?? 3 ? 1 ? ? 0 ?y ? 2 Ae ? ? ?? 1 ?? 3 ? ? ?x ? ?

Where:
? ?? 1 ? ? ?x [ B] ? ? 0 ? ? ?? ? 1 ? ? ?y 0 ?? 1 ?y ?? 1 ?x ?? 2 ?x 0 ?? 2 ?y 0 ?? 2 ?y ?? 2 ?x 0

?2

0

?3
0

?1 0 ?1 ? 2

?2 ?2

?3

0? ?3 ? ? ?3 ? ?

13.3 Approximating Functions for TwoDimensional Linear Triangular Elements (二维线性三角形单元的近似函数)
description of element stress filed (单元应力场的表达)

?? ? ? [C]{? } ? [C][B]{?}
Where:
?C11 [C ] ? ? ?C12 ? ? 0 C12 C22 0 0 ? 0 ? ? C66 ? ?
E ? 2 …plane stress ? ? 1 ? ? C11 ? C22 ? ? E (1 ?? ) ? …plane strain ? ? (1 ? ? )(1 ? 2? ) E? ? …plane stress 2 ? ? 1 ? ? C12 ? ? E? ? ? …plane strain ? (1 ? ? )(1 ? 2? ) C66 ? ? ? E 2(1 ? ? )

…plane stress and plane strain

13.4 Virtual Work Formulation for Plane Elasticity (平面弹性虚功方程)
Virtual work formulation (虚功原理)

?

Ve

? ij?? ij dV ? ? Ti n?ui dS ? ? Fi?ui dV
Se Ve

Which can be written in compact form: (简洁写之)

{??}T ([K ]{?} ? {F} ? {Q}) ? 0 [ K ] ? he ? [ B]T [C ][B]dxdy Stiffness matrix ?
e

[ K ]{?} ? {F} ? {Q}

? Fx ? {F } ? he ? [? ]T ? ? dxdy Body force ?e ? Fy ? vector T T ? x? {Q} ? he ? [? ] ? ? ds Loading vector ?e ?Ty ?

13.5 A Simple Example (简单例子)
Example:
例题:write the stiffness matrix of the element show in figure. (求如图所示的单元在该坐标系下的刚度矩阵)

It is obviously :

xi ? a, x j ? 0, xm ? 0 yi ? 0, y j ? a, ym ? 0

13.5 A Simple Example (简单例子)

xj ai ? xm

yj , ym

1 yj bi ? ? , 1 ym (i, j, m)

1 xj ci ? ? , 1 xm

bi ? a, b j ? 0, bm ? ?a
1 xi 1 A ? 1 xj 2 1 xm yi yj 。 ym

ci ? 0, c j ? a, cm ? ?a a2 A? 2

13.5 A Simple Example (简单例子)
? k ii ? k ? ? k ji ?k ? mi k ij k jj k mj
1? ? 1? ? ? ? b b ? c c ? b c ? c b k im ? ? r s r s r s r s? Et ? k ? 2 2 ? ?, rs 2 k jm ? 4(1 ? ? )A ? ?c b ? 1 ? ? b c c c ? 1 ? ? b b ? ? r s r s r s r s ? ? 2 2 ? ? k mm ?

?r ? i, j, m; s ? i, j, m?
? ? ? ? ? ? ?。 ? ? ? 3? ? ? ? 2 ?

? 1 ? ? 0 ? ? ? 0 Et k? ? 2(1 ? ? 2 ) ? ? ? ?1 ? ? ?? ? ?

1? ? 2 1? ? 2 0 1? ? ? 2 1? ? ? 2

对 1? ? 2 0 1? ? ? 2 1? ? ? 2

1 ?? ?1

称 3? ? 2 1? ? 2

13.6 FEM Code Application (有限元程序的应用)
Ansys application (Ansys 的应用)
Example 1: Circular and circle Holes in a Plate Under Uniform Tension (中心 有孔的矩形板受拉)

FEM Mesh and load condition (网格划分和载荷条件)

Distribution of x-stress (X 向应分布)

13.6 FEM Code Application (有限元程序的应用)
Example 2: Circular Disk Under Diametrical Compression (受压圆盘)

FEM Mesh and load condition (网格划分和载荷条件)

Distribution of x-stress (X 向应分布)

13.6 FEM Code Application (有限元程序的应用)
Abaqus application (Abaqus 的应用)
Example 1: Circular and ellipse Holes in a Plate Under Uniform Tension (中 心有椭圆孔的矩形板左端固定,右端受拉)

FEM Mesh and load condition (网格划分和载荷条件)

Distribution of x-stress (X 向应分布)

Vocabulary(词汇)
Finite dfifference method (FDM) Finite element method (FEM) Boundary element method (BEM) 有限差分法 有限元法 边界元法

Triangular Elements
Virtual Work Node interpolation functions Ansys,Abaqus

三角单元
虚功 节点 插值函数 有限元分析软件


更多相关文档:

计算材料学_图文.ppt

计算材料学 - 粒子物理与核物理实验中的数 据分析 杨振伟 清华大学 第四讲:蒙

计算材料学 之 材料设计、计算及模拟_图文.ppt

计算材料学 之 材料设计、计算及模拟_电力/水利_工程科技_专业资料。计算材料学概述第四章材料设计、计算及模拟 1 主要内容 ?计算材料学的起源 ?计算材料学的...

计算材料学14-5-3_图文.pdf

计算材料学14-5-3 - 计算材料学 计算材料学第五章 电子结构理论 陈效双

计算材料学.doc

计算材料学 - 计算材料学 (Computational Materials S

计算材料学的现实重要性.doc

计算材料学的现实重要性 - 计算材料学的现实重要性 【摘要】随着科学技术的发展,

计算材料学的五个程序.pdf

计算材料学的五个程序_计算机软件及应用_IT/计算机_专业资料。三维行走 INT

计算材料学在复合材料中的应用.doc

计算材料学在复合材料中的应用 - 计算材料学在复合材料中的应用 蒋雯 3120140417 摘要: 随着计算机技术的发展,计算材料学成为复合材料领域越来越 重要的研究分析手段...

计算材料学-1_图文.ppt

获得实验难以实现的极端条件下(如高温、高压) 的材料结构与物性 计算材料学范围和层次 计算材料学 计算材料学主要有两种分类方法:一是按理论模型和方法分类, 二是...

计算材料学1_图文.pdf

计算材料学1 - 计算材料学 Computational Materials S

计算材料学-1._图文.ppt

计算材料学-1. - 计算材料学 计算材料学 Computational mat

计算材料学_图文.ppt

计算材料学 - 计算材料学 杨振华 第一性原理计算方法 ? 第一性原理方法是一种

计算材料学_图文.doc

计算材料学 - 摘要:本文通过第一性原理方法,从理论角度研究掺杂 ZnO 的电子

计算材料学的现状与发展前景.pdf

计算材料学的现状与发展前景 - 维普资讯 http://www.cqvip.co

计算材料学-1_图文.ppt

计算材料学-1 - 计算材料学 计算材料学 Computational mate

计算材料学基础136-200_图文.pdf

计算材料学基础136-200 - 页码,1/6 http://www.lishi

计算材料学_图文.ppt

计算材料学 - 计算材料学 第一章 引言 Performance Composi

计算材料学(第一性原理_密度泛函理论_分子动力学)-md_图文.ppt

计算材料学(第一性原理_密度泛函理论_分子动力学)-md_材料科学_工程科技_专

计算材料学2._图文.ppt

计算材料学2. - 计算材料学 Computational Materials

计算材料学.doc

计算材料学 - 得分 评卷人 沈阳师范大学 硕士研究生期末考试答题册 考试科目: 姓学名: 号: 计算材料学 刘芳慧 12545002 沈阳师范大学 材料物理与化学 2013...

计算材料学Fortran程序汇总.doc

计算材料学Fortran程序汇总_工学_高等教育_教育专区。用Fortran模拟

更多相关标签:
网站地图

文档资料共享网 nexoncn.com copyright ©right 2010-2020。
文档资料共享网内容来自网络,如有侵犯请联系客服。email:zhit325@126.com