扩散过程

扩散过程（精选九篇）

扩散过程 篇1

网络化制造是现代制造业发展趋势之一, 并针对特定需求衍生出多种网络化制造组织形式, 如虚拟企业[1]、基于ASP的网络化制造模式[2]、协同制造链[3]等。但上述制造模式多为“松耦合”管理, 强调合作企业的自主性和自治性, 往往存在生产过程安全保密措施不到位、进度无法及时了解、质量问题无法跟踪等问题, 不能满足战时大型复杂武器装备快速响应、制造过程可控、可测、可靠的要求。扩散制造[4]联盟正是基于解决上述问题而建立的新型网络化制造模式。本文从资源配置的角度提出一种全过程可控的扩散联盟方法, 将耦合性概念引入扩散制造模式, 分析了扩散制造节点间的耦合特征, 提出了基于降耦策略和升耦策略的资源配置解决方案, 亦即相关关系紧密的制造任务由相关关系紧密的资源进行制造, 最终建立扩散制造联盟。

1 基于耦合性的扩散制造联盟

1.1 相关概念

目前国内外研究基于耦合性的复杂产品设计任务规划的人员较多[5,6,7], 而涉足制造耦合性研究的人员则较少, 为了更清晰、明确地论述扩散制造耦合特性及扩散制造联盟结盟过程, 本文首先定义以下相关概念:

定义1 网络化制造耦合性 网络化制造的生产任务由异地制造资源协同完成, 由于各企业的组织结构、生产习惯、信息化平台以及地域上的差异, 使得原本在单独企业内的连续制造过程, 移植到网络上后被拆分为若干离散的制造节点。两个 (或两个以上的) 制造节点间相互作用, 而彼此影响整个制造系统的现象即为网络化制造耦合性。

定义2 扩散资源 分布在扩散企业, 能够承担某项扩散制造任务的生产组织和生产要素, 按照任务粒度大小的不同, 扩散资源分为企业级、单元级、设备级等不同层次。

定义3 扩散节点 实施扩散制造过程中, 复杂产品的制造被分解为若干制造任务, 并分配到具有加工能力的外部企业进行生产, 外部资源与其制造任务相互匹配的称为制造节点。由定义可知, 制造节点由制造任务、制造资源两部分组成, 可形式化定义为

MNi::={MTi, MRi}

式中, MNi为制造节点;MTi为制造任务;MRi为制造资源。

1.2 扩散制造耦合性分析

网络化制造耦合性的产生是由制造任务的离散分布导致的, 扩散制造作为面向大型复杂武器装备的网络化制造模式, 制造节点间的相互影响程度更强, 耦合特性表现得更为突出, 其原因有以下两点:

(1) 制造过程复杂。由于军方对保密、产品质量、生产进度方面的要求, 各级供应商都应纳入扩散制造体系, 构建全过程可控的扩散联盟。各扩散资源承担不同粒度的扩散任务, 原材料、半成品、成品在制造节点间流动, 层层交付, 相互交织, 形成庞大而复杂的扩散制造网络。

(2) 制造任务复杂。大型复杂产品的研制是一个巨大、复杂且技术密集的系统工程。如构成飞机的零件数目多达数百万, 零件的制造不仅要保证自身质量达到要求, 而且还要满足不同粒度制造任务间的相互协调关系, 如同一零件不同加工工序的加工协调关系、零部件间的装配协调关系。

复杂的制造过程和任务使得复杂产品在扩散制造过程中的物料、信息交互量要远远大于一般的网络化制造模式, 如图1所示。

首先主企业向制造资源发布任务数据包, 包括生产起止时间、工艺数据、质量要求等;资源在生产过程中实时向主企业反馈生产进度、质量检测数据;主企业对数据进行分析, 并迅速发出调整指令, 以保证整个扩散系统的正常、顺利运行。再者, 在扩散制造实施过程中, 为保证物料、信息及时准确地在节点间传递, 扩散节点间需要相互沟通协调, 确保生产进度协调一致及产品质量满足要求。庞大的信息、物料传输对主企业的制造过程控制以及节点间的配合造成了诸多困难。倘若某一中间环节出现扰动, 必定会影响其他节点的任务进度及完成质量, 甚至影响整个扩散制造的正常实施。

因此, 关系紧密的扩散任务有必要匹配关系紧密的制造资源形成扩散制造节点, 使节点间的耦合性更强。信息沟通、物料流转更加顺畅, 有利于产品的制造、装配、功能调试以及后续的维修工作, 也有利于降低不确定扰动对制造过程的影响, 保证扩散制造在运行过程中的稳定状态和持续生产能力, 提升扩散制造过程的可靠及可控性。

然而, 扩散任务间相关关系复杂, 相关强度也参差不齐, 如何合理处理任务间的复杂关系, 以及如何进行相关资源配置则是构建扩散联盟的两个关键问题。扩散降耦策略亦即对制造任务弱相关关系进行分析, 对相对弱的关系进行撕裂, 形成“内紧密、外松散”的任务集, 以降低制造节点间的耦合度, 降低控制策略复杂性, 提高任务执行的独立性和并行性。扩散资源升耦配置则是在降耦完成的基础上实现基于资源相关性的扩散资源配置, 为同一任务集的扩散任务选择关系紧密的制造资源, 构建紧耦合的扩散制造联盟。

2 扩散制造结盟过程

2.1 任务-资源匹配

零部件种类繁多、制造过程复杂, 使得一个制造企业没有能力单独承担一个零件或部件的加工, 这也决定了扩散制造任务呈现多粒度 (部件、组件、零件、工序级任务) 的特点, 如图2所示。本文用MTi (i=1, 2, …, nT) 表示扩散任务, nT为扩散任务的总数。根据扩散任务的具体加工需求, 通过扩散制造系统检索相应的制造资源, 具备该任务加工能力的资源可能有若干个, 分布在不同的扩散企业内, 任务-资源形成1对多的匹配关系, 如图2中的①所示。本文用MRi_p (p=1, 2, …, nR) 表示第p个具备MTi制造能力的扩散资源, nR为扩散资源的总数。

2.2 扩散降耦策略

降耦策略包括两个步骤, 分别为扩散任务相关性分析和扩散任务分解, 如图2中的②所示。

2.2.1 任务相关性分析

扩散任务间的复杂相关关系是使扩散系统产生耦合性的原因, 为了明确说明任务相关性, 本文将相关关系分为两类进行说明。

第I类为物料传递关系。即不同粒度任务完成后在扩散企业间传递以进行下一步零件加工、部件装配任务。fi j∈{0, 1}表示物料间的传递关系, fi j=0, 表示任务i、j无流动关系;fi j =1, 表示任务i完成后转向下一任务j。

第II类为除物料传递关系外的任务相关关系。

协调配合关系 各种产品或部件都是许多零件有条件地装配在一起的, 需考虑定位、连接、尺寸大小和装配顺序等, 有装配关系的零部件间相关关系紧密。

加工工序关系 由于复杂零件的加工复杂性, 在很多情况下需要异地转包加工, 同一零件的各工序级加工任务间存在紧密的相关关系。

功能实现关系 共同实现产品某一子功能的零部件生产任务间相互协调配合, 以完成能量转换、信号转化、力的传播, 以及其他功能干涉, 关系紧密。

ITmi j表示扩散任务i、j第II类相关度评价标值, m=1, 2, 3, 分别为协调配合关系、加工工序关系和功能关系, 视其重要或密切程度, 并参考文献[8]中的6分划分标准, 设为ITmi j∈{无关系, 弱, 较弱, 中, 较强, 强}={0, 0.2, 0.4, 0.6, 0.8, 1}, 其中, 第Ⅱ类相关值可表示为

式中, rTi j为扩散任务i、j间第Ⅱ相关度值;ωTm为相关关系权重, m=1, 2, 3。

由第Ⅰ维及第Ⅱ维相关关系共同建立扩散任务i、j间的相关性模型, 第一象限内是fij, 表示物料在任务间的流向;第二象限内是ITi j, 表示协调配合、加工工序以及功能相关关系的强度值, 如图3所示。

例如有7个任务, 0为主企业总装任务;1为部件级任务, 2、3为该部件的零件, 受交货期限或加工能力限制需继续扩散到下一级扩散资源进行生产, 完成后转移到部件生产企业, 以完成部件的装配;4为零件级任务, 5、6为零件4扩散出去的工序级任务, 交由外部扩散资源进行加工。

扩散过程可由有向图表达, 如图4所示。图4中, 实线箭头方向表示任务完成后的物料传递流向, 而具有第Ⅱ类相关关系的任务间则用虚线连接。相关性模型矩阵如图5所示。

2.2.2 扩散任务分解

本节对任务耦合模型中第Ⅱ类相关关系进行模糊聚类[9], 即关系较强的任务划为一类, 而关系较弱或没有关系的任务则作为独立任务进行处理, 形成“内聚集、外松散”的任务集, 实现了基于任务相关性的扩散任务分解。第Ⅱ类相关值也由数字型转化为布尔型, 即r′Ti j∈{0, 1}, 任务间关系紧密, r′Ti j=1;关系松散, 则r′Ti j=0。

如图4a所示的7个扩散任务经过降耦优化后, 撕裂了2-4, 4-7, 1-7之间的相关关系, 而任务7则为独立任务, 最终完成扩散任务的分解, 形成3个任务集{1, 2, 3}, {4, 5, 6}, {7}。其矩阵分解过程如图6所示。

2.3 扩散升耦策略

