镁铁质岩浆体系的矿物-熔体平衡

1.程序功能

程序TEQUIL可根据用户提供的镁铁质岩浆体系的矿物组合及熔体成分，自动检索出可供使用的温度计方法，计算出橄榄石、高钙辉石、低钙辉石、斜长石、钛铁矿、尖晶石5种固溶体矿物-熔体平衡的28个结晶反应的平衡温度（Nielsen＆Dungan,1983）。程序适用于拉斑玄武岩、碱性玄武岩、安山岩等w（SiO2）＜60%的镁铁质火山岩类。

2.方法原理

本程序采用Nielsen和Dungan（1983）关于低压下干镁铁质岩浆体系的矿物-熔体平衡的热力学模型。程序设计的方法是，首先，根据用户提供的镁铁质火山岩的斑晶矿物组合及化学成分，计算出各矿物相的晶体化学式（阳离子系数）。然后，根据用户提供的初始参数，即估计的压力、温度和相对氧逸度，对熔体相的Fe2O3和FeO含量进行修正。在此基础上，按照双晶格熔体结构模型，计算熔体相各组分的活度。按照所选择的结晶反应，计算出矿物-熔体的平衡温度。

3.程序结构

结晶岩热力学软件

4.使用说明

（1）输入格式

程序运行过程中，按照屏幕提示，依次提供下列参数：

IFN/OFN 输入/输出文件名

JFe3（j） 铁镁矿物Fe3+/Fe2+计算选择

TC 估计的平衡温度（℃）

PGPa 估计的平衡压力（GPa）

dlgfo2　相对于FMQ缓冲剂的氧逸度

文件输入格式为：I2（相数），A6（样品号）；

A3,13F6.2。

各变量的排列顺序依次为：Sams（相名称）、SiO2、TiO2、Al2O3、Cr2O3、Fe2O3、FeO、MnO、NiO、MgO、CaO、Na2O、K2O、P2O5。

用户根据实际需要，可增加SrO、BaO、Li2O等3个变量。此时，只需改变子程序MINORM中的文件读入语句及相应的格式语句即可。

每次计算的样品个数不限。

（2）输出格式

全部计算结果输出到文件OFN中。内容包括：各矿物相的化学成分，晶体化学式（阳离子系数），摩尔质量，按照自动检索出的矿物相（结晶反应）（马鸿文，1993a）计算的平衡温度。对于其中那些计算不出合理的lnK值的结晶反应（如某一组分含量或活度为零时），其输出结果表示为Teq（C）=0。

5.程序文本

结晶岩热力学软件

data Meth/'Mg'，'Fe'，'Ni'，'Mn'，3*''，

$　 'En'，'Fs'，'Di'，'Hd'，'Al'，'Ti'，'Cr'，

$　 'Ab'，'An'，5*''，

$　 'Il'，'Mg'，'Cr'，4*''，

$　 'Sp '，'F3'，'Cr'，'Al'，'MF'，2*''/

write（*,*）'Input/Output filename=？’

read（*,10）IFN,OFN

10 format（A）

open（unit=3,file=IFN,status＝'old'）

open（unit=4,file=OFN,status＝'unknown'） ****　Set necessary parameters and calculate mineral formula *************

do j=1,6

if（j.eq.4）goto 18

write（*,15）Mtit（j）

15 format（'Fe3+cal ibration for'，A3，'（Y/N）？'）

read（*,10）JFe3（j）

18　end do

20 read（3,25,ERR=110,END=120）nn,Sample

25 format（12,A6）

write（4,30）Sample

30 format（/5X，'Sample：'，A）

call MINORM（nn,m,CMP,Mtit,Sams,Cpfu,WAlopx,JFe3）

****　Calibrate Fe2O3/FeO ratio of the liquid ****************************

TK=1473.15

write（*,35）Sample

35 format（/'Sample：'，A，'P（GPa）=？'）

read（*,*）PGPa

write（*,*）'dlogfo2（FMQ）=？’

read（*,*） dlgfo2

sum=0

do j=1,m

OXD（j）=Cpfu（14,j）/Ncat（j）

sum=sum+OXD（j）

end do

do j=1,m

OXD（j）=OXD（j）/sum

end do

logfo2=buffo2（3,TK,PGPa）+dlgfo2

call Fe3clb（TK,PGPa,logfo2,m,OXD）

sum=0

do j=1,m

OXD（j）=Ncat（j）*OXD（j）

sum=sum+OXD（j）

end do

do j=1,m

OXD（j）=OXD（j）/sum

end do

****　Calculate lnKd ＆ T（C）of mineral-liq equilibria ********************

call ACTLIQ（m,4,10,OXD,XNF,XNM）

do i=1,nn

if（Sams（i）.eq.'olv'）then

if（Cpfu（4,1）.ge.1）then

XSi=l

elSe

XSi=Cpfu（4,1）

end if