扩散升耦策略是在降耦完成的基础上, 考虑任务相关性及制造资源相关性因素, 设定优化目标, 实现扩散资源配置, 完成扩散联盟的构建。升耦策略也分为两个步骤, 即资源相关性分析及资源优化配置, 如图2中的③所示。

2.3.1 制造资源相关性分析

主企业选择相关性强的制造资源组成联盟, 将其纳入主企业的制造体系, 这是提高制造节点间耦合性、增强制造过程可控性和可靠性的有效手段之一。本文从行业、隶属、历史合作、地理位置及自身制造能力描述资源相关性。

行业关系 行业相近的企业间的专业人员、设备、数字化平台等资源类型和生产习惯会有较大的相似性。接到任务后, 在对原有条件不做太大改动的情况下, 能迅速实现转产, 达到快速响应的目的。根据国家质量监督检验检疫总局发布的《国民经济行业分类》 (GB/T4754-2002) 建立行业结构树, 如图7所示。

参考文献[10]求解扩散资源间行业相关程度, 设完成任务i、j的资源MRi_p与MRj_q所属行业到距离最近的公共父功能节点的距离分别为lp、lq, p、q为扩散资源, 令$x=\frac{1}{2{}^{l{}_p+l{}_q-2}},x＜0.25$时对资源间相关性的影响较弱, 可忽略不计, 所以MRi_p与MRj_q行业相关值为

$r^{\prime }{}_{{\rm I}(i\_p,j\_q)}=\{\begin{array}{l}xx\ge 0.25\\ 0x＜0.25\end{array}$

求出所有资源的r′I (i_p, j_q) 值, 利用

$r{}_{\text{Ι}(i\_p,j\_q)}=\frac{r^{\prime }{}_{\text{Ι}(i\_p,j\_q)}-r^{\prime }{}_{\text{Ι}(i\_p,j\_q)(\min )}}{r^{\prime }{}_{\text{Ι}(i\_p,j\_q)(\max )}-r^{\prime }{}_{{\rm I}(i\_p,j\_q)(\min )}}$

对其进行规格化处理, 使得相关值rI (i_p, j_q) ∈[0, 1]。

隶属关系 存在隶属关系的资源间有相同的上级主管, 便于在实施扩散制造的过程中沟通协调, 有利于问题的快速查找、解决, 保证扩散过程的顺畅, 图8所示为隶属关系结构树。

参照行业相关性值的求解方法, 求得MRi_p与MRj_q隶属相关值rS (i_p, j_q) 。

历史合作关系主企业对以往合作过的扩散企业完成的历史任务从交付时间、完成质量及价格等方面进行评价, 历史合作相关性反映了扩散企业间合作关系的良好程度及潜在的再次合作的可能性, 按照专家打分的方法, 与历史合作关系的评价值由低到高依次为{0, 0.2, 0.4, 0.6, 0.8, 1};

地理位置关系地理位置也是影响企业间合作的因素, 两个扩散企业间的地理位置越近, 在一定程度上越能缩短运输时间, 节约运输成本;反之, 越增加运输时间和运输成本。制造资源到的运输时间和运输成本分别为tTr (i p, j q) 、cTr (i p, j q) 。

制造资源自身能力是确定合作意愿强弱的关键因素之一, 能力越强, 主企业选择其进行合作的机会越大;反之, 合作机会越小。本文给出3个评价企业自身能力的指标:IC={tM (i_p) , cM (i_p) , qM (i_p) }, 其中, tM (i-p) 、cM (i_p) 、qM (i_p) 分别为承担扩散任务i的制造资源MRi_p的加工时间、加工成本和加工产品合格率。

2.3.2 资源升耦配置

传统的网络化制造资源配置过程多考虑地理位置{tTr (i_p, j_q) 、cTr (i_p, j_q) }、企业制造能力{tM (i_p) , cM (i_p) , qM (i_p) }等因素, 在本文中将这些因素也列为制造资源的相关特性, 并在此基础上增加了资源行业、隶属以及历史合作相关性的指标, 使得升耦优化配置成为一个更为复杂的多目标优化问题。

主企业选择扩散企业时, 需要考虑扩散制造资源的产品加工质量、时间和成本等指标;此外, 主企业与制造资源的横向相关程度也是需考虑的重要因素。承担所有任务MTi的扩散资源MRi_p与主企业的横向相关程度可表示为

$\begin{array}{l}R{}_{R1}={\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}(}\omega {}_{(R1)1}\lambda {}_{(R1)1}r{}_{\text{Ι}(i\_p)}+\\ \omega {}_{(R1)2}\lambda {}_{(R1)2}r{}_{\text{S}(i\_p)}+\omega {}_{(R1)3}\lambda {}_{(R1)3}r{}_{\text{C}(i\_p)})(2)\end{array}$

自身能力评价值、加工时间及成本可表示为

$Q{}_{\text{Μ}}={\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}q}{}_{\text{Μ}(i\_p)}$ (3)

${\rm T}{}_{\text{Μ}}={\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}t}{}_{\text{Μ}(i\_p)}$ (4)

$C{}_{\text{Μ}}={\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}c}{}_{\text{Μ}(i\_p)}$ (5)

式中, rI (i_p) 、rS (i_p) 、rC (i_p) 分别为具备任务i制造能力的扩散资源p与主企业的行业、隶属、历史合作相关性相关值;ω (R1) 1、ω (R1) 2、ω (R1) 3为各项权重, ${\displaystyle \sum _{m=1}^3\omega }{}_{(R1)m}=1$;λ (R1) 1、λ (R1) 2、λ (R1) 3为标准化系数, 用以消除量纲上的差异。

经过降耦优化后, 形成多个内耦合、外松散的扩散任务集, 对于集内任务 (r′Tij=1) 进行升耦优化, 提高制造节点间的耦合程度。承担相关任务的资源间需有较强的相关性, 便于资源间相互协调, 保证产品的质量和进度, 所有相关任务选择的扩散资源MRi_p、MRj_q的横向相关性可表示为

$\begin{array}{l}R{}_{R2}={\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}{\displaystyle \sum _{j=1}^{n{}_{\text{Τ}}}\frac{1}{2}}}{r}^{\prime }{}_{\text{Τ}ij}(\omega {}_{(R2)1}\lambda {}_{(R2)1}r{}_{\text{Ι}(i\_p,j\_q)}+\\ \omega {}_{(R2)2}\lambda {}_{(R2)2}r{}_{\text{S}(i\_p,j\_q)}+\omega {}_{(R2)3}\lambda {}_{(R2)3}r{}_{\text{C}(i\_p,j\_q)})i\ne j(6)\end{array}$

式中, ω (R2) 1、ω (R2) 2、ω (R2) 3为资源相关性值的权重;λ (R2) 1、λ (R2) 2、λ (R2) 3为标准化系数, ${\displaystyle \sum _{m=1}^3\omega }{}_{(R2)m}=1$。

(3) 对于存在物料传递关系 (fi j=1) 的制造资源MRi_p、MRj_q间的运输时间、运输成本可表示为

${\rm T}{}_{\text{Τ}\text{r}}={\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}{\displaystyle \sum _{j=i}^{n{}_{\text{Τ}}}f}}{}_{ij}t{}_{\text{Τ}\text{r}(i\_p,j\_q)}i\ne j$ (7)

$C{}_{\text{Τ}\text{r}}={\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}{\displaystyle \sum _{j=1}^{n{}_{\text{Τ}}}f}}{}_{ij}c{}_{\text{Τ}\text{r}(i\_p,j\_q)}i\ne j$ (8)

综合上述分析建立的扩散制造资源配置优化数学模型为

$\begin{array}{l}\max F=\omega {}_1(\omega {}_{11}\lambda {}_{11}R{}_{R1}+\omega {}_{12}\lambda {}_{12}\frac{1}{{\rm T}{}_{\text{Μ}}}+\\ \omega {}_{13}\lambda {}_{13}\frac{1}{C{}_{\text{Μ}}}+\omega {}_{14}\lambda {}_{14}Q{}_{\text{Μ}})+\omega {}_2\lambda {}_{2}R{}_{R2}+\\ \omega {}_3\lambda {}_3\frac{1}{{\rm T}{}_{\text{Τ}\text{r}}}+\omega {}_4\lambda {}_4\frac{1}{C{}_{\text{Τ}\text{r}}}i\ne j(9)\end{array}$

${\displaystyle \sum _{m=1}^4\omega }{}_m=1$ (10)

${\displaystyle \sum _{m=1}^4\omega }{}_{1m}=1$ (11)

式中, ωm (m=1, 2, 3, 4) 分别为主企业选择指标、制造资源相关性指标、运输时间和运输成本的权重;ω1m分别为制造资源与主企业的相关性评价、加工时间、加工成本和质量的权重;λ1m为标准化系数。

时间约束为

ti_p<Ei_p-Si_p (12)

扩散资源p承担任务i的时间ti应小于任务的计划最早开始时间Si和最晚交付时间Ei。成本约束可表达为

${\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}c}{}_{\text{Μ}i\_p}+{\displaystyle \sum _{i=1}^{n{}_{\text{Τ}}}{\displaystyle \sum _{j=1}^{n{}_{\text{Τ}}}c}}{}_{\text{Τ}\text{r}(i\_p,j\_q)}f{}_{ij}\le C$ (13)

扩散任务加工成本cMi_p和运输成本cTr (i_p, j_q) 的总和应小于项目计划成本总和C。