Xoc t=Cpfu（4,6）+Cpfu（4,7）+Cpfu（4,8）+Cpfu（4,9）

XMg=Cpfu（4,9）/Xoct

XFe=Cpfu（4,6）/Xoct

XNi=Cpfu（4,8）/Xoct

XMn=Cpfu（4,7）/Xoct

lnKd（1）=alog（XMg*XSi**0.5）/（XNM（9）*XNF（1）**0.5）

lnKd（2）=alog（XFe*XSi**0.5）/（XNM（6）*XNF（1）**0.5）

if（XNi.le.0.or.XNM（8）.le.0）then

lnKd（3）=0

else

lnKd（3）=alog（XNi*XSi**0.5）/（XNM（8）*XNF（1）**0.5）

end if

lnKd（4）=alog（XMn*XSi**0.5）/（XNM（7）*XNF（1）**0.5）

do j=1,4

Teq（j）=Ai（j）/（lnKd（j）-Bi（j）-273.15

if（lnKd（j）.eq.0）Teq（j）=0

end do

write（4,50）（Meth（j,1）,Teq（j）,j=1,4）

50 format（'olv（'，A2，'）-liq Teq（C）＝'，F6.1）

else if（Sams（i）.eq.'cpx'）then

do j=1,m

OXD（j）=Cpfu（6,j）

end do

call PYXoct（m,OXD,XMgMl,XMgM2,XFeMl,XFeM2,XAlMl,AlMlts）

XCaM2=OXD（10）/（XMgM2+XFeM2+OXD（7）+OXD（10）+OXD（13）

XTi=OXD（2）/4

XAl=OXD（3）/4

XCr=OXD（4）/4

XSi=OXD（1）/2

if（XSi.ge.1）XSi=1

if（XMgM2.le.0）then

lnKd（5）=0

else

lnKd（5）=alog（XMgMl*XMgM2*XSi**2）/（XNM（9）**2*XNF（1）**2）

end if

if（XFeM2.le.0）then

lnKd（6）=0

else

lnKd（6）=alog（XFeMl*XFeM2*XSi**2）/（XNM（6）**2*XNF（1）**2）

end if

lnKd（7）=alog（XMgMl*XCaM2*XSi**2）$　 /（XNM（9）*XNM（10）*XNF（1）**2）

lnKd（8）=alog（XFeMl*XCaM2*XSi**2）

$　/（XNM（6）*XNM（10）*XNF（1）**2）

if（XAl.le.0.or.XNM（3）.le.0）then

lnKd（9）=0

else

lnKd（9）=alog（XAl/XNM（3）

end if

if（XTi.le.0.or.XNM（2）.le.0）then

lnKd（10）=0

else

lnKd（10）=alog（XTi/XNM（2）

end if

if（XCr.le.0.or.XNM（4）.le.0）then

lnKd（11）=0

else

lnKd（11）=alog（XCr/XNM（4）

end if

do j=5,11

Teq（j）=Ai（j）/（lnKd（j）-Bi（j）-273.15

if（lnKd（j）.eq.0）Teq（j）=0

end do

write（4,60）（Meth（j-4,2）,Teq（j）,j=5,11）

60 format（'cpx（'，A2，'）-liq Teq（C）＝'，F6.1）

else if（Sams（i）.eq.'opx'）then

do j=l,m

OXD（j）=Cpfu（5,j）

end do

call PYXoct（m,OXD,XMgMl,XMgM2,XFeMl,XFeM2,XAlMl,AlMlts）

XCaM2=OXD（10）/（XMgM2+XFeM2+OXD（7）+OXD（10）+OXD（13）

XTi=OXD（2）/4

XAl=OXD（3）/4

XCr=OXD（4）/4

XSi=OXD（1）/2

if（XSi.ge.l）XSi=l

lnKd（12）=alog（XMgMl*XMgM2*XS i**2）/（XNM（9）**2*XNF（1）**2）

lnKd（13）=alog（XFeMl*XFeM2*XSi**2）/（XNM（6）**2*XNF（1）**2）

lnKd（14）=alog（XMgMl*XCaM2*XSi**2）

$　 /（XNM（9）*XNM（10）*XNF（1）**2）

lnKd（15）=alog（XFeMl*XCaM2*XS i**2）

$　 /（XNM（6）*XNM（10）*XNF（1）**2）

lnKd（16）=alog（XAl/XNM（3）

lnKd（17）=alog（XTi/XNM（2）

if（XCr.le.0.or.XNM（4）.le.0）then

lnKd（18）=0

else

lnKd（18）=alog（XCr/XNM（4）

end if

do j=12,18

Teq（j）=Ai（j）/（lnKd（j）-Bi（j）-273.15

if（lnKd（j）.eq.0）Teq（j）=0

end do

write（4,70）（Meth（j-11,2）,Teq（j）,j=12,18）

70 format（'opx（'，A2，'）-liq Teq（C）＝'，F6.1）

else if（Sams（i）.eq.'plg'）then

Xoct=Cpfu（11,10）+Cpfu（11,13）

lnKd（19）=alog（Cpfu（11,13）/Xoct）/（XNF（2）*XNF（1）**3）

if（XNM（3）.le.0）then

lnKd（20）=0

else

lnKd（20）=alog（Cpfu（11,10）/Xoct）/（XNM（10）*（XNM（3）**2）

$　 *XNF（1）**2））

end if

do j=19,20

Teq（j）=Ai（j）/（lnKd（j）-Bi（j）-273.15

if（lnKd（j）.eq.0）Teq（j）=0

end do

write（4,80） （Meth（j-18,3）,Teq（j）,j=19,20）

80 format（'plg（'，A2，'）-liq Teq（C）＝'，F6.1）

else if（Sams（i）.eq.'ilm'）then

Xoct=Cpfu（3,4）+Cpfu（3,6）+Cpfu（3,9）

lnKd（21）=alog（Cpfu（3,6）/Xoct）/（XNM（2）*XNM（6））

lnKd（22）=alog（Cpfu（3,9）/Xoc t）/XNM（9）

lnKd（23）=alog（Cpfu（3,4）/Xoct）/XNM（4）

do j=21,23

Teq（j）=Ai（j）/（lnKd（j）-Bi（j）-273.15

end do

write（4,90） （Meth（j-20,4）,Teq（j）,j=21,23）

90 format（'ilm（'，A2，'）-liq Teq（C）＝'，F6.1）

else if（Sams（i）.eq.'spn'）then

Xoct=0

do j=2,5

Xoct=Xoct+Cpfu（2,j）

end do

lnKd（24）=alog（Cpfu（2,2）/Xoct）/（XNM（2）*XNM（9）**2）

if（Cpfu（2,5）.le.0.or.XNM（5）.le.0）then

lnKd（25）=0

else

lnKd（25）=alog（Cpfu（2,5）/Xoct）/XNM（5）

end if

lnKd（26）=alog（Cpfu（2,4）/Xoct）/XNM（4）

lnKd（27）=alog（Cpfu（2,3）/Xoct）/（OXD（3）**2）

lnKd（28）=alog（Cpfu（2,9）/Xoct）/（XNM（9）/（XNM（6）+XNM（9））

do j=24,28

Teq（j）=Ai（j）/（lnKd（j）-Bi（j）-273.15

if（lnKd（j）.eq.0）Teq（j）=0

end do

write（4,100）（Meth（j-23,5）,Teq（j）,j=24,28）

结晶岩热力学软件

6.计算实例

甘肃北祈连九个泉地区蛇绿岩套，拉斑玄武岩的单斜辉石-熔体平衡温度计算（马鸿文等，1994，未发表资料）。

某地玄武安山岩的橄榄石、斜长石-熔体平衡温度计算。

输入文件：exam47.dat

结晶岩热力学软件

输出文件：exam48.dat

结晶岩热力学软件

结晶岩热力学软件

结晶岩热力学软件

有限元商业软件ABAQUS提供了多种模拟裂纹开裂的方法，其中以cohesive element和XFEM最为有效。XFEM只需要添加预制裂纹，即能模拟裂纹的扩展，裂纹路径不受限制，可穿过单元传播。XFEM模拟开裂时，认为材料各处的断裂强度是一致，这很难模拟本身具有离散性的脆性材料裂纹扩展受强度影响。将cohesive单元批量嵌入单元网格之间，并赋予cohesive element服从weibull分布的断裂强度，即可实现裂纹沿任意路径扩展，此时裂纹的扩展受强度的影响。虽然用cohesive element模拟裂纹开裂，裂纹只能演单元边开裂，但其对脆性材料开裂的模拟更接近真实情况。ABAQUS CAE中只提供了简单添加cohesive element的方法，要实现cohesive element的批量添加，需要通过对inp文件的处理来实现。

Python是ABAQUS的脚本语言，可很方便的实现文本 *** 作。cohesive element的添加流程如下：

（1）在ABAQUS CAE中生成inp文件，并将需要添加cohesive element的单元建立集合，为描述方便，下文把此集合称作CO_IN_SET；

（2）用python程序读取inp文件，分别获得节点信息及单元信息，以及CO_IN_SET；

（3）获得CO_IN_SET中单元对应的节点，并找出其中节点的重复次数，重复节点大于2的即需生成新的节点，每一重复节点生成的新节点比其重复次数少1，记录新节点对应的单元；

（4）获得CO_IN_SET中单元对应的边的重复次数，其单元的边由单元的节点按逆时针连接形成，单元的重复边即为需要嵌入cohesive单元的地方；

（5）替换CO_IN_SET单元的节点为新的节点，并按照重复边形成cohesive单元；

（6）输出新的包含cohesive element的inp文件。

---