3 扩散联盟构建平台应用

利用VB.Net进行了系统实现, 开发了扩散制造资源配置系统, 包括主企业对任务的录入、扩散企业录入资源信息、扩散任务分解以及升耦配置功能。现以某型号飞机零部件实施扩散制造生产说明资源配置系统运行过程。如图9所示, 某型号机翼扩散生产的部分零部件有蒙皮、机翼、摇臂、弹簧、副翼、伺服补偿器、拉杆、发动机、发动机机架、襟翼、发动机叶片。

首先, 对待扩散制造任务的相关性进行分析, 根据任务相关性准则建立相关性矩阵, 如图10所示;然后利用降耦分解方法, 消除弱耦合关系, 解耦后的扩散任务聚集为{1, 2, 5, 10}, {3, 4, 6}, {8, 11}, {7}, {9}。

扩散资源配置过程包括两个部分:①录入已注册资源的相关信息, 除基本T、C、Q之外, 还有行业、隶属、合作及地理位置等, 以确定资源间的相关度值, 表1为承担任务1、2的2个资源间的相关关系;②通过在配置平台上进行扩散制造资源的升耦配置, 实现扩散任务与制造资源的紧耦合匹配, 建立了面向某飞机型号的扩散制造联盟。扩散联盟构建成功后, 各参与资源可使用平台所提供的项目协作管理工具, 进行进度监控、工艺协同, 完成整个扩散联盟的正常运行。

4 结束语

本文研究了现有的网络化制造模式后, 针对武器装备等复杂产品的异地协同制造特征, 给出了网络化制造耦合性概念, 分析了扩散制造耦合性的特征;通过扩散制造降耦分解、升耦配置策略, 完成了基于耦合性的扩散制造联盟构建, 并开发了资源优化配置系统, 本项目已在上海某航天所进行了验证, 实现了13类零部件的快速扩散制造。扩散制造耦合模型的建立及其演化, 为复杂产品网络化制造过程分析提供了更为丰富的途径, 可以指导扩散制造资源优化配置, 提高其合理性和科学性。

参考文献

[1]Martinez M T, Fouletier P, Park K H, et al.VirtualEnterprise-organisation, Evolution and Control[J].Int.J.Production Economics, 2001, 74:225-238.

[2]Su Yan, Liao Wenhe, Guo Yu, et al.An ASP-basedProduct Customization Service Systemfor SMEs:aCase Study in Construction Machinery[J].International Journal of Enterprise InformationSystems, 2008, 4 (1) :1125-1130.

[3]吉锋, 何卫平, 王东城, 等.网络化制造环境下面向复杂零件的协同制造链研究[J].计算机集成制造, 2006, 12 (1) :72-77.

[4]俞烽, 廖文和, 郭宇.工艺快速扩散系统的研究[J].中国机械工程, 2007, 18 (19) :2322-2325.

[5]郭峰.机械产品设计过程的建模、评价与优化方法研究[D].杭州:浙江大学, 2007.

[6]Chen Shijie, Lin Li.Decomposition of InterdependentTask Group for Concurrent Engineering[J].Computer&Industrial Engineering, 2003, 44:435-459.

[7]Sun Chungyu, Chen Shijie, Lin L.A StructuredApproach to Measuring Functional Dependency andSequencing of Coupled Tasks in Engineering Design[J].Computer&Industrial Engineering, 2003, 45:195-214.

[8]朱剑英.智能系统非经典数学方法[M].武汉:华中科技大学出版社, 2004.

[9]王海军, 魏小鹏.面向规模化产品族的数值规划方法[J].计算机辅助设计与图形学学报, 2005, 17 (3) :473-478.

惰性气体输运过程中的扩散性质研究 篇2

惰性气体输运过程中的扩散性质研究

系统计算了惰性气体及其二元混合体系在不同温度下的.扩散性质,给出了计算结果,并与J.Kestin等人所计算的最佳值进行了比较,结果符合较好,但本文的方法简单而且使用的势函数解析性好.这就进一步从宏观扩散性质的计算证明了Tang-toennies势的准确性.

作 者：郭建军 石建 杨继先 作者单位：四川工业学院,四川,成都 610039刊 名：四川大学学报(工程科学版) ISTIC EI PKU英文刊名：JOURNAL OF SICHUAN UNIVERSITY(ENGINEERING SCIENCE EDITION)年，卷(期)：34(6)分类号：O552.32关键词：二元混合体系 扩散系数 对应态原理

基于跳扩散过程的欧式交换期权定价 篇3

关键词 跳扩散过程;交换期权;随机微分方程

中图分类号 O211．6 文献标识码 A

Pricing of European Exchange Options on Jump Diffusion Process Model

HUANG Shuangshuang1,HE Zhanbing2

( College of Mathematics, Physics and Information Engineering,Jiaxing University,Jiaxing 314001,China;

2.Hunan Mass Media Vocational Technical College,Changsha,Hunan 410100,China)

Abstract The problem of pricing european exchange options on a jumpdiffusion model was considered. This paper assumes that the stock price is driven by jumpdiffusion process, and the jump process is a homogenous poisson process. The differential equation of option value was derived with noarbitrage theory. By using the method of numeraire conversion, the exact formula for pricing exchange option was obtained.

Key words jumpdiffusion process , exchange option , stochastic differential equation

1 引 言

交换期权是一种期权持有人在到期日有权但不必须以一种资产交换另一种资产的合约［1］. Margrabe在1978年首次给出了在扩散模型中交换期权的闭式解［2］，他的工作是在BlackScholes模型的假设下完成的，而大量的金融统计数据表明，这样的假设与实际情形有很大偏差. 为了减小偏差，较之Margrabe的扩散模型，Merton的跳扩散模型［3］更符合实际，该模型把股票价格的运动过程分为两部分：其一是正常的连续价格波动，即因一些细小的信息的到达使得股票价格进行一些小波动，用布朗运动来刻画；其二是“非正常”的不连续的价格波动，即因一些比较重大的信息的到达使得股票价格进行较大的波动，用泊松过程来刻画.在服从泊松过程的跳扩散模型基础上，金融衍生资产定价问题是一个热点问题，文献［4］基于此模型用期权定价的鞅方法得出了障碍期权的定价公式. 本文假定股票价格过程遵循带有强度参数都为λ的时齐泊松过程的跳扩散过程，在文献［3］的基础上进行了扩展，借鉴了文献［3］和［5］中利用无套利理论的方法推导出期权满足的微分方程，然后运用文献［6］中变换计价单位的方法求解无套利条件下期权满足的随机微分方程，并与扩散模型进行了比较.

2 基本模型及假设

考虑一个具有3项资产（B(t),S1(t),S2(t)）的无摩擦金融市场，假定市场无套利，其中B(t)为无风险债券的价格过程，满足方程dB(t)=rB(t)dt，常数r为银行利率；Si(t)(i=1,2)为第i项风险资产（股票）的价格过程，其不确定性包含扩散和跳跃，假定不支付红利，其在t时刻的价格S1(t)和S2(t)分别满足随机微分方程［3］：

dS1(t)S1(t)=(μ1－λk1)dt+σ1dW1(t)+X1dN(t),(1)

dS2(t)S2(t)=(μ2－λk2)dt+σ2dW2(t)+X2dN(t),(2)

其中,常数μi(i=1,2)是第i个股票的期望收益率，常数σi(i=1,2)是没有跳跃发生时第i个股票收益的波动率；Wi(t)(i=1,2)是标准布朗运动，其相关系数为ρ；N(t)是强度参数为常数λ的时齐泊松过程；常数ki≡E(Xi)(i=1,2)，其中Xi（Xi＞－1，否则会出现非正价格）是第i个泊松过程发生跳跃时第i个股票价格的相对跳跃高度且服从独立同分布；W1(t)，W2(t)，N(t)，X1j，X2j相互独立，Xij(i=1,2)是第i个股票第j次的相对跳跃高度，是独立同分布的，且Xi0＝0；E(•)是无条件期望算子.

引理1随机微分方程（1）和（2）的解［7］分别为：

S1(t)=S1(0)exp ((μ1－λk1－12σ21)t

+σ1W1(t))∏Ntj=0(1+X1j),(3)

S2(t)=S2(0)exp ((μ2－λk2－12σ22)t

+σ2W2(t))∏Ntj=0(1+X2j),(4)

即在服从泊松跳过程的跳扩散模型下，股票价格的显示公式.

3 期权价值微分方程

记联系两种股票的交换期权在t时刻的价值为V(t)，设V(t)＝F(S1,S2，t)，其中F关于t一次连续可微，关于S1,S2二次连续可微，则由交换期权的定义，有

F(S1,S2,T)=(S2－S1)+，(5)

其边值条件为

F(S1,0,t)=0 .(6)

由式(1)和式(2)及It引理，期权的收益率可写成：

dV(t)V(t)=(μv－λk1v－λk2v)dt+σ1vdW1(t)

+σ2vdW2(t)+X1vdN(t)+X2vdN(t), (7)

其中,μv是期权的期望收益率；(σ1v,σ2v)是没有跳跃发生时期权收益的波动率； kiv=∑(Xiv) (i=1,2)，其中Xiv是服从独立同分布的第i个过程发生跳跃时期权的相对跳跃高度.

引理2(广义It公式)［8］设有跳扩散过程dx=adt+ bdW+ydq, 另有函数f(t,x)关于t一阶连续, 关于x二阶连续可导, 则

df=(ft+afx+12fxxb2)dt+bfxdW+Ydq,

其中Y=f(t,x+y)－f(t,x).

定理1 设由跳产生的风险为非系统风险，F(S1,S2,t)是联系于股票S1和股票S2在t时刻的交换期权价值，S1和S2分别满足式(1)和式(2)，则F(S1,S2,t)满足微分方程组：

Ft+12σ21S21F11+12σ22S22F22+ρσ1σ2S1S2F12 +(r－λk1)S1F1+(r－λk2)S2F2－rF+

λE(F(S1(1+X1),S2(1+X2),t)－

F(S1,S2,t))=0，F(S1,S2,T)=(S2－S1)+. (8)

证明：由It引理，有：

μv=［Ft+12σ21S21F11+12σ22S22F22+

ρσ1σ2S1S2F12+(μ1－λk1)S1F1+(μ2－

λk2)S2F2+λE(F(S1(1+X1),S2(1+X2),t)－

F(S1,S2,t))］/F(S1,S2,t).(9)

σ1v=σ1S1F1(S1,S2,t)F(S1,S2,t),(10)

σ2v=σ2S2F2(S1,S2,t)F(S1,S2,t), (11)

X1v=F(S1(1+X1),S2,t)－F(S1,S2,t)F(S1,S2,t), (12)

X2v=F(S1,S2(1+X2),t)－F(S1,S2,t)F(S1,S2,t),(13)

其中F的下标表示偏微分，E(•)是无条件期望算子.

考虑一个包含两种股票S1,S2和期权V的资产组合，令其比例分别为Δ1、Δ2和Δ3，Δ1+Δ2+Δ3=1，记P(t)为组合在t时刻的价值，那么组合的期望收益率可以写成：

dP(t)P(t)=(μp－λk1p－λk2p)dt+σ1pdW1(t)+

σ2pdW2(t) +X1pdN(t)+X2pdN(t) , (14)

其中，μp是组合的期望收益率；(σ1p,σ2p)是没有跳跃发生时组合收益的波动率；kip=∑(Xip)(i=1,2)，Xip是服从独立同分布的第i个过程发生跳跃时组合的相对跳跃高度.

由式(1)、(2)及式(9)式，有：

μp=Δ1μ1+Δ2μ2+Δ3μv , (15)

σ1p=Δ1σ1+Δ3σ1v,(16)

σ2p=Δ2σ2+Δ3σ2v ,(17)

X1p=Δ1X1+Δ3F(S1(1+X1),S2,t)－F(S1,S2,t)F(S1,S2,t),(18)

X2p=Δ2X2+Δ3F(S1,S2(1+X2),t)－F(S1,S2,t)F(S1,S2,t). (19)

选取Δ1=Δ*1，Δ2=Δ*2和Δ3=Δ*3,使得Δ*1σ1+Δ*3σ1v=0和Δ*2σ2+Δ*3σ2v=0. 记此时组合的价值为P*，把式(16)、(17)代入式(14)，得此组合的期望收益率：

dP*(t)P*(t)=(μ*p－λk*1p－λk*2p)dt+X*1pdN(t)+

X*2pdN(t). (20)

如果由跳产生的风险为非系统风险，由资本资产定价理论，组合的期望收益率等于无风险利率r，即μ*p=r，因此得到方程组:

Δ*1μ1+Δ*2μ2+Δ*3μv=r,Δ*1σ1+Δ*3σ1v=0,Δ*2σ2+Δ*3σ2v=0.(21)

将式(9)~(11)代入方程组(21)，得到F满足微分方程：

Ft+12σ21S21F11+12σ22S22F22+ρσ1σ2S1S2F12+

(r－λk1)S1F1+(r－λk2)S2F2－rF+

λE(F(S1(1+X1),S2(1+X2),t)－

F(S1,S2,t))=0. (22)

证毕.

注1 如果λ=0，微分方程组(8)即股票价格服从连续过程的资产交换期权价值方程组：

Ft+12σ21S21F11+12σ22S22F22+

ρσ1σ2S1S2F12+rS1F1+rS2F2－rF=0,

F(S1,S2,T)=(S2－S1)+.(23)

4 交换期权定价公式

定理2 在股票价格过程服从式(1)~(2)的跳扩散模型中，由式(5)定义的交换期权的定价公式为

F(S1,S2,t)

=∑+

n=0e－λ(T－t)(λ(T－t))nn!E［S2e－λk2(T－t)Φ(d1)

•∏nj=0(1+X2j)－S1e－λk1(T－t)Φ(d2)∏nj=0(1+X1j)］, (24)

其中E(•)是期望算子，

d1=1σ2(T-t)［ln（S2S1e-λ(k2-k1)(T-t)•

∏nj=01+X2j1+X1j)+12σ2(T-t)］,

d2=d1－σ2(T－t),(25)

其中σ2=σ21－2ρσ1σ2+σ22.

证明 为了求解微分方程组(8)，做一个变换，令Z=S2S1,

U(Z,t)=F(S1,S2,t)S1.(26)

即引进新的概率测度Q，满足dQdP=S1(t)S1(0)B(t)，则有

FS1=U－ZUZ，

FS2=UZ，

2FS21=Z2S1•2UZ2，

2FS22=1S1•2UZ2，

2FS1S2=－ZS12UZ2.

为了简化，记

F1=U－ZUz，

F2=Uz，

F11=Z2S1Uzz，

F22=1S1Uzz，

F12=－ZS1Uzz,(27)

其中F的下标与U的下标表示偏微分. 将式(26)~(27)代入方程组(22)，容易验证U(Z,t)满足一维随机微分方程组：

Ut+σ22Z2UZZ－λ(k2－k1)ZUz+

λ(1+k1)E(U(1+X21+X1Z,t)－U(Z,t))=0,U(Z,T)=(Z－1)+, (28)

其中，σ2=σ21－2ρσ1σ2+σ22. 将方程组(28)与文献［3］中关于跳扩散模型中的欧式期权定价方程比较，易知可以把U(Z,T)看作是以虚拟资产Z(t)为标的资产的、敲定价格为1的欧式期权的价格，且在该虚拟市场中利率r=0，虚拟资产Z(t)的波动率为σ，泊松过程的强度参数为λS1=λ(1+k1)，且发生跳跃时组合的相对跳跃高度为X2－X11+X1.则由文献［3］中相应的定价公式，即知微分方程(28)的解为：

U(Z,T)=∑+

n=0e－λ(1+k1)(T－t)(λ(1+k1)(T－t))nn!•

EQ［Ze－λ(k2－k1)(T－t)∏nj=01+X2j1+X1jΦ(d1)－Φ(d2)］

=∑+

n=0e－λ(T－t)(λ(T－t))nn!E［Ze－λk2(T－t)Φ(d1)•

∏nj=0(1+X2j)－e－λk1(T－t)Φ(d2)∏nj=0(1+X1j)］,(29)

其中,Φ(x)=12π∫x－

et22dt，约定∏0j=1Xij=1，EQ(•)是在以S1(t)为新的计价单位的概率测度Q下的期望算子，E(•)是原给定概率测度P中的期望算子，且d1，d2分别为

d1=1σ2(T－t)［ln (S2S1e－λ(k2－k1)(T－t)•

∏nj=01+X2j1+X1j)+12σ2(T－t)］,

d2=d1－σ2(T－t),

其中σ2=σ21－2ρσ1σ2+σ22.

综合式(26)、式(29)，得到结论：在股票价格过程服从式(1)、式(2)的跳扩散模型中，由式(5)定义的交换期权的定价公式为

F(S1,S2,t)=∑+

n=0e－λ(T－t)(λ(T－t))nn!•

E［S2e－λk2(T－t)Φ(d1)∏nj=0(1+X2j)－

S1e－λk1(T－t)Φ(d2)∏nj=0(1+X1j)］,

其中E(•)是期望算子，d1,d2由式（25）定义.

证毕.

注2 如果λ=0，即模型简化为扩散模型，则期权定价公式简化为

F(S1,S2,t)=S2Φ(d1)－S1Φ(d2)，(30)

其中

d1=1σ2(T－t)［ln S2S1+12σ2(T－t)］，

d2=d1－σ2(T－t).

此即Margrabe在文献［2］中得到的结果.

4 结 论

本文在不支付红利的前提下，求解了股票价格服从带时齐泊松跳的跳扩散模型的交换期权定价问题，运用了无套利理论推导出期权价值的微分方程，利用变换计价单位的方法把交换期权的定价问题转化成单个的期权定价问题，从而得到交换期权的显示定价公式.但在现实金融市场中，股票价格可能是支付红利的，股票价格的跳可能并不一定是泊松跳过程，波动率也可能不是常数，很多更为复杂的情形还有待于进一步去研究.

参考文献

［1］ 姜礼尚. 期权定价的数学模型和方法［M］. 北京：高等数学教育出版社，2003.

［2］ William Margrabe. The value of an option to exchange one asset for another［J］. The Journal of Fiance, 1978, 33(1): 177-186.

［3］ Robert C Merton. Option pricing when underlying stock returns are discontinuous［J］. Journal of Financial Economics, 1976,(3):125-144.

［4］ 王莉, 杜雪樵.跳扩散模型下的欧式障碍期权的定价［J］. 经济数学, 2008, 25(3): 248-253.

［5］ 姚小义，邹捷中，陈超. 股票价格服从跳-扩散过程的资产交换期权定价模型［J］.长沙铁道学院学报，2002, 20(1)：4-9.

［6］ 钱晓松. 跳扩散模型中交换期权的定价［J］. 扬州大学学报:自然科学版，2004, 7(1)：9-12.

［7］ Knut K Aase. Contingent claims valuation when the security price is combination of an ito process and a random point process［J］. Stochastic processes and their Applications, 1988, 28(2):185-220.

［8］ 张利兵, 潘德惠. 标的股票价格服从跳跃-扩散过程的期权套期保值率确定［J］.系统工程理论方法应用, 2005, 14(1): 23-27.

摩擦焊接过程中的扩散现象 篇4

关键词：摩擦焊接,扩散行为,冶金现象,扩散

在材料连接过程中, 原子扩散行为是十分普遍的, 但是扩散的快慢等方面是差异很大的, 摩擦焊接过程是十分迅速的, 在几秒内就可以完成, 而扩散焊接过程却十分的缓慢, 需要几个小时的时间。因此, 对摩擦焊过程和扩散焊过程中的原子扩散行为做了对比, 从而研究其扩散过程的不同。

首先对摩擦焊过程中的原子扩散行为进行模拟, 选用的材料为面心立方结构Ti。通过对摩擦焊过程的分析, 无论是线性摩擦焊、惯性摩擦焊还是搅拌摩擦焊, 在摩擦过程中, 界面附近生成一层塑性流动层, 在高温下有着无序化的倾向, 因此, 直接采用理想的面心立方、体心立方或密排六方结构来模拟顶端阶段的构型是不合理的, 因此, 首先施加周期性振动周期为20 PS, 振幅为2 nm的相对振动作用产生塑性流动层, 如图1所示。

在产生塑性流动层的基础上, 对模型施加一定的顶锻压力进行扩散。在此过程中, 环境温度设为1200 K, 压力为150 MPa, 如图2所示。

如图2所示, 随着扩散时间的增加, 塑性层中原子扩散距离逐渐增加, 扩散原子数目也随之增加。然后计算了界面附近原子的均方位移曲线, 计算时只考虑振动结束之后阶段的原子扩散, MSD曲线如图3所示, 可见在振动之后, 塑性层中的原子具有长程扩散行为。

如图3所示, 塑性层中原子的均方根位移随着扩散时间的增加逐渐增加, 曲线出现一定的上下波动, 表明原子除了做晶格振动之外还具有扩散行为。为了研究摩擦焊中原子的扩散行为与普通扩散焊中原子扩散行为的异同, 下面构建了扩散焊过程的分子动力学模型, 材料仍然是Ti, 初始构型如图4所示, 温度为1200 K, 两端平均压力为150 MPa, 持续600 PS。通过扩散焊分子动力学模拟显示, 随着扩散时间的增加, 界面附近发生原子扩散现象, 最终构型如图4所示。

如图4所示, 随着扩散时间的增加, 界面附近出现了原子扩散现象, 左方模型中的原子进入右方模型中, 右方中的部分原子通过扩散进入左方模型之中。随着扩散时间的增加, 原子的扩散距离随之增加;同时, 发生扩散现象的原子数目也相应增加。与摩擦焊中的原子扩散相比, 可以看出, 普通扩散焊中原子扩散的时间虽然增加, 然而原子扩散的距离以及数目均相对较小。下面对扩散焊中原子扩散焊中原子的均方根位移进行计算, 然后通过与摩擦焊中原子扩散的均方根位移曲线进行比较, 从而对摩擦焊中原子扩散特点进行解释。普通扩散焊原子扩散均方根位移曲线如图5所示。

如图5所示, 在扩散焊中, 原子的均方根位移曲线随着扩散时间的增加而缓慢增加。在100 PS之前, 由于左、右模型靠近时界面原子扩散剧烈, 原子的均方根位移曲线上升较快, 当扩散一定时间之后, 原子均方根位移曲线上升缓慢, 而且出现较强的波动现象, 表明原子既有晶格附近的振动现象, 也出现了一定程度的扩散现象。根据爱因斯坦公式, 可以得出均方位移曲线的斜率与扩散系数成正比关系。

上述计算结果表明:摩擦焊中原子扩散速度远大于普通扩散中原子的扩散速度, 摩擦焊中原子发生了超扩散现象, 能够在短时间进行扩散, 从而有助于界面消失以及缺陷的弥合, 从而达到可靠的连接。

参考文献

[1]张昭, 刘会杰.搅拌头形状对搅拌摩擦焊材料变形和温度场的影响[J].焊接学报, 2011 (3) :5-8, 113.

扩散过程 篇5

期权是金融市场中一种常见而重要的金融衍生工具。期权的持有者拥有可以在合同规定的期限内按照事先约定的敲定价格K处置 (购买或者出售) 标的资产的权利。障碍期权是奇异期权的一种, 它的最终收益结构不仅依赖于标的资产在期权到期日的价格ST, 而且还要考虑标的资产价格在期权生命期内是否达到某个特定值, 即障碍价格H。障碍期权依标的资产价格达到障碍价格H后期权合约的状态, 可分为两种, 当资产价格S到达障碍价格H时, 期权合约生效, 这种期权称为敲入障碍期权;当资产价格S到达障碍价格H时, 期权合约失效, 这种期权称为敲出障碍期权。根据障碍价格H和敲定价格K的关系, 可以将敲入或者敲出障碍期权分为上升、下降两种:当H>K时, 称为上升;当H

本文以上升敲入看涨平方障碍期权为例进行研究, 即未定权益为 (ST2-K) +I{τH≤T}, H>K, 其中τH=inf{t>, 0St=H}。

二、模型假设

1.假设市场上只有两种有价资产, 分别为无风险资产B (t) 和期权的标的资产St;且其价格分别满足, , 其中r是无风险利率, μ是标的资产不发生跳扩散时的期望收益率, σ是标的资产的波动率, 是标的资产价格发生跳跃时的相对跳跃率, 为了避免资产价格出现负值, 规定ξ≥-1;{Nt, t∈[0, ∞) }是 (Ω, F, P) 上的标准Wiener过程, {Nt, t∈[0, ∞) }是一个Poisson过程, 强度参数为λ, Wt与Nt相互独立;其中r, µ, σ, ξ都设为常数。

2.假定标的资产价格St对无风险资产Bt的贴现过程在P中为鞅, 令, 根据Girsanov定理, 是P*测度下的标准Wiener过程。再令, , 其中σσα222-r=, 根据Girsanov定理, 是Q测度下的标准Wiener过程。则标的资产价格过程为

三、公式推导

对上升敲入看涨平方期权的价格进行贴现, 然后在P*测度下求期望, 既得该期权的价格, 因此

将I1和I2带入 (1) 式即得。

参考文献

[1]陈盛双, 杨云霞.连续平方障碍期权的定价[J].统计与决策, 2006 (14) .

[2]王莉, 杜雪樵.跳扩散模型下的欧式障碍期权的定价[J].经济数学, 2008 (03) .

[3]邢晓芳, 周圣武, 江龙, 王前.基于跳扩散模型的欧式下降敲入期权定价研究统计与决策, 2010 (17) .

[4]杨云锋, 夏小刚, 杨秀妮.跳跃扩散过程的期权定价模型[J].数学的实践与认识, 2010 (06) .

扩散过程 篇6

材料变形过程中组织结构发生剧烈变化, 其直接影响到原子的扩散行为。晶界作为原子扩散的通道之一, 对原子的扩散有着重要影响, 因此, 下面对晶界的作用进行了分子动力学模拟。

材料选用的是面心立方结构Ti, 为了构造晶界, 我们采用将模型按照一定的方向旋转的方法。初始构型如图1 (a) 所示, 然后将构型分为左右两组, 分别旋转一定的方向, 为满足右手定则, 将左侧部分按x方向 (-1 1 0) , y方向 (1 1-2) , z方向 (-1-1-1) 旋转, 而右侧部分按x方向 (1-1 0) , y方向 (-1-1 2) , z方向 (-1-1-1) 旋转, 相对旋转角度约为15°并由此构造了晶界, 构型如图1 (b) 所示, 图中不同的颜色代表不同的原子配位数。

2 晶界对原子扩散的影响

基于之前构建的晶界模型, 在1 000 K下研究晶界的演化过程, 并对晶界附近三层原子的扩散行为进行研究, 以此研究晶界处原子的扩散行为与完整晶体中原子扩散行为的差异。不同时刻晶界的演化过程如图2所示。

在1 000 K条件下, 通过对晶界处原子扩散行为的模拟结果可知, 相对于空位以及位错缺陷, 晶界缺陷消失所需时间更长, 约250 PS才能完全消失;随着扩散时间的增加, 大尺寸的晶界逐渐发生弯曲, 局部弥合, 形成多个小尺寸的局部缺陷, 并逐渐消失;靠近模型边缘处晶界消失速度大于模型中部晶界消失的速度。为了研究晶界对原子的扩散行为的影响作用, 对原子扩散的均方根位移进行计算, 如图3所示。

通过计算表明晶界处原子的扩散均方根位移曲线斜率约为0.8, 而相同条件下完整晶体中原子的扩散均方根位移斜率约为0.05, 即晶界处原子的扩散系数比完整晶体中原子扩散系数大一个数量级, 表明晶界处的原子扩散速率远远大于完整晶体中原子的扩散速度。

参考文献

[1]邹银辉.煤岩体声发射传播机理研究[D].山东:山东科技大学, 2007.

扩散过程 篇7

焰熔法是从熔体中生长蓝宝石、钛酸锶、金红石等光学单晶体的一种方法[1,2,3,4]。由于晶体生长是通过氧化物粉末在高温的氢氧火焰中熔化后结晶所得,其生长环境为氧化或还原气氛。晶体在高温氧化或还原气氛中生长,容易形成负离子空位而产生色心[5,6];生长过程中晶体周围的温度梯度和生长结束后结晶炉内快速降温都使得晶体内部产生严重的热应力,因此利用焰熔法生长的晶体一定要进行退火以排除应力和色心[7,8,9]。

在用焰熔法生长钛酸锶单晶体过程中,生长室内的气氛为富氢状态,即还原气氛。生长过程中,晶体周围的氧分压远小于钛酸锶晶体内的平衡氧分压,在晶体内部产生大量的氧空位而产生色心,因此,需要通过在氧化气氛条件下退火消除色心。实践表明,钛酸锶晶体在退火过程中,氧原子会通过晶体表面缓慢地向内部进行扩散。温度对氧原子在晶体内的扩散速率的影响非常大,氧原子在晶体内的扩散深度和浓度分布不仅与温度有关,而且还与退火的保温时间有关[10,11,12,13]。晶体的升、降温速率对晶体的热应力的去除和产生都会有很大的影响。而晶体内热应力、氧原子的浓度和扩散深度都决定着晶体的质量。因此,本研究采用有限元数值分析方法研究退火过程中保温时间和温度对氧原子在晶体中扩散行为的影响规律,为进一步研究晶体退火工艺和退火机理提供有价值的理论依据。

1有限元模型建立

1.1晶体结构模型

采用焰熔法生长的钛酸锶单晶体如图1所示,其中图1 (a)、(b)分别为晶体退火前和退火后状态,晶体尺寸总长为70mm,其中扩肩(Φ4~20 mm)长度为30 mm,等径生长(Φ20~30mm)长度为40mm。针对晶体具有的轴对称结构特征,考虑模型的对称性,选取如图2所示的二维计算模型。

1.2基本方程及边界条件

菲克定律是描述气体扩散现象的宏观规律。菲克第一定律指出,在单位时间内通过垂直于扩散方向的单位截面积的扩散通量J与该截面处的浓度梯度成正比,数学表达式为:

式中:D为扩散系数,m2/s;C为扩散物质的体积浓度,kg/m3;dC/dx为浓度梯度,“-”号表示扩散由高浓度区向低浓度区扩散;J为扩散通量,kg/(m2·s)。

在数学形式上,与菲克第一定律具有很好相似性的是傅里叶定律,即:

菲克第一定律只适用于稳态扩散。实际上,大多数扩散过程都是在非稳态条件下进行的,而菲克第二定律指出在非稳态扩散过程中,在距离x处,浓度随时间的变化率等于该处的扩散通量随距离变化率的负值,对于钛酸锶晶体退火过程中氧原子的扩散就是非稳态扩散过程,根据图2建立二维模型,晶体内氧原子扩散浓度可表示为:

式中:C为扩散物质的体积浓度,kg/m3;t为扩散时间,s;x为距离,m。

对于瞬态热传导方程:

当ρ=1,c=1时,瞬态热传导方程可转换为:

在数学形式上,与非稳态质量扩散方程具有很好的相似性。通过对式(1)、(2)、(3)、(5)进行对比分析,变量C与t、D与λ具有相似性。因此,对非稳态质量扩散方程的求解可以采用瞬态热传导方程的求解方法。本研究采用Ansys有限元软件的热分析模块进行质量扩散方程的数值计算,热分析模块计算的导热系数用质量扩散系数替代,比热容和密度分别设置为1。

晶体在炉内进行退火时,由于向退火炉内通入氧气,氧气对晶体表面的质量扩散过程可以类比成晶体在高温流体中的对流传热过程,而且质量扩散系数D很小,相当于对流传热过程中的导热系数λ很小。从非稳态传热的角度分析, 即描述晶体内部导热热阻与界面上对流热阻之比的Bi数很大,可将边界上的温度当作常数处理。对于质量扩散方程来说,可以将边界上的扩散物质的浓度设为常数,而且在炉内氧气气氛条件下,晶体表面的氧原子浓度远大于晶体内氧原子空位浓度。因此,对于非稳态质量扩散方程来说,晶体外表面的求解边界条件可确定为氧原子的浓度:C=C0=1,晶体中心线的边界条件为对称轴;在退火开始时,由于还没有氧原子扩散进入晶体内而占据氧原子空位,因此,晶体扩散的初始条件可确定为C=0。

2计算结果与分析

2.1氧浓度分布特征

根据实验条件,当晶体直径为30mm时,在900K条件下的退火时间为50h。通过数值计算得到50h后晶体内氧浓度分布,计算结果如图3所示。

从图3中可以看出,当晶体经过50h退火后,晶体表面浓度最大,为100%C0,内部离表面距离最大点的氧浓度最小,为95.7%C0。在晶体径向方向浓度差为0.5%时的浓度变化深度基本相等,表明氧浓度梯度为常数,说明晶体内氧浓度分布基本上达到了稳定状态,认为晶体退火已经完成。 这时晶体内氧扩散服从菲克定律,氧扩散通量与浓度梯度成正比。图4为退火过程中不同时刻晶体内氧的浓度分布。 从图4中可以看出,随着退火时间的延长,氧从晶体表面逐渐向内部扩散,晶体细小的根部和表面的氧扩散速度很快。 当退火时间为10h时,仅仅只在晶体根部和上部角位置氧的扩散比较明显。在30h内晶体中心氧浓度随时间的延长而增加得很快,但在晶体表面的扩散速度却较慢。结合图5中可以分析得到,在退火过程中由于氧很快地通过表面扩散到晶体表面A点使其浓度增加,导致其浓度梯度和扩散通量急剧下降,然后一直维持在一个很小的值基本不变。而对于B点,在退火初期很短的时间内,氧还没有扩散到此点,其浓度为0,没有扩散通量,而随着周围浓度逐渐增加,浓度梯度也增加而使扩散通量增大。当退火时间达到11h后,B点浓度梯度开始逐渐下降而使得扩散通量逐渐减小。

从图6可以很明显地分析得到B点氧浓度随时间的变化规律和退火深度(指从A点沿轴线向晶体内部方向,氧浓度达到96.5%C0的晶体厚度)。在退火前4h内,中心B点基本上还没有氧的扩散,5~30h内中心B点氧浓度增加迅速,30h后氧扩散速度有所减缓。这是因为退火前期,晶体表面和中心之间的浓度差大,扩散驱动力较强,所以中心B点氧浓度增加较快。对于扩散深度而言,在退火前期,从晶体表面扩散进来的氧要迅速地扩散到大面积的中心低浓度区域中,导致表面上氧积聚起来较难。随着退火时间的延长,晶体内部氧浓度逐渐增大,从表面往中心的氧扩散通量逐渐减小,导致表面浓度急剧增大,因此扩散深度在退火后期增大非常快。

2.2温度对氧扩散的影响

由于温度对扩散系数的影响较大(700K、800K、900K下的扩散系数分别为2.28×10-10m2/s、3.35×10-10m2/s、 5.15×10-10m2/s),而扩散系数又影响着晶体内氧的扩散。 因此,温度对氧扩散的影响非常明显。图7为不同温度条件下退火50h后晶体内氧浓度分布计算结果图。

从图7中可以看出,在700K、800K、900K条件下退火晶体内最小浓度分别为65%C0、83.9%C0和95.8%C0,700 K退火50h后只在晶体表面很浅的一层退透,厚度约为1mm,800K条件下的退火深度约为3mm,而900K时则基本上全部退透。结合图8(a)可以分析得到,在700K、800K、 900K条件下退火时,氧从晶体表面扩散到中心的时间分别为8h、5h和3h,而且退火温度越高,中心氧浓度的增大速度越快,这是由于晶体退火温度越高,氧在晶体内的扩散系数越大而引起的。从图8(b)中分析得到,要使晶体在700K条件下完成退火,退火时间至少在130h以上。

从图9中可以分析得知,退火温度对晶体表面A点的影响很小,都是氧在很短的时间内以大扩散通量从表面扩散进来,然后再以小扩散通量比较稳定地向晶体内部缓慢扩散。 只是温度越高,由A点向晶体扩散氧通量越大,并且在退火过程中扩散通量的下降速率也越大。这主要是因为温度越高,扩散系数越大,表明氧在高温条件下具有更好的扩散能力。对于B点而言,由于温度越高,扩散系数越大,氧能更快地从表面扩散到中心B点,浓度快速增大,降低了浓度差,使氧扩散驱动力降低,所以之后扩散通量又快速地减小而逐渐趋于稳定扩散。

3结论

(1)退火时,氧在短时间内以大通量从晶体表面扩散,然后以小通量比较稳定地向内部缓慢扩散,温度越高,氧扩散到晶体中心的时间越短,降低中心点的浓度梯度,使晶体内部扩散通量快速减小而逐渐趋于稳定扩散。

(2)退火前期氧迅速地扩散到大面积的中心低浓度区域,退火后期,晶体内部氧浓度逐渐增大,从表面往中心的氧扩散通量逐渐减小,而表面浓度急剧增大,扩散深度增大非常明显。

(3)在700K、800K、900K条件下退火50h后晶体内最小浓度分别为65%C0、83.9%C0和95.8%C0,氧从晶体表面扩散到中心的时间分别为8h、5h和3h,700K、800K退火深度分别为1mm、3mm,而900K时则退透。

摘要：针对钛酸锶晶体退火过程中氧原子的扩散问题建立了二维有限元模型,研究了温度、晶体尺寸、退火时间对氧原子扩散状态的影响规律。结果表明:氧在短时间内以大通量从晶体表面扩散,然后以小通量比较稳定地向内部缓慢扩散,温度越高,氧扩散到中心时间越短,在700K、800K、900K条件下退火50h,氧从晶体表面扩散到中心的时间分别为8h、5h和3h,退火深度分别为1mm、3mm和退透。

扩散过程 篇8

随着我国工业化的推进, 国内的环境污染问题越来越严峻。而当前最紧迫又很少受到其它环境科学涉及的是土壤和地下水污染的威胁[1]。随着地球物理技术的发展, 更多的地球物理方法被用于土壤和地下水污染的调查, 地球物理方法与地质方法相比具有成本低、速度快、原位无损等特点[2]。董路等[3]使用高密度电法对某场地内无机酸污染进行调查, 成功的圈定了被硫酸废液侵蚀的区域;程业勋等[4]对北京某垃圾填埋场的地下污染渗漏液进行综合地球物理方法调查, 较好的揭示了污染液体的渗漏趋势;王丽娟等[5]使用综合地球物理方法对某化工厂的氧化–蒸发塘的污水渗漏进行调查探测, 成功的查明渗漏点;杨进[6]等根据电阻率测深方法得到的视电阻率拟断面等值线图, 定性地反映了海水入侵的通道位置和断面形态。然而, 要加深对地下无机物污染液体渗漏过程的研究, 仅仅对地下无机污染液体进行静态的地球物理观测显得越来越不够, 要了解地下无机污染液体在土壤中的渗透过程, 就需要对无机污染液体在土壤中的渗透过程进行动态的监测。

延时性高密度电阻率法勘探的核心思想是通过对同一地点不同时间的地下介质的电阻率进行测量, 针对目的区域进行实时性跟踪勘探, 对比不同时间的地下介质的电阻率变化, 结合工程地质、环境分析以及其他物探资料, 判断引起电阻率变化的地质工程与环境意义, 从而确定真正的隐患部位。国外这方面的应用研究包括将延时性高密度电阻法应用于监测水在地下包气带中的渗流、地下水的抽取对地下水位的影响、化学污染物的溢出和堤坝的渗漏 (Johansson and Dahlin 1996) [7]。国内孙红亮等人利用延时性高密度电法研究评价水对某矿山排土场地质结构的影响, 通过对比雨季前和雨季后延时性高密度电法的成果资料, 说明了延时性高密度电法对于解决未饱和区域水的流动的可行性, 填补了高密度电法延时性勘探理论在国内的空白[8]。为了进一步研究和验证延时性高密度电法的勘探理论, 非常有必要进行理论研究和实验实践。在前期理论研究基础上, 本次在野外进行了延时性的高密度电法探索性实验。实验目的是详细研究温纳装置高密度电法对土壤浇盐水前、浇盐水后视电阻率的时空变化特征, 揭示盐水溶液在地下的渗透过程, 从而为研究无机物液体在土壤中的渗漏过程提供理论指导。

1 方法简介及物性基础

高密度电法是以地下被探测目标体与周围介质之间的电性差异为基础, 利用人工建立的稳定地下直流电场, 依据预先布置的若干道电极可灵活选定装置排列方式进行扫描观测, 研究地下一定范围内大量的丰富的空间电性特征, 从而查明和研究有关地质问题的一类直流电法勘探方法。

高密度电法的原理与普通电阻率法相同, 但是与普通电法相比, 高密度电法具有成本低、效率高、信息丰富、解释方便的特点, 因此被广泛的应用于水利水电工程、环境工程地质、工程地质勘查、城市工程勘察和工程质量检测等领域[9]。

高密度电法将多种常规直流电法的跑极方式集合起来, 工作中可根据实际需要选择不同装置采集数据。本次实验采用的是温纳装置, 温纳装置是一种四极装置, 在数据采集时有2个供电电极A、B和两个测量电极M、N, 4个电极的排列顺序依次为A、M、N、B, 且AM=MN=NB, 通过A、B两个电极向地下供固定电压或电流, M、N两极测量其两点之间的电流或者电压, 从而得到相应点的视电阻率值。

假设一条测线上所铺设的总电极数为n, 从第一个电极到第n个电极分别编号为1、2、3……n-1、n, 相邻2个电极之间的距离 (最小电极距) 为a。则在开始采集第1个数据时, A、M、N、B 4个电极分别对应1、2、3、4号电极, 即AM=MN=NB=a, 接下来4个电极依次同步的以距离a向后推移, 直到第n-3个数据A、M、N、B 4个电极分别对应n-3、n-2、n-1、n号电极, 这样所采集的n-3个数据深度一致, 它们为第1层的数据;接下来, A再回到1号电极点, 而使AM=MN=NB=2a, 即B、M、N分别对应3、5、7号电极, 测得第2层的第1个数据, 再同样的以距离a依次向后推移, 直到A、M、N、B 4个电极分别对应n-6、n-4、n-2、n, 得到了第2层的n-6个数据;再分别使得AM=MN=NB=3a、4a、5a……N·a (N为小于n/3的最大整数) , 即可得到1个倒梯形的全部N层的数据。工作中可以根据实际的需要选择数据采集到第几层。图1为温纳装置数据采集过程的示意图。

美国对2000多个垃圾填埋场渗透液进行分析, 发现44种无机物中普遍存在6种离子 (Cl-、Ca2+、Na+、Fe2+、NH4+和SO42-) [1]。因此, 为研究无机物污染液体在地下渗漏过程, 可用2.5%的氯化钠溶液模拟无机物污染液体。由于高离子浓度的液体具有良好的导电性, 这为高密度电法对地下的无机物污染液体渗漏过程进行监测提供了物性基础。

2 实验测试及分析

本次实验场地选择在武汉大学的一个草坪上, 地势较高, 土壤比较松散, 场地周围无明显水沟。

高密度电法采用重庆奔腾仪器厂的WDJD-3型高密度电阻率仪器, 60根铜电极, 最小电极距为0.3m, 剖面长度18m。测量采用温纳装置, 测量3层, 第1层对应深度 (采用AB/6估计) 约0.15m, 第2层对应深度为0.30m, 第3层对应深度为0.45m。

在距离测线约0.30m处, 挖了一条深约0.25m的水沟。该沟用于灌溉2.5%的盐水。

在灌溉盐水之前, 首先用温纳装置测量了土壤的视电阻率。然后开始灌溉盐水并计时, 以开始灌水的瞬间作为0时刻, 分别在0、6、12、20、30min时进行温纳装置的高密度电法观测, 每次温纳装置观测的时间约为2min左右, 在30min时刻停止灌溉盐水。然后分别在35、40、50、55、60、70、80、104、120min时进行了温纳装置的观测实验。图2是未加盐水之前的高密度电法剖面, 从图可以看出, 电阻率大致在30~60Ω·m, 表层电阻率较高, 在40~60Ω·m, 随着深度增大视电阻率逐渐减小, 说明表层土壤比较干燥, 深部含有一定水份电阻率略低。

图3是开始加盐水以后30min时的高密度电法剖面, 分析可知:

(1) 视电阻率在注水后0~30min内急剧下降, 30min内平均下降幅度大约为15Ω·m。在未加盐水时, 剖面内视电阻率集中在30~60Ω·m, 在加盐水后30min时刻, 整个剖面内视电阻率迅速下降, 低阻区域扩大, 剖面内视电阻率集中在15~45Ω·m, 剖面内出现了视电阻率最低值小于18Ω·m的区域。

(2) 从同一时刻剖面的视电阻率来看, 视电阻率仍然随着深度增大视电阻率逐渐减小。在加盐水后30min的时刻, 剖面表层的大部分测点视电阻率介于30~45Ω·m, 而表层以下的视电阻率大约为15~30Ω·m。

(3) 在注水后0~30min内, 随着时间的推移, 剖面逐渐显现出两个明显的低阻区域, 两个低阻中心的横坐标分别大约在x=3.9 m和x=11.85 m处。在未加盐水的时候, 在x=11.85m处的视电阻率比较明显的低于与其同深度不同横坐标的其它测点的视电阻率值, 其最低的视电阻率小于30Ω·m。注水30min后, 出现了另一低阻区域, 该区域中心的横坐标大约为3.9m, 区域内最小视电阻率小于22Ω·m, 而x=11.85 m处的低阻区的面积大幅扩大, 两低阻区内视电阻率小于26Ω·m的区域扩大到从中间相连, x=3.9m处的视电阻率最低值小于22Ω·m, x=11.85m处的的低阻区最小视电阻率降为小于18Ω·m。

图4是加盐水以后120min时刻的高密度电法剖面, 从图4分析可知:

(1) 视电阻率在急剧下降后开始趋于平稳, 并有恢复回升的趋势。例如, 在开始注水后30min时, 视电阻率小于22Ω·m的范围大约覆盖了3m<x<5m、-0.6m<z<-0.3m和10m<x<14m、-0.5<z<0m的两个区域, 而在开始注水120min后, 视电阻率小于22Ω·m的范围缩小至3.5m<x<4.5m、-0.6m<z<-0.3m和11m<x<13m、-0.4<z<0m。

(2) 剖面内两个低阻区域仍然存在, 并且很好的反映了整个剖面的视电阻率值随时间变化的趋势, 具有代表性。

图5是x=3.9m、AB/6=0.3m测点的视电阻率随时间变化曲线图。

从图5可以看出, 在开始注盐水0min时刻, x=3.9m、AB/6=0.3m测点处的视电阻率为43.377Ω·m, 在开始注盐水后6min时刻, 该测点处视电阻率急剧下降到约31.115Ω·m, 平均每分钟的下降值约达到2Ω·m, 在开始注盐水后12min时刻, 该测点处视电阻率下降到约27.867Ω·m, 平均每分钟下降值大约为0.53Ω·m。在20、30、35、40min时刻, 其视电阻率值分别大约为24.249、20.646、19.773、19.169Ω·m, 平均下降的速率估算分别为0.45、0.36、0.18、0.12Ω·m/min。这说明, 在开始注盐水后40min以内, 视电阻率急剧下降, 下降的速率逐渐降低。在开始注盐水后40~120min内, 该测点处的视电阻率基本稳定在20Ω·m左右, 在注水后50、55、60、70、80、104、120min时, 所测得该点的视电阻率分别为18.51、18.373、18.371、18.795、19.772、20.344、20.486Ω·m, 因此可以看出, 在注水后40~60min内, 其视电阻率仍在非常缓慢的下降, 在注水后60~120min内, 其视电阻率开始非常缓慢的上升。因此, 在经历了注水后前40min的视电阻率快速下降期以后, 视电阻率值基本稳定在20Ω·m左右, 这段稳定期内分别有一段视电阻率继续缓慢下降和一段视电阻率转而缓慢回升的过程。

图6的变化趋势与图5相似, x=11.85m、AB/6=0.15m测点处的视电阻率同样经历了快速下降期和稳定期。在t为0、30、120min时, 该测点相应的视电阻率分别为30.299、15.312、16.883Ω·m, 这是图6曲线上的3个极值。因此从图6可以看出, 其视电阻率快速下降期在注水后30min结束, 然后进入稳定期, 视电阻率在稳定期内缓慢的上升。

图7分别是注水后40min的视电阻率变化值剖面图, 其变化值为40min时剖面上所有测的视电阻率值减去注水前所测剖面上所有相应点的视电阻率值所得。

从图7可以看出, 在注盐水后40min时, 视电阻率减小量最大处位于x=3.9m、AB/6=0.3m和x=11.85m、AB/6=0.15m周围的两个低阻区域的中心部分, 其减小量最大值达到20~25Ω·m, 推测这两个低阻区域可能存在过水通道。

图8是注盐水后40min时的视电阻率变化值相对未加水时视电阻率值的比值剖面图, 其比值为注盐水后40min时剖面上所有点的视电阻率变化值与注盐水前所测剖面上所有对应点的视电阻率值的百分比值。

(△η= (ρ注水后-ρ注水前) /ρ注水前) 等值线剖面

从图8也可以看出, 注水后40min时视电阻率变化比值最大的区域与视电阻率剖面内两个低阻区域吻合, 该时刻视电阻率相对初始背景的视电阻率变化比值最大可达50%。

3 结论

通过本次实验研究表明, 高密度电法用于动态监测地下无机物溶液扩散的过程具有较好的效果。在地表浇注2.5%的盐水溶液, 可以导致浅部土壤的视电阻率下降15~25Ω·m, 电阻率相对变化值最高可高达50%。地表浇注盐水溶液后, 地下浅部土壤的视电阻率会随着土壤含盐水量的变化而变化, 视电阻率的整个变化过程可以分为快速下降期和平稳恢复期两个阶段。由此可见, 延时性高密度电法对盐水溶液的扩散有很好的探测效果, 揭示了延时性高密度电法这种新的地球物理方法可有效检测无机类化学污染物的扩散与渗漏过程。

参考文献

[1]程业勋, 杨进, 赵章元.环境地球物理学的现状与发展[J].地球物理学进展, 2007, 22 (4) :1364~1369.

[2]Weller A., Franges W., Seichter M..Three-dimensionalinversion of induced polarization data from simulated waste[J].Journal of Applied Geophysics, 2000, 40 (2~3) :67~83.

[3]董路, 叶腾飞, 能昌信等.ERT技术在无机酸污染场地调查中的应用[J].环境科学研究, 2008, 21 (6) :67~71.

[4]程业勋, 刘海生, 赵章元.城市垃圾污染的地球物理调查[J].工程地球物理学报, 2004, 1 (1) :26~30.

[5]王丽娟, 王支农, 郝书军.地球物理方法在环保工作中的应用实例[J].勘察科学技术, 2003, (2) :58~63.

[6]杨进, 刘庆成, 程业勋等.环境地球物理方法在地下水污染监测中的应用[J].环境科学研究, 1998, 11 (6) :43~46.

[7]M.H.Loke.Electrical imaging surveys for environmental andengineering studies[M].Colorado:Colorado U.S.Geological Survey, 1999, 31~34.

[8]孙红亮, 肖宏跃, 雷宛.高密度电法的延时性勘探研究[J].成都理工大学学报 (自然科学版) , 2008, 35 (2) :158~161.

扩散过程 篇9

关键词：Fe-Cr合金,扩散,分子动力学,五频率模型

铁素体/马氏体 (Ferritic/Martensitic, F/M) 钢具有优异的几何稳定性、抗辐照肿胀和优良的耐腐蚀性能等特性, 可以作为未来聚变示范堆和第四代核反应堆候选结构材料之一[1]。研究表明, Cr的加入对材料的辐照性能产生了很大的影响[2]。F/M钢辐照过程中会产生富Cr的α’相, 导致材料脆化、应力腐蚀和化学腐蚀等性能变化[3]。目前α’相的形核过程的动力学细节还不清楚, 但是肯定应与空位扩散机制有关。核反应堆中材料受到中子辐照, 产生大量空位。因此要想弄懂形核机理, 必须了解合金中原子扩散过程。目前国内外关于Cr在FeCr合金中扩散的研究尚未完善, 该文采用计算机模拟方法来研究Cr在Fe-Cr合金中扩散过程的原子尺度的相关问题。首先通过原子扩散的五频率模型, 计算了Cr在Fe-Cr稀合金中的扩散系数, 再通过分子动力学模拟计算得到了Cr在Fe-Cr合金中的扩散系数, 对比两者的结果, 验证计算的准确性, 从而为相场模拟、动力学蒙特卡罗模拟等研究工作提供关于扩散的基础数据。

1 势函数

常用的合金体系势函数模型是嵌入势EAM (embedded-atom method) [4]。对于FeCr体系, Olsson等人发展了双带模型 (two band model, 2BM) 势, 该模型在EAM势基础上加入了s电子的贡献, 体系中原子i能量表示为:

V (rij) 是原子在i和j位置电子的相互作用。

代表d和s电子的密度, Fb为与原子间距离有关的函数。F (ρ) 可表示为

系数Ai是参数化的结合力和排斥力的相对强度。

该文采用Olsson等给出的2BM势函数进行计算。利用Lammps程序, 建立一个包含1024个Fe原子的bcc-Fe模型, 用能量最小化方法计算得到bcc-Fe的晶格常数a和总能量E, 并通过 (4) 式和 (5) 式分别计算了Fe中的空位形成能Efv和Fe-Cr合金中Cr与空位V之间的结合能:

N是原子总数, V代表空位。 (4) 式中E[ (N-1) Fe+1V]表示包含N-1个Fe和1个空位的总能量, E[Fe]表示纯bcc-Fe体系中Fe的能量, 类推可知 (5) 式中各项含义。对 (5) 式计算时, 分别使Cr和V位于彼此的第一近邻和第二近邻位置, 可相应得到Cr原子与空位的第一和第二近邻结合能bE1, nnCr-V和bEC, 2nnr-V, 结果如表1所示。

2 Cr原子扩散系数的MD模拟

该文选取Fe-1at%Cr合金体系, 进行分子动力学模拟计算。所建立MD模型与3.2节中所述的相同, Cr的扩散系数的计算方法同上。对含有杂质的bcc-Fe来说, 空位平衡浓度Cve (T) 的计算式为:

通过拟合公式 (22) , 可得出Cr原子在Fe-Cr合金系中的扩散系数表达式为:

对比式 (25) 与 (27) 可看出, 五频率模型的计算结果与分子动力学模拟计算的结果比较相符, 说明五频率模型的计算结果较为准确。

3 结语

该文通过原子模拟方法和2B M势函数, 计算了Cr在Fe-Cr合金中扩散系数。首先通过五频率模型预测Cr原子的扩散系数, 另外用MD模拟直接计算Cr原子扩散系数。研究表明, 2BM势适用于Cr原子在FeCr合金中扩散的研究, 分子动力学和五频率模型两种方法的计算结果比较吻合, 说明五频率模型计算结果较准确, 可用于研究杂质在合金中的扩散。

参考文献

[1]贺新福, 杨文樊, 胜论.Fe Cr合金幅照损伤的多尺度模拟[J].物理学报, 2009, 58 (12) :8657-8669.

[2]Miller M K, Hyde J M, Cerezo A, et al.Comparison of low temperature decomposition in Fe Cr and duplex stainless steels[J].Applied surface science, 1995, 87:323-328.

[3]Bonny G, Terentyev D, Malerba L.Identification and characterization of Cr-rich precipitates in Fe Cr alloys:An atomistic study[J].Computational Materials Science, 2008, 42 (1) :107-112.