Disentanglement by Dissipative Open System Dynamics
非线性动力学入门-西安交通大学教师个人主页
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
另一方面梁的轴向应变的表达式也会因变形大小的不同而采用不同的表达式比如小变形时应变而当考虑大变形时可能采用的应变表达式就是进而得到的梁的振动方程将会是一个含有高度非线性项的偏微分方程组
非线性动力学入门
张新华
西安交通大学 工程力学系 2011 年 07 月
前 言
─1687 年，牛顿(Isaac Newton, 1643 ～ 1727)发表了《自然哲学之数学原 理》(Mathematical Principles of Natural Philosophy)，标志着经典力学(亦即牛 顿力学)的正式诞生。牛顿力学主要研究自由质点系的宏观运动规律。 ─1788 年，拉格朗日(Joseph Louis Lagrange, 1736 ～ 1813)发表了分析力 学教程(Analytical Mechanics)，标志着拉格朗日力学的诞生。Lagrange 力学属 于分析力学的主要内容之一，在位形空间中研究带有约束的质点系动力学。 ─1833 年，哈密尔顿(William Rowen Hamilton, 1805 ～ 1865)对 Lagrange 力学进行了改造，引进了相空间(2n 维空间)，对系统内在的对称性(辛对称， Symplectic)进行了刻画。狭义上的哈密尔顿力学只适用于保守系统，而广义 的哈密尔顿力学在适用于非保守系统。哈密尔顿力学也属于分析力学的主要 组成部分。在此后发展起来的量子力学中 Hamilton 力学发挥着巨大的作用。 目前在天体力学、计算 Hamilton 力学，量子力学，甚至弹性力学(即所谓的 辛弹性力学)中哈密尔顿力学依然发挥着重要作用。 ─1927 年，Birkhoff(George David Birkhoff, 1844 ～ 1944)发表了“动力系 统”(Dynamical Systems)，标志着 Birkhoff 动力学的正式问世。Birkhoff 动力 学建立了研究非完整力学的框架。 ─1892 ～ 1899， 彭加莱(Henri Poincaré, 1854 ～ 1912)发表了三卷本的“天 体力学中的新方法”(New Methods of Celestial Mechanics)，系统性地提出了 研究动力学系统的定性方法，即几何方法。经典力学的目标之一就是设法求 得系统的解析解，而 Poincaré意识到对于大多数非线性系统而言，求其解析 解是不可能的，而必须发展新的研究方法。他超越了他的时代，极富远见地 预测到了非线性系统混沌现象(系统的解对初始条件具有极端敏感依赖性)的 存在。更为重要的是，Poincaré开创了研究非线性动力系统的几何方法，当之 无愧地被誉为非线性科学之父，其影响是划时代的。 ─1892 年，李亚普诺夫(Aleksandr Mikhailovich Lyapunov, 1857 ～ 1918)在 他的博士论文“运动稳定性的一般问题”(General problem of the stability of motion )中，系统地探讨了非线性动力学系统的稳定性问题。他提出了两种研 究稳定性的方法：李亚普诺夫第一方法(间接方法)和李亚普诺夫第二方法(直 接方法)。他从代数角度出发，对动力学系统的研究开创了一个崭新的领域。 彭加莱与李亚普诺夫，前者从几何角度，后者从代数角度，开拓了非线 性科学的研究疆域和研究手段。 ─1963 年，Lorenz(Edward Norton Lorenz, 1917 ～ 2008)发表了“确定性 非周期流”(Deterministic Nonperiodic Flow)的论文，认为大气系统的性态对 初值极为敏感，从而导致准确的长期天气预报是不可能的。该文标志着人类 首次借助于计算机发现了混沌(Chaos)现象的存在。 ─1757 年，欧拉(Leonhard Euler, 1707 ～ 1783)发表了压杆稳定性的论 文，首次探讨了力学系统的分岔现象。作为分岔理论重要分支的突变理 论(Catastrophe Theory)则主要由法国数学家托姆(René Thom, 1923 ～ 2002)于 上个世纪 60 年代创立，由齐曼(Christopher Zeeman，1925 ～)在 70 年代大力 推广普及。 ─1834 年，英国的罗素(John Scott Russell, 1808 ～ 1882)骑着马在 Union 运河上散步时，发现了现在称之为孤立波(又称作孤波，Solitary wave)的 i
时间尺度上非迁移Birkhoff系统的Mei对称性定理
时间尺度上非迁移Birkhoff系统的Mei对称性定理*张毅†(苏州科技大学土木工程学院, 苏州 215011)(2021 年2 月25日收到; 2021 年9 月9日收到修改稿)研究并证明时间尺度上非迁移Birkhoff系统的Mei对称性定理. 首先, 建立任意时间尺度上Pfaff-Birkhoff原理和广义Pfaff-Birkhoff原理, 由此导出时间尺度上非迁移Birkhoff系统(包括自由Birkhoff系统、广义Birkhoff系统和约束Birkhoff系统)的动力学方程. 其次, 基于非迁移Birkhoff方程中的动力学函数经历变换后仍满足原方程的不变性, 给出了时间尺度上Mei对称性的定义, 导出了相应的判据方程. 再次, 建立并证明了时间尺度上非迁移Birkhoff系统的Mei对称性定理, 得到了时间尺度上Birkhoff系统的Mei守恒量.并通过3个算例说明了结果的应用.关键词：Birkhoff系统, Mei对称性定理, 时间尺度, 非迁移变分学PACS：45.20.Jj, 11.30.Na, 02.30.Xx DOI: 10.7498/aps.70.202103721 引 言Birkhoff力学起源于Birkhoff[1]的著作《动力系统》. Santilli[2]首次提出Birkhoff力学一词, 并详细地讨论了Birkhoff方程的构造、变换理论及其对强子物理的应用. 梅凤翔等[3]和Galiullin等[4]从各自角度分别独立地研究了Birkhoff系统动力学, 他们的研究各具特色且更侧重于分析力学. 文献[5]构建了广义Birkhoff系统动力学. 梅凤翔先生[6]指出Birkhoff力学是分析力学发展的第4个阶段. 近年来, Birkhoff力学在对称性理论[7−13]、几何动力学[14,15]、全局分析与稳定性[16,17]、数值计算[18−22]等研究方向上都取得了重要进展.时间尺度, 即实数集的任意非空闭子集, 最早是由Hilger博士[23]引进的. 由于实数集和整数集本身就是一类特殊的时间尺度, 因而在时间尺度上不仅可以统一地处理连续系统和离散系统, 而且可以处理既有连续又有离散的复杂动力学过程. 近20年来, 时间尺度分析理论不仅在理论上不断完善[24−26], 其应用领域也在不断拓广[27−34]. 文献[35]最早提出并研究了时间尺度上基于delta导数的自由Birkhoff系统动力学及其Noether对称性. 文献[36]利用对偶原理将文献[35]的结果拓展到nabla导数情形. 文献[37]给出了时间尺度上非迁移Birkhoff系统的Noether定理. 但是, 这些研究尚限于: 1)自由Birkhoff系统; 2) Noether对称性;3)守恒量是Noether型的. 文献[38, 39]初步研究了时间尺度上Birkhoff系统的Lie对称性和Mei 对称性, 但是其守恒量的证明基于第二Euler-La-grange方程, 而数值计算表明该方程并不成立[34].此外, 根据Bourdin[33]的研究, 在离散层面非迁移情形的结果是保变分结构及其相关性质的, 尽管迄今时间尺度上非迁移变分问题研究还很少. 本文研* 国家自然科学基金(批准号: 11972241, 11572212)和江苏省自然科学基金(批准号: BK20191454)资助的课题.† 通信作者. E-mail: zhy@© 2021 中国物理学会 Chinese Physical Society 究时间尺度上非迁移Birkhoff 系统的Mei 对称性,包括自由Birkhoff 系统、广义Birkhoff 系统和约束Birkhoff 系统, 建立并证明上述3类Birkhoff 系统的Mei 对称性定理, 给出时间尺度上新型守恒量, 称之为Mei 守恒量.2 时间尺度上非迁移Birkhoff 方程关于时间尺度上微积分及其基本性质, 读者可参阅文献[24, 25].2.1 Pfaff-Birkhoff 原理及其推广在时间尺度上, 非迁移Pfaff 作用量为R β:T ×R 2n →R B :T ×R 2n →R a ∆βa βC 1,∆rd(T )β,γ=1,2,···,2n a γσρ其中 是时间尺度上Birkhoff 函数组, 是时间尺度上Birkhoff 函数, 是Birkhoff 变量 对时间的delta 导数. 设所有函数都是 函数. .非迁移是指作用量(1)中的变量 没有经过前跳算子 或后跳算子 的作用而发生跃迁[33].等时变分原理且满足端点条件以及互易关系原理(2)称为时间尺度上非迁移Pfaff-Birkhoff 原理.等时变分原理(2)可推广为Φβ=Φβ(t,a γ)式中 表示附加项[5]. 原理(5)式可称为时间尺度上非迁移广义Pfaff-Birkhoff 原理.2.2 自由Birkhoff 系统由原理(2), 容易导出σ(t )δa ∆β其中 是前跳算子. 考虑到 的独立性,由时间尺度上Dubois-Reymond 引理[24], 得到C β其中 为常数. 因此有方程(8)为时间尺度上非迁移Birkhoff 方程.2.3 广义Birkhoff 系统由原理(5), 可导出类似于方程(8), 有方程(10)可称为时间尺度上非迁移广义Birkhoff-方程.2.4 约束Birkhoff 系统约束方程为将(11)式取变分, 得由(6)式和(12)式, 容易导出λj =λj (t,a β)λj 其中 为约束乘子. 假设约束(11)式相互独立, 则由(11)式和(13)式可解出 . 于是方程(13)可写成P β=λj∂f j∂a β其中 . 方程(14)可视作与约束Birk-hoff 系统(13)和(11)相应的自由Birkhoff 系统.只要初始条件满足约束方程(11), 那么方程(14)的解就给出约束Birkhoff 系统的运动.3 Mei 对称性3.1 自由Birkhoff 系统引进无限小变换t →ϑ(t )=t +υξ0+o (υ)C 1,∆rd υ∈R ϑ(t )¯T¯σ¯∆其中映射 是1个严格递增 函数, 是无限小参数, 是一个新的时间尺度 , 前跳算子为 , delta 导数为 .B R β¯B ¯Rβ在变换(15)下, 动力学函数 和 变换为 和 , 有υ=0将(16)式在 处Taylor 级数展开, 得到Y (0)=ξ0∂/∂t +ξβ∂/∂a β其中 .定义1 对于时间尺度上非迁移Birkhoff 系统(8), 如果成立, 则变换(15)称为Mei 对称性的.判据1 如果变换(15)满足如下判据方程:则变换相应于时间尺度上非迁移Birkhoff 系统(8)的Mei 对称性.3.2 广义Birkhoff 系统B R βΦβ¯B¯R β¯Φβ设时间尺度上动力学函数 , 和 经历变换(15)后, 成为 , 和 , 有于是有下述定义2和判据2.定义2 对于时间尺度上非迁移广义Birkhoff 系统(10), 如果成立, 则变换(15)称为Mei 对称性的.判据2 如果变换(15)满足如下判据方程:则变换相应于时间尺度上非迁移广义Birkhoff 系统(10)的Mei 对称性.3.3 约束Birkhoff 系统B R βP βf j ¯B ¯R β¯P β¯f j 设时间尺度上动力学函数 , 和 , 以及约束 经历变换(15)后, 成为 , , 和 , 有于是有下述定义3和判据3.定义3 对于时间尺度上与约束Birkhoff 系统(13)和(11)相应的自由Birkhoff 系统(14), 如果成立, 则变换(15)称为Mei 对称性的.判据3 如果变换(15)满足如下判据方程:则变换相应于时间尺度上相应自由Birkhoff 系统(14)的Mei 对称性.定义4 对于时间尺度上约束Birkhoff 系统(13)和(11), 如果方程(24)以及如下方程成立, 则变换(15)称为Mei 对称性的.判据4 如果变换(15)满足判据方程(25)和如下限制方程:则变换相应于时间尺度上约束Birkhoff 系统(13)和(11)的Mei 对称性.4 Mei 对称性定理4.1 自由Birkhoff 系统定理1 假设变换(15)满足判据方程(19), 则时间尺度上非迁移Birkhoff 系统(8)存在新型守恒量G M 其中 是规范函数, 满足因此, (28)式是系统的守恒量. 证毕.定理1可称为时间尺度上非迁移Birkhoff 系统(8)的Mei 对称性定理, (28)式称为Mei 守恒量.4.2 广义Birkhoff 系统定理2 假设变换(15)满足判据方程(22), 则时间尺度上非迁移广义Birkhoff 系统(10)存在新型守恒量G M 其中 是规范函数, 满足证明∇∇tI M =0将方程(22)和方程(33)代入(34)式, 得到, 于是(32)式是系统的守恒量.定理2可称为时间尺度上非迁移广义Birkhoff 系统(10)的Mei 对称性定理, (32)式称为Mei 守恒量. 证毕.4.3 约束Birkhoff 系统定理3 假设变换(15)满足判据方程(25), 则时间尺度上与约束Birkhoff 系统(13)和(11)相应的自由Birkhoff 系统(14)存在新型守恒量G M 其中 是规范函数, 满足G M 定理4 假设变换(15)满足判据方程(25)和限制条件(27)式, 则时间尺度上约束Birkhoff 系统(13)和(11)存在新型守恒量(35), 其中规范函数 满足方程(36).定理3为时间尺度上与约束Birkhoff 系统(13)和(11)相应的自由Birkhoff 系统(14)的Mei 对称性定理.定理4为时间尺度上非迁移约束Birkhoff 系统的Mei 对称性定理, (35)式是Mei 守恒量.5 算 例例1 研究时间尺度上Birkhoff 系统, 设Birk -hoff 函数和Birkhoff 函数组为试研究该系统的Mei 对称性与守恒量.由方程(8)得到T =R 如取 , 则方程(38)成为这是著名的Hojman-Urrutia 问题[3,4]. 该问题本质上不是自伴随的, 因此没有Lagrange 结构或Hami-lton 结构.下面来计算Mei 对称性. 经计算, 有取生成函数为则生成函数(41)满足判据方程(19), 因此它相应于系统的Mei 对称性. 将(41)式代入方程(29), 可解得由定理1, 系统有Mei 守恒量, 形如(44)式表明, 对于任意的时间尺度, (44)式都是Birkhoff 系统(37)的守恒量. 如取生成函数为那么生成函数(45)也是Mei 对称的, 由方程(29)得由定理1, 得到Mei 守恒量T =R σ(t )=t 对于守恒量(47), 如果系统是通常的Birkhoff 系统, 即取 , 则 , 从而(47)式给出T =h Z h>0σ(t )=t +h 这是通常意义下Hojman-Urrutia 问题的守恒量[3].如果是离散情形, 即取 , 这里 , 则 , 从而(47)式成为h 这是步长为 的离散版本的Mei 守恒量.例2 研究时间尺度上广义Birkhoff 系统的Mei 对称性与守恒量.广义Birkhoff 方程(10)给出计算Mei 对称性, 由于将(52)式代入判据方程(22), 有解(53)式和(54)式相应于系统的Mei 对称性. 将(53)式代入方程(33), 解得由定理2, 系统有Mei 守恒量, 形如G M =−2t 同理, 相应于生成函数(54), , 由定理2得(56)式和(57)式是由Mei 对称性(53)和(54)导致的Mei 守恒量.例3 研究时间尺度上约束Birkhoff 系统约束为g φ试研究其Mei 对称性与守恒量,其中 和 是常数.方程(13)给出由方程(59)和方程(60),解得因此有做计算取生成函数为则µ(t )=σ(t )−t ν(t )=t −ρ(t )其中 为向前互差函数, 为向后互差函数. 由判据4, 生成函数(64)相应于系统的Mei 对称性. 将(65)式代入方程(36),解得由定理4, 系统有Mei 守恒量, 形如6 讨 论T =R σ(t )=t µ(t )=0如果取时间尺度 , 则前跳算子 ,互差函数 , 因此上述结果退化为通常意义下Birkhoff 系统、广义Birkhoff 系统和约束Birkh-off 系统连续版本的变分原理、Birkhoff 方程和Mei 对称性定理.T =R 例如, 对于自由Birkhoff 系统, 当取 时,原理(2)成为方程(8)成为由判据方程(19)容易得到于是, 定理1退化为下述推论1.推论1 假设变换(15)满足判据方程(19),则自由Birkhoff 系统(69)的Mei 对称性导致如下G M 其中 是规范函数, 满足推论1是通常意义下自由Birkhoff 系统连续版本的Mei 对称性与守恒量定理[7]. 而方程(68)、方程(69)和方程(71)就是通常意义下自由Birk-hoff 系统连续版本的Pfaff-Birkhoff 原理、Birk-hoff 方程和Mei 守恒量.T =h Z h >0σ(t )=t +h µ(t )=h 如果取时间尺度 , 常数 , 则前跳算子 , 互差函数 . 此时, 原理(2)成为方程(8)成为则定理1退化为下述推论2.推论2 假设变换(15)满足判据方程(19), 则自由Birkhoff 系统(74)的Mei 对称性导致如下形式的守恒量:G M 其中 是规范函数, 满足h 推论2是通常意义下自由Birkhoff 系统离散版本的Mei 对称性与守恒量定理. 而方程(73)—(75)就是通常意义下自由Birkhoff 系统离散版本步长为 的Pfaff-Birkhoff 原理、Birkhoff 方程和Mei 守恒量.7 结 论对称性和守恒量一直是分析力学研究的一个重要方面. 经典的对称性主要有Noether 对称性和Lie对称性. Mei对称性是本质上不同于前两种对称性的一种不变性, 它可以导致Mei守恒量. Mei守恒量不同于Noether守恒量, 是一种新的守恒量. 本文提出并研究了时间尺度上非迁移Birkhoff系统的Mei对称性定理.一是提出了时间尺度上非迁移Pfaff-Birkhoff 原理和广义Pfaff-Birkhoff原理, 导出了时间尺度上Birkhoff系统, 包括自由Birkhoff系统、广义Birkhoff系统和约束Birkhoff系统的动力学方程.主要结果是原理(2)和(5), Birkhoff方程(8), (10)和(13).二是定义了时间尺度上非迁移Birkhoff系统的Mei对称性, 并导出了它的判据方程. 主要结果是4个定义和4个判据.三是提出并证明了时间尺度上非迁移Birkhoff 系统、非迁移广义Birkhoff系统和非迁移约束Birkhoff系统的Mei对称性定理. 主要结果是4个定理, Mei守恒量(28), (32)和(35).T=R T=h Z当取时间尺度和时, 文中定理给出通常意义下自由Birkhoff系统、广义Birkhoff 系统和约束Birkhoff系统的连续版本和离散版本的Mei对称性与守恒量定理. 由于除了实数集和整数集以外, 时间尺度还可以有很多选择, 因此时间尺度上Birkhoff系统的Mei对称性定理具有一般性.参考文献B irkhoff G D 1927 Dynamical Systems (Providence: AMSCollege Publ. ) pp59–96[1]S antilli R M 1983 Foundations of Theoretical Mechani cs II (New York: Springer-Verlag) pp1–280[2]M ei F X, Shi R C, Zhang Y F, Wu H B 1996 Dynamics of Birkhoffian System (Beijing: Beijing Institute of Technology Press) pp1–228[3]G aliullin A S, Gafarov G G, Malaishka R P, Khwan A M1997 Analytical Dynamics of Helmholtz, Birkhoff and Nambu Systems (Moscow: UFN) pp118–226[4]M ei F X 2013 Dynamics of Generalized Birkhoffian Systems (Beijing: Science Press) pp1–206[5]M ei F X, Wu H B, Li Y M, Chen X W 2016 J. Theor. Appl.Mech. 48 263 (in Chinese) [梅凤翔, 吴惠彬, 李彦敏, 陈向炜2016 力学学报 48 263][6]M ei F X 2004 Symmetries and Conserved Quantities of [7]Constrained Mechanical Systems (Beijing: Beijing Institute of Technology Press) pp1–482 (in Chinese) [梅凤翔 2004 约束力学系统的对称性与守恒量 (北京: 北京理工大学出版社) 第1—482页]W ang P, Xue Y, Liu Y L 2012 Chin. Phys. B 21 070203 [8]Z hang Y, Zhai X H 2015 Nonlinear Dyn. 81 469[9]Z hang H B, Chen H B 2017 J. Math. Anal. Appl. 456 1442 [10]Z hang Y 2018 Int. J. Non-Linear Mech. 101 36[11]X u X X, Zhang Y 2020 Acta Phys. Sin. 69 220401 (in Chinese) [徐鑫鑫, 张毅 2020 物理学报 69 220401][12]Z hang L J, Zhang Y 2020 Commun. Nonlinear Sci. Numer.Simul. 91 105435[13]G uo Y X, Liu C, Liu S X 2010 Commun. Math. 18 21[14]L iu S X, Liu C, Guo Y X 2011 Chin. Phys. B 20 034501[15]C hen X W, Li Y M 2013 Acta Mech. 224 1593[16]L uo S K, He J M, Xu Y L 2016 Int. J. Non-Linear Mech. 78 105[17]L iu S X, Liu C, Guo Y X 2011 Acta Phys. Sin. 60 064501 (in Chinese) [刘世兴, 刘畅, 郭永新 2011 物理学报 60 064501] [18]L iu S X, Liu C, Hua W, Guo Y X 2016 Chin. Phys. B 25 114501[19]K ong X L, Wu H B, Mei F X 2013 Appl. Math. Comput. 225 326[20]K ong X L, Wu H B 2017 Acta Phys. Sin. 66 084501 (in Chinese) [孔新雷, 吴惠彬 2017 物理学报 66 084501][21]H e L, Wu H B, Mei F X 2017 Nonlinear Dyn. 87 2325[22]H ilger S 1990 Results Math. 18 18[23]B ohner M, Peterson A 2001 Dynamic Equations on TimeScales (Boston: Birkhäuser) pp1–353[24]B ohner M, Georgiev S G 2016 Multivariable DynamicCalculus on Time Scales (Switzerland: Springer International Publishing AG) pp1–600[25]G eorgiev S G 2018 Fractional Dynamic Calculus andFractional Dynamic Equations on Time Scales (Switzerland: Springer International Publishing AG) pp1–357[26]A tici F M, Biles D C, Lebedinsky A 2006 Math. Comput.Modell. 43 718[27]B ohner M 2004 Dyn. Syst. Appl. 13 339[28]B artosiewicz Z, Torres D F M 2008 J. Math. Anal. Appl. 3421220[29]B enkhettou N, Brito da Cruz A M C, Torres D F M 2015Signal Process. 107 230[30]D ryl M, Torres D F M 2017 Springer Proceedings inMathematics & Statistics 195 223[31]H an Z L, Sun S R 2014 Oscillation Theory of DynamicEquations on Time Scales (Jinan: Shandong University Press) pp1–232 (in Chinese) [韩振来, 孙书荣 2014 时间尺度上动态方程振动理论 (济南: 山东大学出版社) 第1—232页][32]B ourdin L 2014 J. Math. Anal. Appl. 411 543[33]A nerot B, Cresson J, Belgacem K H, Pierret F 2020 J. Math.Phys. 61 113502[34]S ong C J, Zhang Y 2015 J. Math. Phys. 56 102701[35]S ong C J, Zhang Y 2017 J. Nonlinear Sci. Appl. 10 2268 [36]S ong C J, Cheng Y 2020 Appl. Math. Comput. 374 125086 [37]Z hang Y 2019 Chaos, Solitons Fractals 128 306[38]Z hang Y, Zhai X H 2019 Commun. Nonlinear Sci. Numer.Simul. 75 251[39]Mei’s symmetry theorems for non-migrated Birkhoffiansystems on a time scale*Zhang Yi †(College of Civil Engineering, Suzhou University of Science and Technology, Suzhou 215011, China)( Received 25 February 2021; revised manuscript received 9 September 2021 )AbstractThe Mei symmetry and its corresponding conserved quantities for non-migrated Birkhoffian systems on a time scale are proposed and studied. Firstly, the dynamic equations of non-migrated Birkhoffian systems (including free Birkhoffian systems, generalized Birkhoffian systems and constrained Birkhoffian systems) on a time scale are derived based on the time-scale Pfaff-Birkhoff principle and time-scale generalized Birkhoff principle. Secondly, based on the fact that the dynamical functions in the non-migrated Birkhoff’s equations still satisfy the original equations after they have been transformed, the definitions of Mei symmetry on an arbitrary time scale are given, and the corresponding criterion equations are derived. Thirdly, Mei’s symmetry theorems for non-migrated Birkhoffian systems on a time scales are established and proved, and Mei conserved quantities of Birkhoffian systems on a time scale are obtained. The results are illustrated by three examples.Keywords: Birkhoffian system, Mei’s symmetry theorem, time scale, non-migrated variational calculus PACS: 45.20.Jj, 11.30.Na, 02.30.Xx DOI: 10.7498/aps.70.20210372* Project supported by the National Natural Science Foundation of China (Grant Nos. 11972241, 11572212) and the Natural Science Foundation of Jiangsu Province, China (Grant No. BK20191454).† Corresponding author. E-mail: zhy@。
TKX-50热分解氮气形成机理的分子动力学模拟
文章编号：1006-9941 (2018)01-0075-05T K X-50热分解氮气形成机理的分子动力学模拟余一\张蕾M,姜胜利\王星\赵寒月\陈军1>2(1.中国工程物理研究院高性能数值模拟软件中心，北京100088; 2.北京应用物理与计算数学研究所，北京100088)摘要：为研究新型富氮含能化合物5,5^联四唑-1，r-二氧二羟铵（T K X-50)高能钝感背后的微观机制，采用从头算分子动力学方法模拟了T K X-50在不同压力及温度下的分解过程，通过分析主要产物N，的生成路径，揭示了T K X-50热分解随温度与压力变化的规律。
模拟显示T K X-50分解的主要产物为叱0和其中叱存在三条主要的生成路径，两条来源于唑环环裂过程，另一条与铵盐和唑环的相互作用相关联。
唑环环裂直接生成叱的过程受温度影响较大，温度越高，断裂速度越快，对压力不敏感。
铵盐与唑环相互作用生成\的过程则依赖于扩散，扩散速率与温度呈正相关，与压力呈负相关。
三条反应路径的共同作用使得TKX-50的反应速率宏观上呈现随温度升高而升高，随压力升高而下降的趋势。
关键词：5,5，-联四唑-1,1，-二氧二羟铵（T K X-50);热分解；反应路径；分子动力学模拟中图分类号：TJ55;0643.12文献标志码：A D O I：10.11943/j.is s n.1006-9941 .2018.01 .0091引言基于碳骨架氧化放能的设计思路，高能炸药的边界已经从黑索今（RDX)、奥克托今（HM X)等传统炸药拓展到了八硝基立方烷（0NC)[1]，六硝基六氮杂异伍兹烷（CL-20)[2]等新型的高能化合物。
这些新型化合物充分利用了笼状结构中蕴藏的张力，使之分解时释放出更多的能量，从而提升爆炸威力。
由于这些新型化合物普遍具有敏感度高，合成步骤复杂，成本高等缺点，所以人们继而将目光转向了高氮材料，由于N帒N三键的键能远高于单键或双键，这些材料在分解形成氮气时能放出大量能量，又因为产物N2对环境友好，故高氮材料被认为是含能材料设计与合成的未来[3]。
原子分子光子系统的耗散相互作用和退相干
原子与分子物理学报JOURNAL OF ATOMIC AND MOLECULAR PHYSICS第37卷第6期2020年12月Vol. 37 No. 6Dec. 2020doi : 10.19855/j.l000-0364.2020.065001原子分子光子系统的耗散相互作用和退相干景俊(浙江大学物理学系光学研究所#杭州310027$摘要：原子分子系统与量子化的电磁场或光子模式耦合的系统是非相对论量子力学理论研究和实验研 究的主要对象和模型•现实系统必然与外界环境耦合,且即便原子隔绝较好、光学腔壁品因子足够高,原子系统也不等价于少数几个能级构成的简单模型：它仍然有不为零的几率跃迁到不可控的能级空间、与原 子相互作用的自由空间真空场的量子效应也必须考虑.本文将结合开放量子系统理论的基本要素与与子光学学基本模型，对原子分子系统在电磁场中发生的耗散以及量子子相干过程做简单综述，并重点介绍描 述量子系统退相干过程简主流理论工具—— 主方程.关键词：原子子学；量子子相干；量子子方程中图分类号：O65 文献标识码：A 文章编号：1000-0364(2020)06-0935-05Relaxation and decoherence process of the atom - photon systemJING Jun(Institute of Optics , Department of Physics , Zhejiang University , Hangzhou 310027 # China )Abstract : Atomic and molecular systems coupled te the quantized electro 一 magnetic field os photons constitute the most part of the interested systems os models in the field of non - relativistio quantum mechanics. Even fot the welS 一 isolated dtomio system as well as the high Q - factoo covity , the quantum mechanico. system connot be considered as a simpSe modd consisted of sevvraS leve's. In practico , C is affected by the surrounding environ ment. Fot instances , the quantum system has nonzero probability te jump out of the subspaco undet contco and it is subjeC te tie quantum effecO raised by the interaction with tie vvcuum state of the free spaco. In the frame work of open 一 quantum 一 system dynamico combining with the fundamentae modes of atomio optico, we wiH briefly review the relevvnt basio concepts of the dissipation and decoherenco processes of O v atomio systems cou pled t the exOrna felds. The focus of the paper is to introduco the quantum master equation as the mainstream tod tr theof quantum decoherenco.Key wores : Atomio optico ; Quantum decoherenco ; Quantum mastee equation1引言原子分子系统与量子化的电磁场或光子模式耦合的系统⑴无论在理论还是实验方面都是非相 对论量子力学研究的主要对象和模型[2-4].如果将系统拓展到人工原子、超导回路、离子阱、光子晶格、量子点等新量子平台.i],那么它不仅涵盖了传统的原子分子物理和量子光学的主要内 容，而且构成量子信息和量子调控的物理基础[12,13].正如在单个二能级原子和单模激光场的耦合模型 Jaynes - Cummings 模型中所观察到的，如果不考虑能量耗散和量子退相干效应，那收稿日期：2020-5-19基金项目：国家自然科学基金( 11974311, U1801661 )；浙江省自然科学基金(LD18A040001 )；中央高校基础研究基金(2018QNA3004)作者简介：景俊(1979—),男，江苏省泰州市人，博士，主要致力于开放量子系统动力学、非微扰量子调控、几何量子计算、量子退相干、量子光学系统、固态量子调控、以及量子力学基本理论方面的研究.E-maf : **************.cn936原子与分子物理学报第37卷么它的动力学是相干且可逆的•这正是薛定定方程所描述的结果•对人工原子，包括量子点、超导量子比特、里德堡原子，以及新颖电磁辐射源，包括激光、微波激射、同步辐射、微波源的广泛研究涵盖了从无线电波到远红外的频段与原子分子系统的耦合，从而在相当程度上再次引起人们对原子与光相互作用的兴趣.新的研究方法能够获得有关原子和分子结构以及动力学的更确切的信息、从而有助于控制原子和分子的内部和外部自由度、并且产生新的辐射形式.这些进展已经引来越来越多的物理学家、化学家、其他研究人员和工程师对物质和低能辐射之间发生的作用［14'17］产生兴趣（本文的目标是为已经具有一定量子力学和原子分子物理基础的研究人员介绍在开放量子力学框架下研究原子与光相互作用的基础理论，特别是关于量子系统退相干动力学的描述理论一主方程.主方程及其等价方法构成了整合涉及物质—辐射相互作用的几乎所有物理方面的基本工具［18-23］2开放量子系统与量子退相干在量子力学基本问题方面，含时薛定量方程以及冯诺伊曼方程已经在原则上提供了闭合量子系统的动力学方程.如果不计较开放系统的环境或者噪声对它的影响和作用，或者拥有对全部希尔伯特空间的计算能力，就没有必要在开放量子系统的框架下处理问题•但实际上对于大多数实际问题这两者都是无法做到的•首先定义何为环境.在全量子框架下，环境是与量子系统相互作用的另一个量子系统，拥有足够多、足够复杂的自由度，以至于相对人们感兴趣的量子系统而言，是不受控制的.具有可操作性的处理方案的前提是通过测量或物理建模给出环境自由度的统计性质.这些性质包括且不限于玻色统计或费米统计、谱密度函数、关联函数、温度一如果它处在热力学平衡态下的话.然后通过理论推演，使得这些统计性质反映在受研究的量子系统的动力学方程中•其次何为噪声.布朗运动中的涨落可逆当作是一种经典噪声，按照涨落—耗散定理的要求，它具有一定的关联函数和谱密度函数，它们两者互为傅里叶变换关系.量子系统面对的噪声，可以是经典噪声，也可可是量子噪声.在物理上，两者都要符合涨落—耗散定理；在特定情况下,符合一定关联函数的经典噪声和量子噪声对系统的影响在系综统计后不可区分.因此有时可可用噪声取代环境自由度的全量子化处理.在纯量子框架下，噪声可区直接用环境算符以及新合系数直接表达出来.就其效应而言，有时需要特别指出，哪些噪声造成了系统的退极化和能量耗散,哪些噪声不涉及能量损耗而直接反映出同一系综内不同量子系统在相位上的脱散.前者反映了原子分子系统的布居数在高低能态上的转移；后者代表了相位关联的消失.对量子系统而言，前者当然也同时弓来系统间失去相位上的固有联系，但后者在许多文献中才被称为“纯量子的退相干”.因此量子退相干本质上是一种量子力学效应，与经典耗散或随机涨落有所区别，虽然后者也会引起量子效应.在哥本哈根学派定义的量子测量假设以及后来由冯诺伊曼严格化后的量子测量理论中，不难发现退相干是量子测量的效应.测量仪器扮演了量子化的环境自由度的角色，系统根据与待测物理量配套的测量仪器指针态展开的本征基矢做随机投影.无论是读出结果的测量，还是不读结果的测量，指针态对应的本征基矢之间的相干项在测量后全部消失整在广泛关注量子信息处理的前沿研究中，退相干显然是任何量子处理，特别是量子计算，或任何必须利用量子线性叠加原理优势的技术不可回避的问题.它在绝大多数时候显然是难以逾越的障碍.根据DiVincenzo为可行的量子计算机定出的必须满足的若干条件，量子系统的相干时间除以量子门的运行时间必须足够大，否则量子线路模型就没有实际价值.所以无论是量子物理的基本问题，还是量子前沿技术的发展，都需要人们重视退相干过程的基础研究，从而设计出提高量子相干时间的方案,或者要么回避、要么正面利用退相干效应.这里不对有到已放3主方程的推导方法现在在入技术性环节一如何建立开放的原子分子系统的退相干动力学主方程•“主”这个词意思就是关注某一少自由度或少体系统，而把其环境或不关心的希尔伯特子空间对该系统的作用考虑到方程的结构和参量中去，最后得到方程仅含有少自由度或少体系统的算符.假设初始时刻系统与环境是没有有合的，即*（o）=*（0）第6期景俊：原子分子光子系统的耗散相互作用和退相干9379*(+),且整体哈密顿量可写为：*S*9D+D9*+(*(1)其中*是系统的自由哈密顿量，*是环境的自由哈密顿量，h,S&>+9G是系统与环境相互作用哈密顿量，这里的+和G分别表示系统和环境的算符，s>代表第第种量子退相干通道.(用来衡量相互作用哈密顿量相对系统自由哈密顿量的相对幅度，一般符合微扰论的要求.为简单起见，以下在相互作用表象下进行推导，即对所有算符o，都有0(:=9h+h”0(2)首先从整体系统的冯诺伊曼方程出发(它总是成立的):%P t(：dt二一*(：,pT((3)并将其形式解*(:二P t(+)-加+s[H(s)，p T(7/代入方程的右边对易子.可以得到：%pT(：dt二-ia[*(：,p(0)/-+s.h(：,[h(s)，P t(s)]](4)然后对该方程做第一步马尔可夫近似以及玻恩近似，即将⑷式要p(s)替换为p(t)，物理上抛弃了从初始时刻到当下时刻之前的系统演化历史对系统当前变化量的贡献，而是只关心当下时刻的贡献，数学上带来的误差为0((3).并让方程等号两边都对环境部分自由度进行部分求迹，从而得到关于系统密度算符的二阶微扰的动力学方程(注意在下面这个方程后，不再使用约等号)：df p t2—Jc.h(：，p j(0)/-(2JdsJc.h(：，[h(s)，p T(：//(5) 0其中第一项正比于Tc.G(t)p E(0)/，也就是环境算符G的期望值，可可将其设为零.这本质上不是一个假设，因为总可下令H"&，+9 (G-2G〉)、H"H+(&〈G3s>.在环境弛豫时间远小于系统动力学时间的假设下，可可认为环境状态基本不变(一般可认为环境处在热力学平衡态)，这就得到了微扰论下的玻恩-马尔可夫主方程(Bom-Markov master equation):—^%(—二-J%T C.h(：，[H(s)，p s(t)9pg(0)//(6)它在二阶微扰成立的前提下可可包含量子环境的结构参量，也就是各模式与系统的相对该合强度、各模式自身的本征频率的统计信息，因而在许多文献中也被称为“非马尔可夫”主方程.注意它不是“精确主方程”.23/，因为没有包含结构化的环境对系统动力学的所有阶次的贡献.从微扰论物理的角度需要进一步明确两点：1，主方程的建立基于系统-环境相互作用哈密顿量与系统自由哈密顿量之间存在能量/时间的尺度分离,否则对全部希尔伯特空间的计算就不可避免从2,环境的状态一般默认为热力学平衡态，那么系统在环境噪声下的随机量子跃迁过程和耗散应该受到涨落-耗散定理的限制.4真空环境下的Lindblad主方程下下用相对具体的一个例子进一步推演，为简单起见，将(吸收进方程.假设相互作用哈密顿量H,=SB e+S e B.系统跳跃算符满足S(t): e iHst Se~iHst=Se-，0t,其中,是原子分子系统本征频率.环境自由哈密顿量为h#=&宀；e；,环境集体湮灭算符为B=&冲；.这里；[是玻色环境的湮灭算符，g[是是合系统与第第个模式的贡合系数.所以在相互作用表象下，H(t)=SB e(t)+S e B(t),B(t):&t g[e-(心0)t(7)于是玻恩-马尔可夫主方程可具体表达为P s(t)=-F(t)[S s p S(t)-S p S(t)S]-'[p(t)SS e-S e p s(t)S/+h.C.(8)其中两个含时系数的定义是tF(j)二+乃#.b(：#(s)P b],t'=\dsTi B.B#(s)B(t)*#/(9)这里暗含环境状态p为热力学平衡态的假设.另外对于许多本征频率为微波波段的人工原子而言，光频波段的电磁环境因为能量与之远不匹配而总可以认为处在真空或零温.所下可合理假定电磁环境处在真空态或零温，则有938原子与分子物理学报第37卷F(-=\&」g k*9,[_,0)Q-s)\d,](,)(%-9s(,_，o)0几(10)其中=(,)=&」*+(,-,[)是环境态密度或模式密度函数.而'(-=0.略加计算，即将F(-的积分上限推至无穷，消除其含时特征，也就是执行第二步马尔可夫近似，就可以得到相互作用表象下的Lindblad型主方程：5[S t S,P s(-/(11)其中退相干因子％=,)8(,-，o)与真空自发辐射理论中的费米黄金法则吻合；而5=-!>+，=，是著名的兰姆位移，反应了原0，-，0子分子系统因为受到电磁场环境的作用而产生的能级偏移.注意教科书里的兰姆位移一般发生在基态，但实际上它也可夫发生在激发态.5讨论需要再次强调的是，主方程是量子力学框架下的动力学方程一不难发现它在幺正变换下保持形式不变.一个合理的主方程必然得到合法的系统密度矩阵，也就是每一时刻的系统密度矩阵必然是半正定的.动力学方程所对应的系统状态的映射必须是全正定且迹守恒(Complete Positive Trace-preserving).要得到特定物理体系的主方程，原则上应从特定全量子模型出发，仿照前文步步推演，不可夫做过分的非物理的近似•比如将玻恩-马尔可夫主方程的积分上限推至无穷有时候在短时动力学中会带来非正定的问题，其物理原因就是这种近似违背了物理定律的因果律.发Lindblad方程为代表的现象学方程，比如量子光学领域熟悉的光学布洛赫方程，只要保证系数为正常数，就不会带来非物理的结果•但是它不能反映强耦合情况下的物理，比如非马尔可夫效应.另外跳跃算符应用系统哈密顿算符的本征态基矢展开，从而在物理上对应真实的量子跳跃(quantum jump)过程.只有本征态基矢之间的瞬时量子跳跃才能正确地描述微观量子过程•值得指出的是，只有极其少数特定主方程可以得到解析解.比如纯退相位(pure dephasing)模型和零温环境下的少数量子跳跃模型.有时不得不采用量子跳跃方法(quantum jump method)、量子轨迹方法"quantum trajectoa approach)或量子郎之万方程替代量子Lindblad方程进行数值求解.比如可以证明：发下的量子态扩散方程(一种量子轨迹方程或随机薛定谔方程).*4，*5]与Lindblad 方程完全等价：*-1(z”)〉=(—*+Sz：-1S e S)I1(z”)〉(12)其中：是满足系综平均条件A[::]=%(t-s)的高斯型白噪声，这里町]是指对系综求平均值.求解随机薛定定方程后还要对噪声轨迹求系综平均，从而得到密度矩阵*(-=A[I1(:*)3〈1(:*)I]-6总结和展望量子退相干理论在当前研究中的兴起和应用伴随着单量子系统的实验实现.通常布洛赫方程或者爱因斯坦几率方程被用来描述受到电磁场驱动的原子分子系综的动力学•但离子阱、单模腔、量子点、超导回路等新兴量子平台的发明提供了观察和操控单个粒子的可能性.对辐射场信息的获取，比如通过对系统发射到环境内光子的探测,会导致系统波函数的突变•无论是探测到光子还是没有探测到，也就是零结果，都会导致信息的增加.这样新的洞察便会注入到原子光子动力学及耗散过程中，这就发展出新的退相干动力学方法.除了对物理的新洞察，这些方法也会使得对复杂问题的模拟成为可能，比如用主方程方法也可处理激光冷却问题•参考文献：[1]Cohen-Tannoudji C，Dupont一Roc J，Grynbere G.Atom,一photon interactions:Basic processes and applications[M].New York:Wiley，1992.[2]Cohen-Tannoudji C.Manipulating atoms with photons[J].Rev.Mod.Phys.，1998,70：707.[3]Chu S.The manipulation of neutral particles[J].Rev Mod.Phys.，1998，70:685.[4]Phillips W ser cooling and trapping of neutral atoms[J].Rev.Mod.Phys.，1998，70:721.[5]CywiskiL，Witzel W M，Das Sarma S.Pure quantumdephasing of a solid-state electron spin qubit in a第6期景俊：原子分子光子系统的耗散相互作用和退相干939large nuclear spin bath coupled by long-range hyperfine-mediated interactions.J].Phys.Rev.B#2009,79：245314.[6]You J Q,Nori F.Atomio physics and quantum opticsusing superconducting circuits[J].Nature,2011,474：589.[7]Ashhab S.Superradiance transition in a system with asingle qubit end a single osciiatoc[J].Phys.Rev.A,2013,87：013826.[8]Jing J,Wu L-A.Decoherencc and control of a qubitin spin baths:an exact master equation study[J].Sci.Rep.,2018,8：1471.[9]Ferraro E,Breuee H P,Napoli A,Jivulescu M A,Messina A.Non-Markovian dynamics of a single e-lechon spin coupled to a nucleac spin bath[J].Phys.Rp.B,2008,78：064309.[10]Leibfried D,Blat t R,Monroe C,Wineland D.Quantum dynamics of single trapped ions[J/.Rp.Mod.Phys.,2003,75：281.[11]Beaudoin F,Gambett o J M,Blais A.Dissipation andultrastrong coupling in circuit QED[J/.Phys.Rp.A,2011,84：043832.[12]Nielsen M A,Chuang I.Quantum,computatioo andquantum information[M/.England:Cambridge Uni-i>asiiyPa s,2000.[13]Wiseman H M,Milburn G J.Quantum measurepentand control[M/.Cambridge:Cambridge UniversityPa s,2009.[14]Breuer H P,Petruccione F.The pheory of opee quantum wystem.[M/.Oxford:Oxford University Press,2002.[15]Gardiner C,doOer P.Quantum nofe[M/.Berlin:Springer-Verlaa,2004.[16]Lexget t A J,Chakravarty S,Dorsey A T,et al.Dynamics of the dissipative two-state system[J/.Rp.Mrd Phy.,1987,59： 1.[17/Breuer H P,Kappler B,Petruccione F.Stochastic wave-function method for non-Markovian quantummaster equations[J/.Phy.Rp.A,1999,59:1633.[18]Lindblad G.On the generators of quantum dynamicaisemigroups[ J/.Commun.Matt.Phys.,1976,48:119.[19]Strunz W T,Diosi L,Gisin N.Open system dynamicswith non-Markovian quantum trajectories[J/.Phys.Rp.L^p.,1999,82：1801.[20]Zurek W H.Decoherencc,einselection,and the quantum origins of the cOssicai[J/.Rp.Mof.Phy.,2003,75:715.[21/Kofman A G,Kurizki G.Unified theory of dynamicaay suppressed qubit decoherencc in thermal baths[J/.Phy.Rp.Lett.,2004,93:130406.[22]Eree N,Gordon G,Nest M,Kurizki G.Thermody-诅]!^^conWoi by Frequent quantum measurements[J].Nature,2008,452：724.[23]Vacchini B,Breuer H P.Exact master equations forthe non-Markovian deesy of a qubit[J/.Phy.Rp.A,2010,81:042103.[24]Jing J,Zhao X,You J Q,et al Time-local quantum-state-diffusion equation for mululevei quantum systems[J].Phy.Rp.A,2012,85:042106.[25]Zhang J,Jing J,Wu L A,et al Measurement-induced cooeingooaqubitin steuctueed eneieonments[J] (Phy.Rp.A,2019,100:022107(。
PhysRevA.87.042115
反应扩散系统中的图灵斑图动力学介绍：———非线性科学专题之十
反应扩散系统中的图灵斑图动力学介绍：———非线性科学
专题之十
欧阳颀
【期刊名称】《物理通报》
【年(卷),期】1999(000)005
【摘要】1 图灵斑图 1952年,被后人称为计算机科学之父的著名英国数学家阿兰·图灵(Alan Turing)把他的目光转向生物学领域。
在他的一篇著名论文《形态形成的化学基础》中,图灵用一个反应扩散模型成功地说明了某些生物体表面所显示的图纹,如斑马身上的斑图是怎样产生的。
想象在生物胚胎发育的某个阶段,生物体内某些被称为'形态子'的生物大分子与其他反应物发生生物化学反应,同时在机体内随机扩散。
图灵表明在适当的条件下,这些原来浓度均匀分布的'形态子'会在空间自发的组织成一些周期性的结构。
【总页数】4页(P4-7)
【作者】欧阳颀
【作者单位】北京大学非线性科学中心
【正文语种】中文
【中图分类】O414
【相关文献】
1.图灵斑图动力学的数学机制 [J], 刘迎东
2.具有非线性发生率的传染病模型图灵斑图研究 [J], 王涛;靳祯;孙桂全
3.斑图动力学：非线性科学专题之九 [J], 夏蒙棼
4.螺旋波的斑图动力学：非线性科学专题之十一 [J], 欧阳颀
5.双层耦合非对称反应扩散系统中的振荡图灵斑图 [J], 刘雅慧;董梦菲;刘富成;田淼;王硕;范伟丽
因版权原因，仅展示原文概要，查看原文内容请购买。
量子动力学中的非马可夫过程及其应用
量子动力学中的非马可夫过程及其应用量子动力学是研究微观粒子在量子力学框架下的运动规律的学科，它对于理解和描述微观世界的行为具有重要意义。
在传统的马可夫过程中，系统的演化是由一个马尔可夫矩阵所决定的，而在量子动力学中，存在着一类非马可夫过程，其演化规律不再遵循马尔可夫性质。
本文将介绍量子动力学中的非马可夫过程及其应用。
非马可夫过程是指系统的演化不仅与当前状态有关，还与之前的历史状态有关。
在量子动力学中，非马可夫过程的出现主要是由于量子纠缠效应的存在。
量子纠缠是指两个或多个微观粒子之间存在一种特殊的关联关系，当其中一个粒子发生测量时，另一个粒子的状态也会瞬间发生变化。
这种关联关系使得系统的演化不再遵循马尔可夫性质，而是呈现出非局域性的行为。
非马可夫过程在量子信息科学中具有重要的应用价值。
量子信息科学是研究利用量子力学原理进行信息处理和传输的学科，它在密码学、通信、计算等领域具有巨大的潜力。
在量子计算中，非马可夫过程可以用来实现量子比特之间的相互耦合，从而实现量子门操作，进而实现量子计算的高效性。
在量子通信中，非马可夫过程可以用来实现量子纠缠的分发和保护，从而实现安全的量子通信。
除了在量子信息科学中的应用，非马可夫过程还在其他领域展现出了巨大的潜力。
在凝聚态物理中，非马可夫过程可以用来描述电子在凝聚态材料中的输运行为。
传统的马可夫过程无法很好地描述电子在强关联体系中的行为，而非马可夫过程可以更加准确地描述这种行为。
在生物物理学中，非马可夫过程可以用来描述生物分子的运动行为。
生物分子的运动往往受到周围环境的影响，而非马可夫过程可以更好地模拟这种复杂的运动行为。
在实际应用中，非马可夫过程的研究面临着许多挑战。
首先，由于非马可夫过程的演化规律与系统的历史状态有关，因此需要对系统的演化进行全程的追踪和记录。
这对于实验的设计和数据处理提出了更高的要求。
其次，非马可夫过程的数学描述相对复杂，需要借助于量子力学的数学工具和方法进行分析和求解。
sic模具高温模压石英玻璃物相接触角的分子动力学模拟
第３９卷第３期２０２０年３月硅㊀酸㊀盐㊀通㊀报ＢＵＬＬＥＴＩＮＯＦＴＨＥＣＨＩＮＥＳＥＣＥＲＡＭＩＣＳＯＣＩＥＴＹＶｏｌ.３９㊀Ｎｏ.３Ｍａｒｃｈꎬ２０２０ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟吴㊀悠１ꎬ２ꎬ３ꎬ邹㊀斌１ꎬ２ꎬ３ꎬ王俊成１ꎬ２ꎬ３ꎬ黄传真１ꎬ２ꎬ３ꎬ朱洪涛１ꎬ２ꎬ３ꎬ姚㊀鹏１ꎬ２ꎬ３(１.山东大学机械工程学院先进射流工程技术研究中心ꎬ济南㊀２５００６１ꎻ２.山东大学高效洁净机械制造教育部重点实验室ꎬ济南㊀２５００６１ꎻ３.山东大学机械工程国家级实验教学示范中心ꎬ济南㊀２５００６１)摘要:针对高温下石英玻璃纳米液滴在ＳｉＣ模具表面接触角难以测量的问题ꎬ采用分子动力学方法ꎬ模拟研究了不同温度和粗糙表面面向模压的ＳｉＯ２/ＳｉＣ高温接触角以及ＳｉＯ２熔体的界面结构ꎮ应用压力张量法发现了ＭＳ￣Ｑ势函数模拟的ＳｉＯ２熔体表面张力较接近实际值ꎬ即ＳｉＯ２高温表面特性模拟可优先采用ＭＳ￣Ｑ势函数ꎮ针对ＳｉＣ模具纳米级表面的粗糙度ꎬ发现当粗糙度因子ｒ>１.５时润湿模式由Ｗｅｎｚｅｌ变为Ｃａｓｓｉｅ￣Ｂａｘｔｅｒꎬ此时Ｒａ的变化对接触角值无明显影响ꎬＲｍｒ值减小使得接触面积分数ｆ减小ꎬ接触角值随之增大ꎮ因此ꎬ保持ｒ大于１.５的同时适当减小Ｒｍｒ值有利于减小固液摩擦ꎬ降低石英玻璃工件和ＳｉＣ模具界面上的脱模力ꎮ随着温度升高ＳｉＯ２表面结构变得松散ꎬ导致其在ＳｉＣ表面接触角减小ꎮ在超过２３００Ｋ时接触角值的变化率增大ꎬ为减小工件￣模具界面的粘附ꎬ模压温度应选择２３００Ｋ以下ꎮ关键词:石英玻璃ꎻ碳化硅ꎻ接触角ꎻ分子动力学ꎻ模压中图分类号:Ｏ６４７.１ꎻＴＳ９４３.６５㊀㊀文献标识码:Ａ㊀㊀文章编号:１００１￣１６２５(２０２０)０３￣０９２３￣０９ＭｏｌｅｃｕｌａｒＤｙｎａｍｉｃｓＳｉｍｕｌａｔｉｏｎｏｆＨｉｇｈＴｅｍｐｅｒａｔｕｒｅＣｏｎｔａｃｔＡｎｇｌｅｆｏｒＭｏｌｄｅｄＳｉｌｉｃａＧｌａｓｓｉｎＳｉＣＤｉｅＷＵＹｏｕ１ꎬ２ꎬ３ꎬＺＯＵＢｉｎ１ꎬ２ꎬ３ꎬＷＡＮＧＪｕｎｃｈｅｎｇ１ꎬ２ꎬ３ꎬＨＵＡＮＧＣｈｕａｎｚｈｅｎ１ꎬ２ꎬ３ꎬＺＨＵＨｏｎｇｔａｏ１ꎬ２ꎬ３ꎬＹＡＯＰｅｎｇ１ꎬ２ꎬ３(１.ＣｅｎｔｅｒｆｏｒＡｄｖａｎｃｅｄＪｅｔＥｎｇｉｎｅｅｒｉｎｇＴｅｃｈｎｏｌｏｇｉｅｓꎬＳｃｈｏｏｌｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇꎬＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙꎬＪｉｎａｎ２５００６１ꎬＣｈｉｎａꎻ２.ＫｅｙＬａｂｏｒａｔｏｒｙｏｆＨｉｇｈＥｆｆｉｃｉｅｎｃｙａｎｄＣｌｅａｎＭｅｃｈａｎｉｃａｌＭａｎｕｆａｃｔｕｒｅꎬＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙꎬＭｉｎｉｓｔｒｙｏｆＥｄｕｃａｔｉｏｎꎬＪｉｎａｎ２５００６１ꎬＣｈｉｎａꎻ３.ＮａｔｉｏｎａｌＤｅｍｏｎｓｔｒａｔｉｏｎＣｅｎｔｅｒｆｏｒＥｘｐｅｒｉｍｅｎｔａｌＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇＥｄｕｃａｔｉｏｎꎬＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙꎬＪｉｎａｎ２５００６１ꎬＣｈｉｎａ)Ａｂｓｔｒａｃｔ:Ｄｕｅｔｏｔｈｅｄｉｆｆｉｃｕｌｔｉｅｓｏｆｍｅａｓｕｒｅｍｅｎｔｏｆｃｏｎｔａｃｔａｎｇｌｅｗｈｅｎｐｒｅｓｓｉｎｇｓｉｌｉｃａｇｌａｓｓｎａｎｏ￣ｄｒｏｐｌｅｔｉｎＳｉＣｄｉｅｓｕｒｆａｃｅａｔｈｉｇｈｔｅｍｐｅｒａｔｕｒｅꎬｔｈｅＳｉＯ２/ＳｉＣｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｃｏｎｔａｃｔａｎｇｌｅａｔｄｉｆｆｅｒｅｎｔｔｅｍｐｅｒａｔｕｒｅａｎｄｒｏｕｇｈｓｕｒｆａｃｅａｓｗｅｌｌａｓｔｈｅｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆＳｉＯ２ｍｅｌｔｗｅｒｅｓｉｍｕｌａｔｅｄｂｙｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓ.ＵｓｉｎｇｔｈｅｐｒｅｓｓｕｒｅｔｅｎｓｏｒｍｅｔｈｏｄꎬｉｔｉｓｆｏｕｎｄｔｈａｔｔｈｅｓｕｒｆａｃｅｔｅｎｓｉｏｎｏｆＳｉＯ２ｍｅｌｔｓｉｍｕｌａｔｅｄｂｙＭＳ￣Ｑｐｏｔｅｎｔｉａｌｆｕｎｃｔｉｏｎｉｓｃｌｏｓｅｔｏｔｈｅｅｘｐｅｒｉｍｅｎｔａｌｖａｌｕｅ.ＴｈｕｓｔｈｅＭＳ￣ＱｐｏｔｅｎｔｉａｌｃａｎｂｅｐｒｅｆｅｒｒｅｄｉｎｓｉｍｕｌａｔｉｏｎｏｆＳｉＯ２ｈｉｇｈ￣ｔｅｍｐｅｒａｔｕｒｅｓｕｒｆａｃｅｃｈａｒａｃｔｅｒ.ＦｏｒｔｈｅｒｏｕｇｈｎｅｓｓｏｆｎａｎｏｓｃａｌｅｓｕｒｆａｃｅｏｆＳｉＣｄｉｅ.Ｗｈｅｎｔｈｅｒｏｕｇｈｎｅｓｓｆａｃｔｏｒｒｅｘｃｅｅｄｓ１.５ꎬｔｈｅｗｅｔｔｉｎｇｍｏｄｅｃｈａｎｇｅｓｆｒｏｍＷｅｎｚｅｌｔｏＣａｓｓｉｅ￣Ｂａｘｔｅｒ.ＡｔｔｈｉｓｔｉｍｅꎬｔｈｅｃｈａｎｇｅｏｆＲａｈａｓｎｏｓｉｇｎｉｆｉｃａｎｔｉｍｐａｃｔｏｎｔｈｅｃｏｎｔａｃｔａｎｇｌｅꎬａｎｄｔｈｅｒｅｄｕｃｅｄＲｍｒｖａｌｕｅｒｅｓｕｌｔｓｉｎｔｈｅｄｅｃｒｅａｓｅｏｆｃｏｎｔａｃｔａｒｅａｆｒａｃｔｉｏｎｆａｎｄｔｈｅｉｎｃｒｅａｓｅｏｆｃｏｎｔａｃｔａｎｇｌｅ.Ｔｈｅｒｅｆｏｒｅꎬｋｅｅｐｉｎｇｒｇｒｅａｔｅｒｔｈａｎ１.５ｗｈｉｌｅａｐｐｒｏｐｒｉａｔｅｌｙｒｅｄｕｃｉｎｇＲｍｒｖａｌｕｅｉｓｃｏｎｄｕｃｉｖｅｔｏｒｅｄｕｃｉｎｇｓｏｌｉｄ￣ｌｉｑｕｉｄｆｒｉｃｔｉｏｎａｎｄｒｅｄｕｃｉｎｇｔｈｅｓｔｒｉｐｐｉｎｇｆｏｒｃｅｏｎｔｈｅｉｎｔｅｒｆａｃｅｂｅｔｗｅｅｎｆｕｓｅｄｓｉｌｉｃａｗｏｒｋｐｉｅｃｅａｎｄＳｉＣｄｉｅ.ＡｓｔｈｅｔｅｍｐｅｒａｔｕｒｅｉｎｃｒｅａｓｅｓꎬｔｈｅｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆＳｉＯ２ｂｅｃｏｍｅｌｏｏｓｅꎬｒｅｓｕｌｔｉｎｇｉｎｉｔｓｃｏｎｔａｃｔａｎｇｌｅｄｅｃｌｉｎｉｎｇｏｎｔｈｅＳｉＣｓｕｒｆａｃｅ.Ｔｈｅｃｈａｎｇｅｒａｔｅｏｆｔｈｅｃｏｎｔａｃｔａｎｇｌｅｉｎｃｒｅａｓｅｓｗｈｅｎｔｅｍｐｅｒａｔｕｒｅｅｘｃｅｅｄｓ２３００Ｋ.Ｉｎｏｒｄｅｒｔｏｒｅｄｕｃｅｔｈｅａｄｈｅｓｉｏｎｏｆｔｈｅｗｏｒｋｐｉｅｃｅｗｉｔｈｔｈｅｄｉｅꎬｔｈｅｍｏｌｄｉｎｇｔｅｍｐｅｒａｔｕｒｅｓｈｏｕｌｄｂｅｂｅｌｏｗ２３００Ｋ.Ｋｅｙｗｏｒｄｓ:ｓｉｌｉｃａｇｌａｓｓꎻｓｉｌｉｃｏｎｃａｒｂｉｄｅꎻｃｏｎｔａｃｔａｎｇｌｅꎻｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓꎻｃｏｍｐｒｅｓｓｉｏｎｍｏｌｄｉｎｇ基金项目:山东省自然科学基金重大基础研究项目(ＺＲ２０１８ＺＢ０５２１)作者简介:吴㊀悠(１９９４￣)ꎬ男ꎬ硕士研究生ꎮ主要从事模具摩擦磨损方面的研究ꎮＥ￣ｍａｉｌ:ｊｏｌｍｕｇｉ＠１６３.ｃｏｍ通讯作者:邹㊀斌ꎬ教授ꎮＥ￣ｍａｉｌ:ｚｂ７８＠ｓｄｕ.ｅｄｕ.ｃｎ９２４㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷０㊀引㊀言近年来ꎬ微结构光学元件因其在信息通信领域的广泛应用受到越来越多的关注ꎮ模压是制造光学玻璃元件常用的方法之一ꎮ模压过程中ꎬ一方面模具上精准的几何结构复印到了玻璃元件上ꎻ另一方面玻璃元件与模压模具之间的摩擦与粘附造成了模具的磨损ꎬ从而导致模具几何结构精度的丧失ꎮ这些均为模具与工件的固液界面相互接触㊁相对运动作用的结果ꎮ因此ꎬ有必要研究它们的接触特性ꎬ从而对模具结构和模压工艺的设计提供一定的指导ꎮ润湿是能够实现模压的先决条件ꎬ而接触角是表征润湿性的主要参数ꎮ对于理想光滑表面的接触角模型ꎬＹｏｕｎｇ[１]从固㊁液㊁气界面张力平衡的角度建立了经典的杨氏方程ꎮ对于粗糙表面和非均质表面ꎬＷｅｎｚｅｌ[２]和Ｃａｓｓｉｅ[３]等学者在杨氏方程的基础上分别建立了Ｗｅｎｚｅｌ和Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ模型ꎮＷｅｎｚｅｌ润湿模式中ꎬ液体始终完全浸渍粗糙结构波谷中ꎬ即出现所谓 钉扎 现象ꎮ而在模压工艺中ꎬ这一现象增大了工件液相在模具表面流动的阻力ꎬ从而增大了固液摩擦ꎬ使得模具更易磨损ꎮ而Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿模式认为由于粗糙结构波谷中存在气相ꎬ液体无法浸渍其中ꎬ因此是模压比较期望的润湿模式ꎮ在光学玻璃材料中ꎬ石英玻璃的光学性能有其独特之处ꎮ它既可以透过远紫外光谱ꎬ且透射率在同类材料中最优ꎬ又可透过可见光和近红外光谱ꎮ同时ꎬ其机械性能也高于普通玻璃[４]ꎮ目前ꎬ石英玻璃光学元件的制造多采用超精密磨削或是激光刻蚀的方法ꎬ生产效率低ꎬ生产成本高ꎮ而模压方法可以提高加工效率ꎬ实现大规模生产ꎮ然而ꎬ由于石英玻璃软化点温度高达１５００~１６００ħꎬ常用的模具材料如ＷＣ等在此温度下无法正常工作ꎬ而ＳｉＣ热稳定性好ꎬ熔点达２８３０ħ[５]ꎬ故较为适合作为模压石英玻璃的模具材料ꎮ由于模压温度高ꎬ超出了目前接触角测量仪器的工作范围ꎬ并且微结构模具表面粗糙度达到纳米级ꎬ要研究粗糙度对润湿性的影响需要测量纳米尺度微液滴的接触角ꎬ因此试验测定比较困难ꎮ分子动力学方法使微纳液滴高温接触角的模拟及研究其成因机理成为可能ꎮ徐威等[６]通过改变ＬＪ势函数中的作用参数模拟了纳米水滴在不同能量表面上的铺展过程和润湿形态ꎬ模拟结果与经典润湿理论计算得到的结果呈现相似变化趋势ꎮ王龙[７]研究了铜和金液滴在石墨烯和碳纳米管等不同结构基底上的润湿和融合过程ꎬ发现金属液滴的融合受液滴形状和基底结构影响ꎮ由于势函数的限制ꎬ目前对于界面润湿的分子动力学模拟多集中于ＬＪ流体ꎬ如水和液氩等ꎬ而多元素流体较少见于报道ꎮ势函数按多体作用的复杂程度可分为对势和多体势ꎮ对于石英玻璃的模拟ꎬ多体势不能较好地计算Ｓｉ￣Ｏ键的断裂ꎬ因此不利于高温动力学性能的研究[８]ꎮＢｅｅｓｔ等[９]基于石英玻璃经典的ＢＭＨ势提出了ＢＫＳ势ꎬ并运用第一性原理方法结合实验数据拟合了参数ꎮＳｕｎｄａｒａｒａｍａｎ等[１０]使用ＢＫＳ势预测了石英玻璃的力学性能ꎬ发现当短程和长程截止半径分别为５.５Å和１０Å时模拟的石英玻璃结构更符合实际ꎮＤｅｍｉｒａｌｐ等[１１]首先将Ｍｏｒｓｅ势与电荷平衡法(ＱＥｑ)相结合ꎬ建立了ＭＳ￣Ｑ力场ꎬ并研究了石英玻璃在压力变化过程中的相变ꎮ丁元法等[１２]比较了ＢＫＳ与ＭＳ￣Ｑ模型下石英玻璃的高温扩散特性ꎬ认为计算高温下石英玻璃的扩散传输性能可优先选择ＭＳ￣Ｑ力场ꎬ但并未比较其表面特性ꎮ从以上文献可以看出ꎬＢＫＳ和ＭＳ￣Ｑ是石英玻璃分子动力学模拟中常用的对势ꎮ其中ꎬＢＫＳ势从低温到高温的较大温度范围均具有较好性能ꎬ而ＭＳ￣Ｑ势在高温下的传输性能要优于ＢＫＳ势ꎮ目前对于石英玻璃高温表面特性的模拟的报道较为少见ꎬ因此仍需要对两种势函数性能的优劣进行比较ꎮＨｏｓｓｅｉｎｉ等[１３]研究了不同形貌的疏水表面上的水滴行为ꎬ发现表面形貌㊁柱高度㊁空隙率和沉积角是影响表面疏水性的主要参数ꎮ因此ꎬ本模拟在ＳｉＣ表面构建了纳米方柱阵列ꎬ并分别使用粗糙度评定参数Ｒａ和Ｒｍｒ表示柱高和空隙率ꎬ研究了纳米级表面粗糙结构对面向高温模压的ＳｉＯ２/ＳｉＣ接触角的影响ꎮ本文使用分子动力学方法ꎬ对比了分别采用ＢＫＳ和ＭＳ￣Ｑ势函数计算的ＳｉＯ２熔体高温表面张力ꎮ模拟了不同模压温度和ＳｉＣ模具表面粗糙结构ＳｉＯ２的高温接触角ꎬ研究了微尺度下的高温界面润湿特征和ＳｉＯ２熔体表面层结构特点ꎬ本研究的结构框图如图１所示ꎮ１㊀模拟方法和模型１.１㊀界面张力的计算为了保证成形精度并减小工件与模具的粘附ꎬ光学玻璃模压温度一般选择在转变点温度到软化点温度第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９２５㊀之间[１４]ꎬ故选择石英玻璃实际软化点附近温度１９００Ｋ作为模拟温度ꎮＳｉＯ２表面张力模拟系统如图２所示ꎮＳｉＯ２熔体表面模型的建立方法是在模拟盒子中随机加入８００个Ｓｉ原子和１６００个Ｏ原子ꎬ使其密度约为２.２ｇ/ｃｍ３ꎮ随后在５０００Ｋ下驰豫２００ｐｓ以消除初始构型的影响ꎬ并以１００Ｋ/２０ｐｓ的速度降温至１９００Ｋꎮ最后在ｘ轴方向模型两边各添加一厚度为２０Å的真空层ꎬ以避免周期性边界的影响ꎮ面向高温模压的ＳｉＯ２/ＳｉＣ界面张力模拟系统如图３(ｃ)所示ꎬ其中ＳｉＣ是由金刚石结构的β￣ＳｉＣ原胞排列而成ꎮ分别将ＳｉＯ２和ＳｉＣ单独在１９００Ｋ下驰豫２００ｐｓꎬ然后将它们添加进同一模拟盒子中ꎬ使ＳｉＯ２与ＳｉＣ{１００}表面接触并对整个体系进行驰豫ꎮ图１㊀基于分子动力学的ＳｉＣ模具高温模压石英玻璃的物相接触角模拟研究关系框图Ｆｉｇ.１㊀ＳｔｕｄｙｒｅｌａｔｉｏｎｓｈｉｐｏｆｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｃｏｎｔａｃｔａｎｇｌｅｏｆｍｏｌｄｅｄｏｐｔｉｃａｌｇｌａｓｓｉｎＳｉＣｄｉｅｂａｓｅｄｏｎｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓ图２㊀ＳｉＯ２熔体表面张力模拟计算体系Ｆｉｇ.２㊀ＳｉｍｕｌａｔｉｏｎｓｙｓｔｅｍｏｆＳｉＯ２ｍｅｌｔｓｕｒｆａｃｅｔｅｎｓｉｏｎ图３㊀ＳｉＯ２/ＳｉＣ高温界面张力模拟计算体系Ｆｉｇ.３㊀ＳｉｍｕｌａｔｉｏｎｓｙｓｔｅｍｏｆＳｉＯ２/ＳｉＣｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｉｎｔｅｒｆａｃｅｔｅｎｓｉｏｎ㊀㊀表面与界面张力的计算均采用压力张量法[１５]ꎬγ＝１２ʏＬｘ２－Ｌｘ２[ＰＮ(ｘ)－ＰＴ(ｘ)]ｄｘ(１)式中ꎬγ表示表面张力ꎻＬｘ为模拟体系在ｘ方向的长度ꎻ因子１/２是由于模拟系统有两个界面ꎻＰＮ与ＰＴ分别为系统界面法向和切向压力张量分量ꎮＰＮ(ｘ)＝ρｋＢＴ－１Ｖðｉ<ｊｘ２ｉｊｒｉｊｄＵ(ｒｉｊ)ｄｒｉｊ[](２)９２６㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷ＰＴ(ｘ)＝ρｋＢＴ－１Ｖðｉ<ｊｙ２ｉｊ－ｚ２ｉｊ２ｒｉｊｄＵ(ｒｉｊ)ｄｒｉｊ[](３)式中ꎬｋＢ为玻尔兹曼常数ꎻＶ为模拟盒子体积ꎻＵ表示势函数ꎻρ为液体密度ꎻＴ为温度ꎻｒ为原子间距ꎬｘ㊁ｙ㊁ｚ为其在三个坐标轴方向的分量ꎮ为了比较不同势函数模拟界面张力的可靠性ꎬ在模拟中ꎬＳｉＯ２原子间的相互作用先后用ＢＫＳ[１０]和ＭＳ￣Ｑ[１１]势函数进行表征如下ꎬ势函数参数见表１ꎮＵＢｉｊ＝ｑｉｑｊｅ２ｒｉｊ＋Ａｉｊｅ－ｒｉｊρ－Ｃｉｊｒ６ｉｊ(４)ＵＭｉｊ＝ｑｉｑｊｅ２ｒｉｊ＋Ｄ０[ｅ－２α(ｒｉｊ－ｒ０)－２ｅ－α(ｒｉｊ－ｒ０)](５)式中ꎬｑ表示电荷量ꎻＡ㊁Ｃ㊁Ｄ均为与相互作用强度有关的参数ꎻα为与平衡距离有关的参数ꎮ表１㊀ＳｉＯ２势函数参数Ｔａｂｌｅ１㊀ＰａｒａｍｅｔｅｒｓｏｆｐｏｔｅｎｔｉａｌｓｆｏｒＳｉＯ２ＰａｉｒＢＫＳＡｉｊ/(ｋｃａｌ ｍｏｌ－１)ρ/ÅＣｉｊ/(ｋｃａｌ ｍｏｌ－１ Å６)ＭＳ￣ＱＤ０/(ｋｃａｌ ｍｏｌ－１)α/Å－１ｒ０/ÅＳｉ￣Ｓｉ￣￣￣０.１７７３２.０４４３.７６０Ｏ￣Ｏ３２０２５.７９９８０.３６２３４０３５.５８７５０.５３６３１.３７３３.７９１Ｓｉ￣Ｏ４１５１７５.６４２９０.２０５２３０７９.４５５４４５.９９７０２.６５２１.６２８㊀㊀ＳｉＣ原子间的相互作用使用ｔｅｒｓｏｆｆ势函数表征ꎮ由于ＳｉＯ２/ＳｉＣ界面间为范德华作用ꎬ因此使用ＬＪ势函数表征ꎬ其参数来自于ＵＦＦ力场[１６]ꎬ如表２所示ꎮ表２㊀ＳｉＯ２/ＳｉＣ势函数参数Ｔａｂｌｅ２㊀ＰａｒａｍｅｔｅｒｓｏｆｐｏｔｅｎｔｉａｌｆｏｒＳｉＯ２/ＳｉＣｉｎｔｅｒａｃｔｉｏｎＰａｉｒε/(ｋｃａｌ ｍｏｌ－１)σ/ÅＳｉ￣Ｓｉ０.４０２３.８２６Ｏ￣Ｓｉ０.１５５３３.４７２３Ｓｉ￣Ｃ０.２０５３.６２８Ｏ￣Ｃ０.０７９３３.２７４５㊀㊀模拟过程中ꎬ模拟盒子三个方向均为周期性边界ꎮＢＫＳ势函数Ｓｉ和Ｏ原子所带电荷分别为＋２.４ｅ㊁－１.２ｅꎬ截断半径为５.５ÅꎻＭＳ￣Ｑ势函数Ｓｉ和Ｏ原子所带电荷分别为＋１.３１８ｅ㊁－０.６５９ｅꎬ截断半径为９Åꎮ静电力的计算采用ＰＰＰＭ(Ｐａｒｔｉｃｌｅ￣ＰａｒｔｉｃｌｅＰａｒｔｉｃｌｅ￣Ｍｅｓｈ)方法且精度为１０－５ꎬ截断半径为１０Åꎮ时间步长为１ｆｓꎮ模拟均在正则系综(ＮＶＴ)下进行ꎬ调温使用Ｎｏｓé￣Ｈｏｏｖｅｒ方法ꎮ１.２㊀接触角的模拟图４㊀面向高温模压的ＳｉＯ２/ＳｉＣ接触角模拟系统Ｆｉｇ.４㊀ＳｉｍｕｌａｔｉｏｎｓｙｓｔｅｍｏｆＳｉＯ２/ＳｉＣｃｏｎｔａｃｔａｎｇｌｅｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ在模拟系统中ＳｉＯ２与ＳｉＣ分别为液相和固相ꎬ其中ꎬＳｉＣ基底为ＳｉＣ单胞排列而成的超晶胞ꎬ对其作出三点假设:(１)假设ＳｉＣ表面为不存在缺陷的理想表面ꎻ(２)由于ＳｉＣ不同晶面表面能相差较小ꎬ假设其{１００}晶面为与ＳｉＯ２相接触的表面ꎻ(３)假设ＳｉＯ２熔体液滴为理想的球形ꎮ接触角的模拟系统如图４所示ꎬ首先在半径２.３４８Å的球体区域随机添加１２００个Ｓｉ原子和２４００个Ｏ原子ꎬ随后在５０００Ｋ下驰豫２００ｐｓ并以１００Ｋ/２０ｐｓ的速度降温至１９００Ｋ得到模拟所用的ＳｉＯ２液滴模型ꎮ最后将ＳｉＯ２和ＳｉＣ添加进同一模拟盒子ꎬｘ㊁ｙ方向使用周期性边界ꎬｚ方向使用固定和镜像边界ꎬ其他第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９２７㊀设置与上一节相同ꎮ图５㊀微结构阵列模具表面多尺度形貌示意图Ｆｉｇ.５㊀Ｍｕｌｔｉｓｃａｌｅｓｕｒｆａｃｅｔｏｐｏｇｒａｐｈｙｏｆｍｉｃｒｏ￣ｓｔｒｕｃｔｕｒｅａｒｒａｙｄｉｅ为了研究纳米级表面粗糙结构对接触角的影响ꎬ采用文献[１３]的方法构建了方柱形阵列ＳｉＣ壁面ꎮ其结构如图５所示ꎮＷｅｎｚｅｌ[２]和Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ[３]模型可用式(６)㊁(７)表示:ｃｏｓθＷ＝ｒｃｏｓθＹ(６)ｃｏｓθＣ＝ｆｃｏｓθＹ＋ｆ－１(７)式中ꎬθＹ为本征接触角ꎻθＷ与θＣ为Ｗｅｎｚｅｌ和Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ模型接触角ꎻｒ为粗糙度因子ꎬ表示表面的实际面积与表观面积之比ꎻｆ为接触面积分数ꎬ表示液滴与基底实际接触的面积和基底表观面积之比ꎮｒ与ｆ可用粗糙度评定参数表示为ｒ＝１＋８ＲｍｒＲａＳｍ(８)ｆ＝Ｒｍｒ２(９)式中ꎬＲａ为轮廓算数平均偏差ꎬＲｍｒ为轮廓支承长度率ꎬＳｍ为轮廓微观不平度的平均间距ꎮ模拟中若Ｒｍｒ㊁Ｓｍ取恒值ꎬ则Ｒａ取０.２５㊁０.５㊁１㊁１.５倍晶格参数ꎬ分别对应ｒ值为１.５㊁２.０㊁３.０㊁４.０ꎻ若Ｒａ取恒值ꎬ改变柱间距为２㊁１.５㊁１㊁０.５倍晶格参数使Ｒｍｒ取０.３３㊁０.４０㊁０.５０㊁０.６７ꎬ分别对应ｆ值为０.１１㊁０.１６㊁０.２５㊁０.４４ꎮ２㊀模拟结果与讨论２.１㊀界面张力图６为１９００Ｋ下ＢＫＳ与ＭＳ￣Ｑ模型中ＳｉＯ２熔体表面张力随时间的演化ꎮ由于使用压力张量法时系统需要较长时间稳定ꎬ模拟连续进行了２０ｎｓꎮ为了降低模拟初期结构不稳定的影响ꎬ在计算累计平均值时采用最后１０ｎｓ数据ꎮＢＫＳ与ＭＳ￣Ｑ势函数的计算结果分别为０.４１０Ｎ/ｍ和０.３９０Ｎ/ｍꎮ由文献[１７]可知ꎬＳｉＯ２在１９００Ｋ下的实际表面张力应为０.３０２Ｎ/ｍ左右ꎮ考虑到模拟本身的准确性以及实际实验过程中ꎬ氧分子等充当表面活性剂降低了测定的表面张力值ꎬ因此模拟值高出实验值０.０７~０.１１Ｎ/ｍ左右是可接受的[１８]ꎮ可以看出ＭＳ￣Ｑ势函数的模拟结果比ＢＫＳ势函数更为接近实验值ꎮ图６㊀ＢＫＳ势函数与ＭＳ￣Ｑ势函数下ＳｉＯ２熔体表面张力演化Ｆｉｇ.６㊀ＥｖｏｌｕｔｉｏｎｏｆＳｉＯ２ｍｅｌｔｓｕｒｆａｃｅｔｅｎｓｉｏｎｕｓｉｎｇＢＫＳｐｏｔｅｎｔｉａｌａｎｄＭＳ￣Ｑｐｏｔｅｎｔｉａｌ图７为两种不同势函数模拟的面向模压的ＳｉＯ２/ＳｉＣ高温界面张力演化图ꎮ界面张力系统达到稳定时间相对较短ꎬ因此模拟时间设为１５ｎｓꎬ取最后４.５ｎｓ时间段的数据计算平均值ꎮ从图中可见ꎬＢＫＳ和ＭＳ￣Ｑ势函数的模拟结果分别为０.８４６Ｎ/ｍ和０.６８２Ｎ/ｍꎮ同时ꎬ计算了两种表面模型的一维密度分布和径向分布函数(ＲａｄｉａｌＤｉｓｔｒｉｂｕｔｉｏｎＦｕｎｃｔｉｏｎꎬＲＤＦ)ꎮ图８(ａ)为模型中心至模拟盒子边缘的密度分布ꎬＢＫＳ与ＭＳ￣Ｑ模型的表面均存在类似汽液共存界面的低密度相ꎬ其厚度分别为４Å与６ÅꎮＢＫＳ模型内部的密度约为９２８㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷２.６９ｇ/ｃｍ３ꎬ大于ＭＳ￣Ｑ的２.２５ｇ/ｃｍ３ꎮ这说明ＢＫＳ表面模型的密度大于ＭＳ￣Ｑ模型ꎬ因此原子间距要小于ＭＳ￣Ｑ模型ꎬ具有更强的原子间作用力ꎬ从而具有更大的表面张力ꎮ图８(ｂ)中ＲＤＦ的计算结果也支持了这一解释ꎮＲＤＦ图中前三个峰代表Ｏ￣Ｏ㊁Ｓｉ￣Ｏ㊁Ｓｉ￣Ｓｉ原子对的第一近邻ꎬ峰值点的横坐标即原子对的键长ꎬ其数值如图８(ｂ)所示ꎬＭＳ￣Ｑ模型中数量较多且强度较大的Ｓｉ￣Ｏ㊁Ｏ￣Ｏ键长均大于ＢＫＳ模型ꎬ从而导致其表面张力小于ＢＫＳ模型ꎮ综上所述ꎬ高温模压时ＳｉＯ２表面性质的模拟可优先选择ＭＳ￣Ｑ势函数ꎮ图７㊀面向高温模压的ＳｉＯ２/ＳｉＣ界面张力演化过程Ｆｉｇ.７㊀ＥｖｏｌｕｔｉｏｎｏｆＳｉＯ２/ＳｉＣｉｎｔｅｒｆａｃｅｔｅｎｓｉｏｎｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ图８㊀使用ＢＫＳ与ＭＳ￣Ｑ势函数的ＳｉＯ２表面结构对比Ｆｉｇ.８㊀ＣｏｍｐａｒｉｓｏｎｏｆＳｉＯ２ｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｕｓｉｎｇＢＫＳａｎｄＭＳ￣Ｑｐｏｔｅｎｔｉａｌ２.２㊀表面粗糙结构对接触角的影响液滴在理想光滑表面上的接触角称为本征接触角ꎮ由于微纳尺度下液滴表面存在一定的密度和压力涨落ꎬ为了获得具有统计学意义的接触角值ꎬ在处理数据时采用等密度拟合曲线法[１９]ꎮ最后得到面向高温模压的ＳｉＯ２/ＳｉＣ的本征接触角为１１９.２５ʎꎮ根据杨氏方程[１]ꎬ可以利用上节模拟得到的界面张力对本征接触角进行检验ꎮｃｏｓθ＝γｓｖ－γｓｌγｌｖ(１０)式中ꎬγｓｖ表示固￣气界面能ꎻγｓｌ与γｌｖ分别表示固￣液界面能和气￣液界面能ꎬ它们在数值上与固￣液界面张力和液体表面张力相等ꎮ由上节可知ꎬγｓｌ和γｌｖ的值分别为０.６８２Ｎ/ｍ和０.３９０Ｎ/ｍꎮ由于本模拟中采用的ＳｉＯ２/ＳｉＣ相互作用参数来自ＵＦＦ力场ꎬ因此选择文献[２０]中使用ＵＦＦ力场算得的ＳｉＣ{１００}晶面表面能数值０.５０２Ｎ/ｍ代入式(１０)ꎮ最终得到的本征接触角理论值为１１７.４９ʎꎬ与模拟数值差别不大ꎮ图９是ＳｉＯ２熔体在具有不同粗糙度因子的ＳｉＣ模具表面上接触角的模拟ꎬ从图中可见ꎬ接触角总体值在１３５ʎ左右波动ꎮ接触角变化不大则表明固液相互作用能较为稳定ꎬ但图９(ｂ)显示ꎬ在ｒ＝１.５时固液相互作用能较大ꎮ由图１０可以看出ꎬ当粗糙度因子ｒ＝１.５时ꎬＳｉＯ２熔体浸润了沟槽ꎬ出现 钉扎 现象ꎬ此时润湿第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９２９㊀接近于Ｗｅｎｚｅｌ状态ꎮ 钉扎 现象增强了固液界面的摩擦ꎬ液滴的铺展需要克服更大的能垒ꎬ因此虽然固液相互作用能较大但并未使液滴接触角减小ꎻｒ>１.５时ꎬＳｉＯ２均处于Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿状态ꎬ因此接触角不再受粗糙度因子变化的影响ꎮ从Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ状态到Ｗｅｎｚｅｌ状态的过渡称为润湿转变ꎮ润湿转变可以通过施加压力实现ꎬ但临界转变压力会随着纳米柱的高度减小而减小[２１]ꎻ当其减小到一定程度ꎬ就可能仅仅依靠系统压力实现润湿转变ꎮ本研究中ꎬ当ｒ>１.５时系统压力无法使液滴保持Ｗｅｎｚｅｌ状态ꎬ因此转变为Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ状态ꎻ而此状态下ꎬ界面的固液相互作用比Ｗｅｎｚｅｌ状态更低ꎬ液滴在固体表面的扩散系数变小[２２]ꎬ说明此时ＳｉＯ２熔体在剪切作用下更易在ＳｉＣ模具表面滑移ꎬ从而导致粘性摩擦作用较小ꎮ图９㊀ＳｉＣ模具表面粗糙度对ＳｉＯ２熔体液滴接触角和两者相互作用能的影响Ｆｉｇ.９㊀ＥｆｆｅｃｔｏｆＳｉＣｄｉｅｒｏｕｇｈｎｅｓｓｏｎｃｏｎｔａｃｔａｎｇｌｅｏｆＳｉＯ２ｍｅｌｔｄｒｏｐｌｅｔａｎｄｉｎｔｅｒａｃｔｉｏｎｅｎｅｒｇｙ图１０㊀ＳｉＣ模具表面粗糙度对ＳｉＯ２熔体液滴接触状态的影响Ｆｉｇ.１０㊀ＥｆｆｅｃｔｏｆＳｉＣｄｉｅｒｏｕｇｈｎｅｓｓｏｎｃｏｎｔａｃｔｓｔａｔｅｏｆＳｉＯ２ｍｅｌｔｄｒｏｐｌｅｔ图１１㊀接触面积分数对面向高温模压的ＳｉＯ２/ＳｉＣ接触状态的影响Ｆｉｇ.１１㊀ＥｆｆｅｃｔｏｆｃｏｎｔａｃｔａｒｅａｆｒａｃｔｉｏｎｏｎＳｉＯ２/ＳｉＣｃｏｎｔａｃｔｓｔａｔｅｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ图１２㊀接触面积分数对面向高温模压的ＳｉＯ２/ＳｉＣ接触角和相互作用能的影响Ｆｉｇ.１２㊀ＥｆｆｅｃｔｏｆｃｏｎｔａｃｔａｒｅａｆｒａｃｔｉｏｎｏｎＳｉＯ２/ＳｉＣｃｏｎｔａｃｔａｎｇｌｅａｎｄｉｎｔｅｒａｃｔｉｏｎｅｎｅｒｇｙｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ９３０㊀陶㊀瓷硅酸盐通报㊀㊀㊀㊀㊀㊀第３９卷图１１为粗糙度因子ｒ>１.５时ＳｉＣ模具表面接触面积分数对ＳｉＯ２熔体液滴接触状态的影响ꎬ此时ＳｉＯ２始终呈现Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿状态ꎬ液滴随接触面积分数的增大逐渐在ＳｉＣ表面铺展ꎮ图１２(ａ)对比了模拟接触角和理论接触角ꎬ虽然两者存在一定误差ꎬ但其都随着接触面积分数的增大而减小ꎻ同时图１２(ｂ)显示接触角越小ꎬＳｉＯ２/ＳｉＣ高温模压界面具有越强的相互作用ꎮ表面粘着力和热应力是脱模力的两个组成部分[２３]ꎬ减小轮廓支承长度率即减小了ＳｉＯ２与ＳｉＣ模具界面的实际接触面积ꎬ从而能减小工件和模具之间的表面粘着力ꎮ而Ｃａｓｓｉｅ￣Ｂａｘｔｅｒ润湿模式无 钉扎 现象ꎬ减小了工件和模具的传热面积ꎬ在一定程度上减小了热应力ꎮ因此适当减小Ｒｍｒ值可以降低以及工件与模具之间的脱模力ꎬ从而减小模具磨损ꎬ提高寿命ꎮ同时也缩减了模具制造过程中抛光工序的工作量ꎮ２.３㊀温度对接触角的影响图１３㊀面向高温模压的ＳｉＯ２/ＳｉＣ接触角与温度的关系Ｆｉｇ.１３㊀ＲｅｌａｔｉｏｎｓｈｉｐｂｅｔｗｅｅｎＳｉＯ２/ＳｉＣｃｏｎｔａｃｔａｎｇｌｅａｎｄｔｅｍｐｅｒａｔｕｒｅｆｏｒｈｉｇｈｔｅｍｐｅｒａｔｕｒｅｍｏｌｄｉｎｇ图１３为接触面积分数为０.２５时温度对面向高温模压的ＳｉＯ２/ＳｉＣ接触角的影响ꎮ从图可见ꎬ接触角值随温度的增大而减小ꎬ说明温度升高ＳｉＯ２的表面张力减小ꎻ当温度高于２３００Ｋ时ꎬ接触角的变化较大ꎬ这是因为使用ＭＳ￣Ｑ势函数模型的ＳｉＯ２自扩散激活温度约为２３００Ｋ左右[１２]ꎮ自扩散激活温度可以视作石英玻璃的软化点温度ꎬ在此温度附近ꎬ石英玻璃的结构发生较大变化ꎬ松动和新生的化学键均大大增加[２４]ꎮ图１４(ａ)是不同温度下ＳｉＯ２熔体表面的结构特点ꎬ由图可见ꎬＳｉＯ２熔体表面存在一个等密度层ꎻ在等密度层外部ꎬ密度随温度升高而增大ꎻ在等密度层内部ꎬ密度随温度升高而减小ꎻ这意味着ＳｉＯ２熔体内部的原子在高温作用下逐渐向外扩散ꎬ这种密度梯度的减小使ＳｉＯ２表面结构更加松散ꎬ减小了其表面张力ꎮ图１４(ｂ)还表明ꎬ温度对ＳｉＯ２熔体表面原子化学键的键长并未产生较大影响ꎬ但随着温度升高ꎬ各化学键对应的峰值依次减小ꎬ表明表面层原子的无序程度增加ꎮ同时ꎬ对比图８(ｂ)可以发现ꎬ表面层中Ｓｉ￣Ｓｉ对应的峰值变得较为微弱甚至消失ꎬ这说明在表面层硅原子数量较少ꎬ而氧原子大量聚集ꎬＯ￣Ｏ键的强度远小于Ｓｉ￣Ｏ键强度ꎬ这可能是ＳｉＯ２表面张力随温度减小的另一个原因ꎮ图１４㊀不同温度下ＳｉＯ２熔体表面结构特点Ｆｉｇ.１４㊀ＳｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆＳｉＯ２ｍｅｌｔｉｎｄｉｆｆｅｒｅｎｔｔｅｍｐｅｒａｔｕｒｅ３㊀结㊀论采用分子动力学方法模拟了ＳｉＯ２熔体的界面结构ꎬ将ＳｉＣ模具纳米级表面理想化为纳米方柱阵列ꎬ研㊀第３期吴㊀悠:ＳｉＣ模具高温模压石英玻璃物相接触角的分子动力学模拟９３１究了粗糙度和温度对面向模压的ＳｉＯ２/ＳｉＣ高温接触角的影响ꎬ得到以下结论:(１)使用ＭＳ￣Ｑ势函数模拟的ＳｉＯ２熔体表面张力比ＢＫＳ势函数计算的结果与实验值更为接近ꎬ因此模拟ＳｉＯ２高温熔体的表面性质使用ＭＳ￣Ｑ势函数更为合理ꎮ(２)在１９００Ｋ的模压温度下ꎬ当粗糙度因子ｒ>１.５时ꎬＲａ的变化对接触角值无明显影响ꎮＲｍｒ值减小使得接触面积分数ｆ减小ꎬ接触角值随之增大ꎮ此时润湿模式从Ｗｅｎｚｅｌ转变为Ｃａｓｓｉｅ￣Ｂａｘｔｅｒꎬ减小了工件￣模具之间的摩擦ꎮ由于热应力和界面粘着力的减小ꎬ石英玻璃光学元件模压后的脱模力将会降低ꎮ同时ꎬ由于降低了Ｒｍｒ值的要求ꎬＳｉＣ模具加工过程中抛光工序的工作量也相应得到了减小ꎮ(３)面向高温模压的ＳｉＯ２/ＳｉＣ的接触角随模压温度升高而减小ꎮ当模压温度超过２３００Ｋ时ꎬ接触角变化率显著增大ꎮ因此模压温度选择在２３００Ｋ以下可以降低模压时模具因玻璃熔体的粘附造成的磨损ꎮ参考文献[１]㊀ＹｏｕｎｇＴＨ.Ａｎｅｓｓａｙｏｎｔｈｅｃｏｈｅｓｉｏｎｏｆｌｉｑｕｉｄｓ[Ｊ].Ｐｈｉｌ.Ｔｒａｎｓ.Ｒｏｙ.Ｓｏｃ.Ｌｏｎｄｏｎꎬ１８０５ꎬ９５:６５￣８７.[２]㊀ＷｅｎｚｅｌＲＮ.Ｒｅｓｉｓｔａｎｃｅｏｆｓｏｌｉｄｓｕｒｆａｃｅｓｔｏｗｅｔｔｉｎｇｂｙｗａｔｅｒ[Ｊ].Ｉｎｄ.Ｅｎｇ.Ｃｈｅｍ.ꎬ１９３６ꎬ２８:９８８￣９４.[３]㊀ＣａｓｓｉｅＡＢＤꎬＢａｘｔｅｒＳ.Ｗｅｔｔａｂｉｌｉｔｙｏｆｐｏｒｏｕｓｓｕｒｆａｃｅｓ[Ｊ].Ｔｒａｎｓ.Ｆａｒａｄａｙ.Ｓｏｃ.ꎬ１９４４ꎬ４０:５４６￣５１.[４]㊀陈㊀光.新材料概论[Ｍ].北京:科学出版社ꎬ２００３:４６￣４８.[５]㊀和丽芳.纳米碳化硅材料的制备及应用[Ｄ].太原:山西大学ꎬ２０１１.[６]㊀徐㊀威ꎬ兰㊀忠ꎬ彭本利ꎬ等.微液滴在不同能量表面上润湿状态的分子动力学模拟[Ｊ].物理学报ꎬ２０１６(２１):２１６８０１.[７]㊀王㊀龙.金属液滴在碳纳米材料表面的润湿与融合[Ｄ].济南:山东大学ꎬ２０１７.[８]㊀ＤｉｎｇＹꎬＺｈａｎｇＹꎬＺｈａｎｇＦꎬｅｔａｌ.ＭｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｔｕｄｙｏｆｔｈｅｓｔｒｕｃｔｕｒｅｉｎｖｉｔｒｅｏｕｓｓｉｌｉｃａｗｉｔｈＣＯＭＰＡＳＳｆｏｒｃｅｆｉｅｌｄａｔｅｌｅｖａｔｅｄｔｅｍｐｅｒａｔｕｒｅｓ[Ｃ].ＭａｔｅｒｉａｌｓＳｃｉｅｎｃｅＦｏｒｕｍ.２００７.[９]㊀ＢｅｅｓｔＢＷＨＶꎬＫｒａｍｅｒＧＪꎬＳａｎｔｅｎＲＡＶ.Ｆｏｒｃｅｆｉｅｌｄｓｆｏｒｓｉｌｉｃａｓａｎｄａｌｕｍｉｎｏｐｈｏｓｐｈａｔｅｓｂａｓｅｄｏｎａｂｉｎｉｔｉｏｃａｌｃｕｌａｔｉｏｎｓ[Ｊ].ＰｈｙｓｉｃａｌＲｅｖｉｅｗＬｅｔｔｅｒｓꎬ１９９０ꎬ６４(１６):１９５５.[１０]㊀ＳｕｎｄａｒａｒａｍａｎＳꎬＣｈｉｎｇＷＹꎬＨｕａｎｇＬ.Ｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆｓｉｌｉｃａｇｌａｓｓｐｒｅｄｉｃｔｅｄｂｙａｐａｉｒｗｉｓｅｐｏｔｅｎｔｉａｌｉｎｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｉｍｕｌａｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆＮｏｎ￣ＣｒｙｓｔａｌｌｉｎｅＳｏｌｉｄｓꎬ２０１６ꎬ４４５￣４４６:１０２￣１０９.[１１]㊀ＤｅｍｉｒａｌｐＥꎬＣａｇｉｎꎬＴａｈｉｒꎬＧｏｄｄａｒｄＷＡ.Ｍｏｒｓｅｓｔｒｅｔｃｈｐｏｔｅｎｔｉａｌｃｈａｒｇｅｅｑｕｉｌｉｂｒｉｕｍｆｏｒｃｅｆｉｅｌｄｆｏｒｃｅｒａｍｉｃｓ:ａｐｐｌｉｃａｔｉｏｎｔｏｔｈｅｑｕａｒｔｚ￣ｓｔｉｓｈｏｖｉｔｅｐｈａｓｅｔｒａｎｓｉｔｉｏｎａｎｄｔｏｓｉｌｉｃａｇｌａｓｓ[Ｊ].ＰｈｙｓｉｃａｌＲｅｖｉｅｗＬｅｔｔｅｒｓꎬ１９９９ꎬ８２(８):１７０８￣１７１１.[１２]㊀丁元法ꎬ张㊀跃ꎬ张凡伟ꎬ等.石英玻璃高温分子动力学模拟中的势函数[Ｊ].物理化学学报ꎬ２００８ꎬ２４(５):７８８￣７９２.[１３]㊀ＨｏｓｓｅｉｎｉＳꎬＳａｖａｌｏｎｉＨꎬＳｈａｈｒａｋｉＭＧ.Ｉｎｆｌｕｅｎｃｅｏｆｓｕｒｆａｃｅｍｏｒｐｈｏｌｏｇｙａｎｄｎａｎｏ￣ｓｔｒｕｃｔｕｒｅｏｎｈｙｄｒｏｐｈｏｂｉｃｉｔｙ:Ａｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓａｐｐｒｏａｃｈ[Ｊ].ＡｐｐｌｉｅｄＳｕｒｆａｃｅＳｃｉｅｎｃｅꎬ２０１９ꎬ４８５:５３６￣５４６.[１４]㊀周书强.小口径非球面玻璃透镜的模压仿真及实验研究[Ｄ].长沙:湖南大学ꎬ２０１６.[１５]㊀ＨｅｒｄｅｓＣꎬＥｒｖｉｋＡꎬＭｅｊíａꎬｅｔａｌ.Ｐｒｅｄｉｃｔｉｏｎｏｆｔｈｅｗａｔｅｒ/ｏｉｌｉｎｔｅｒｆａｃｉａｌｔｅｎｓｉｏｎｆｒｏｍｍｏｌｅｃｕｌａｒｓｉｍｕｌａｔｉｏｎｓｕｓｉｎｇｔｈｅｃｏａｒｓｅ￣ｇｒａｉｎｅｄＳＡＦＴ￣γＭｉｅｆｏｒｃｅｆｉｅｌｄ[Ｊ].ＦｌｕｉｄＰｈａｓｅＥｑｕｉｌｉｂｒｉａꎬ２０１８ꎬ４７６:９￣１５.[１６]㊀ＲａｐｐｅＡＫꎬＣａｓｅｗｉｔＣＪꎬＣｏｌｗｅｌｌＫＳꎬｅｔａｌ.ＵＦＦꎬａｆｕｌｌｐｅｒｉｏｄｉｃｔａｂｌｅｆｏｒｃｅｆｉｅｌｄｆｏｒｍｏｌｅｃｕｌａｒｍｅｃｈａｎｉｃｓａｎｄｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｉｍｕｌａｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆｔｈｅＡｍｅｒｉｃａｎＣｈｅｍｉｃａｌＳｏｃｉｅｔｙꎬ１９９２ꎬ１１４(２５):１００２４￣１００３５.[１７]㊀ＷｕＣꎬＣｈｅｎＧ.ＣｏｍｐｕｔａｔｉｏｎａｌｍｏｄｅｌｏｎｓｕｒｆａｃｅｔｅｎｓｉｏｎｏｆＣａＯ￣ＭｎＯ￣ＳｉＯ２ｓｌａｇｓｙｓｔｅｍ[Ｊ].ＨｏｔＷｏｒｋｉｎｇＴｅｃｈｎｏｌｏｇｙꎬ２０１４ꎬ４３(１９):５８￣６２. [１８]㊀ＢｅｌａｓｈｃｈｅｎｋｏＤＫꎬＯｓｔｒｏｖｓｋｉｉＯＩ.Ｃｏｍｐｕｔｅｒｓｉｍｕｌａｔｉｏｎｏｆｓｍａｌｌｎｏｎｃｒｙｓｔａｌｌｉｎｅｓｉｌｉｃａｃｌｕｓｔｅｒｓ[Ｊ].ＩｎｏｒｇａｎｉｃＭａｔｅｒｉａｌｓꎬ２００２ꎬ３８(９):９１７￣９２１. [１９]㊀王宝和ꎬ李㊀群.接触角的研究现状及其在凝胶干燥中的作用[Ｊ].干燥技术与设备ꎬ２０１４ꎬ１２(１):３９￣４６.[２０]㊀罗晓光.ＺｒＢ２/ＳｉＣ复合材料界面结构的分子动力学模拟[Ｄ].哈尔滨:哈尔滨工业大学ꎬ２００７.[２１]㊀ＬｉｕＢꎬＬａｎｇｅＦＦ.Ｐｒｅｓｓｕｒｅｉｎｄｕｃｅｄｔｒａｎｓｉｔｉｏｎｂｅｔｗｅｅｎｓｕｐｅｒｈｙｄｒｏｐｈｏｂｉｃｓｔａｔｅｓ:ｃｏｎｆｉｇｕｒａｔｉｏｎｄｉａｇｒａｍｓａｎｄｅｆｆｅｃｔｏｆｓｕｒｆａｃｅｆｅａｔｕｒｅｓｉｚｅ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｌｌｏｉｄ＆ＩｎｔｅｒｆａｃｅＳｃｉｅｎｃｅꎬ２００６ꎬ２９８(２):８９９￣９０９.[２２]㊀ＫａｌｉｎＭꎬＰｏｌａｊｎａｒＭ.Ｔｈｅｅｆｆｅｃｔｏｆｗｅｔｔｉｎｇａｎｄｓｕｒｆａｃｅｅｎｅｒｇｙｏｎｔｈｅｆｒｉｃｔｉｏｎａｎｄｓｌｉｐｉｎｏｉｌ￣ｌｕｂｒｉｃａｔｅｄｃｏｎｔａｃｔｓ[Ｊ].ＴｒｉｂｏｌｏｇｙＬｅｔｔｅｒｓꎬ２０１３ꎬ５２(２):１８５￣１９４.[２３]㊀叶㊀伟.强流脉冲电子束处理对热模压脱模性能的影响[Ｄ].重庆:重庆理工大学ꎬ２０１３.[２４]㊀ＷａｎｇＣꎬＫｕｚｕｕＮꎬＴａｍａｉＹ.Ｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓｓｔｕｄｙｏｎｓｕｒｆａｃｅｓｔｒｕｃｔｕｒｅｏｆｈｏｔｗｏｒｋｉｎｇｔｅｃｈｎｏｌｏｇｙα￣ＳｉＯ２ｂｙｃｈａｒｇｅｅｑｕｉｌｉｂｒａｔｉｏｎｍｅｔｈｏｄ[Ｊ].ＪｏｕｒｎａｌｏｆＮｏｎ￣ＣｒｙｓｔａｌｌｉｎｅＳｏｌｉｄｓꎬ２００３ꎬ３１８:１３１￣１４１.。
稳定的高功率激光系统在高级引力波探测器中的应用
Stabilized high-power laser system forthe gravitational wave detector advancedLIGOP.Kwee,1,∗C.Bogan,2K.Danzmann,1,2M.Frede,4H.Kim,1P.King,5J.P¨o ld,1O.Puncken,3R.L.Savage,5F.Seifert,5P.Wessels,3L.Winkelmann,3and B.Willke21Max-Planck-Institut f¨u r Gravitationsphysik(Albert-Einstein-Institut),Hannover,Germany2Leibniz Universit¨a t Hannover,Hannover,Germany3Laser Zentrum Hannover e.V.,Hannover,Germany4neoLASE GmbH,Hannover,Germany5LIGO Laboratory,California Institute of Technology,Pasadena,California,USA*patrick.kwee@aei.mpg.deAbstract:An ultra-stable,high-power cw Nd:Y AG laser system,devel-oped for the ground-based gravitational wave detector Advanced LIGO(Laser Interferometer Gravitational-Wave Observatory),was comprehen-sively ser power,frequency,beam pointing and beamquality were simultaneously stabilized using different active and passiveschemes.The output beam,the performance of the stabilization,and thecross-coupling between different stabilization feedback control loops werecharacterized and found to fulfill most design requirements.The employedstabilization schemes and the achieved performance are of relevance tomany high-precision optical experiments.©2012Optical Society of AmericaOCIS codes:(140.3425)Laser stabilization;(120.3180)Interferometry.References and links1.S.Rowan and J.Hough,“Gravitational wave detection by interferometry(ground and space),”Living Rev.Rel-ativity3,1–3(2000).2.P.R.Saulson,Fundamentals of Interferometric Gravitational Wave Detectors(World Scientific,1994).3.G.M.Harry,“Advanced LIGO:the next generation of gravitational wave detectors,”Class.Quantum Grav.27,084006(2010).4. B.Willke,“Stabilized lasers for advanced gravitational wave detectors,”Laser Photon.Rev.4,780–794(2010).5.P.Kwee,“Laser characterization and stabilization for precision interferometry,”Ph.D.thesis,Universit¨a t Han-nover(2010).6.K.Somiya,Y.Chen,S.Kawamura,and N.Mio,“Frequency noise and intensity noise of next-generationgravitational-wave detectors with RF/DC readout schemes,”Phys.Rev.D73,122005(2006).7. B.Willke,P.King,R.Savage,and P.Fritschel,“Pre-stabilized laser design requirements,”internal technicalreport T050036-v4,LIGO Scientific Collaboration(2009).8.L.Winkelmann,O.Puncken,R.Kluzik,C.Veltkamp,P.Kwee,J.Poeld,C.Bogan,B.Willke,M.Frede,J.Neu-mann,P.Wessels,and D.Kracht,“Injection-locked single-frequency laser with an output power of220W,”Appl.Phys.B102,529–538(2011).9.T.J.Kane and R.L.Byer,“Monolithic,unidirectional single-mode Nd:Y AG ring laser,”Opt.Lett.10,65–67(1985).10.I.Freitag,A.T¨u nnermann,and H.Welling,“Power scaling of diode-pumped monolithic Nd:Y AG lasers to outputpowers of several watts,”mun.115,511–515(1995).11.M.Frede,B.Schulz,R.Wilhelm,P.Kwee,F.Seifert,B.Willke,and D.Kracht,“Fundamental mode,single-frequency laser amplifier for gravitational wave detectors,”Opt.Express15,459–465(2007).#161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 1061712. A.D.Farinas,E.K.Gustafson,and R.L.Byer,“Frequency and intensity noise in an injection-locked,solid-statelaser,”J.Opt.Soc.Am.B12,328–334(1995).13.R.Bork,M.Aronsson,D.Barker,J.Batch,J.Heefner,A.Ivanov,R.McCarthy,V.Sandberg,and K.Thorne,“New control and data acquisition system in the Advanced LIGO project,”Proc.of Industrial Control And Large Experimental Physics Control System(ICALEPSC)conference(2011).14.“Experimental physics and industrial control system,”/epics/.15.P.Kwee and B.Willke,“Automatic laser beam characterization of monolithic Nd:Y AG nonplanar ring lasers,”Appl.Opt.47,6022–6032(2008).16.P.Kwee,F.Seifert,B.Willke,and K.Danzmann,“Laser beam quality and pointing measurement with an opticalresonator,”Rev.Sci.Instrum.78,073103(2007).17. A.R¨u diger,R.Schilling,L.Schnupp,W.Winkler,H.Billing,and K.Maischberger,“A mode selector to suppressfluctuations in laser beam geometry,”Opt.Acta28,641–658(1981).18. B.Willke,N.Uehara,E.K.Gustafson,R.L.Byer,P.J.King,S.U.Seel,and R.L.Savage,“Spatial and temporalfiltering of a10-W Nd:Y AG laser with a Fabry-Perot ring-cavity premode cleaner,”Opt.Lett.23,1704–1706 (1998).19.J.H.P¨o ld,“Stabilization of the Advanced LIGO200W laser,”Diploma thesis,Leibniz Universit¨a t Hannover(2009).20. E.D.Black,“An introduction to Pound-Drever-Hall laser frequency stabilization,”Am.J.Phys.69,79–87(2001).21.R.W.P.Drever,J.L.Hall,F.V.Kowalski,J.Hough,G.M.Ford,A.J.Munley,and H.Ward,“Laser phase andfrequency stabilization using an optical resonator,”Appl.Phys.B31,97–105(1983).22. A.Bullington,ntz,M.Fejer,and R.Byer,“Modal frequency degeneracy in thermally loaded optical res-onators,”Appl.Opt.47,2840–2851(2008).23.G.Mueller,“Beam jitter coupling in Advanced LIGO,”Opt.Express13,7118–7132(2005).24.V.Delaubert,N.Treps,ssen,C.C.Harb,C.Fabre,m,and H.-A.Bachor,“TEM10homodynedetection as an optimal small-displacement and tilt-measurement scheme,”Phys.Rev.A74,053823(2006). 25.P.Kwee,B.Willke,and K.Danzmann,“Laser power noise detection at the quantum-noise limit of32A pho-tocurrent,”Opt.Lett.36,3563–3565(2011).26. A.Araya,N.Mio,K.Tsubono,K.Suehiro,S.Telada,M.Ohashi,and M.Fujimoto,“Optical mode cleaner withsuspended mirrors,”Appl.Opt.36,1446–1453(1997).27.P.Kwee,B.Willke,and K.Danzmann,“Shot-noise-limited laser power stabilization with a high-power photodi-ode array,”Opt.Lett.34,2912–2914(2009).28. ntz,P.Fritschel,H.Rong,E.Daw,and G.Gonz´a lez,“Quantum-limited optical phase detection at the10−10rad level,”J.Opt.Soc.Am.A19,91–100(2002).1.IntroductionInterferometric gravitational wave detectors[1,2]perform one of the most precise differential length measurements ever.Their goal is to directly detect the faint signals of gravitational waves emitted by astrophysical sources.The Advanced LIGO(Laser Interferometer Gravitational-Wave Observatory)[3]project is currently installing three second-generation,ground-based detectors at two observatory sites in the USA.The4kilometer-long baseline Michelson inter-ferometers have an anticipated tenfold better sensitivity than theirfirst-generation counterparts (Inital LIGO)and will presumably reach a strain sensitivity between10−24and10−23Hz−1/2.One key technology necessary to reach this extreme sensitivity are ultra-stable high-power laser systems[4,5].A high laser output power is required to reach a high signal-to-quantum-noise ratio,since the effect of quantum noise at high frequencies in the gravitational wave readout is reduced with increasing circulating laser power in the interferometer.In addition to quantum noise,technical laser noise coupling to the gravitational wave channel is a major noise source[6].Thus it is important to reduce the coupling of laser noise,e.g.by optical design or by exploiting symmetries,and to reduce laser noise itself by various active and passive stabilization schemes.In this article,we report on the pre-stabilized laser(PSL)of the Advanced LIGO detector. The PSL is based on a high-power solid-state laser that is comprehensively stabilized.One laser system was set up at the Albert-Einstein-Institute(AEI)in Hannover,Germany,the so called PSL reference system.Another identical PSL has already been installed at one Advanced LIGO site,the one near Livingston,LA,USA,and two more PSLs will be installed at the second #161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10618site at Hanford,WA,USA.We have characterized the reference PSL and thefirst observatory PSL.For this we measured various beam parameters and noise levels of the output beam in the gravitational wave detection frequency band from about10Hz to10kHz,measured the performance of the active and passive stabilization schemes,and determined upper bounds for the cross coupling between different control loops.At the time of writing the PSL reference system has been operated continuously for more than18months,and continues to operate reliably.The reference system delivered a continuous-wave,single-frequency laser beam at1064nm wavelength with a maximum power of150W with99.5%in the TEM00mode.The active and passive stabilization schemes efficiently re-duced the technical laser noise by several orders of magnitude such that most design require-ments[5,7]were fulfilled.In the gravitational wave detection frequency band the relative power noise was as low as2×10−8Hz−1/2,relative beam pointingfluctuations were as low as1×10−7Hz−1/2,and an in-loop measurement of the frequency noise was consistent with the maximum acceptable frequency noise of about0.1HzHz−1/2.The cross couplings between the control loops were,in general,rather small or at the expected levels.Thus we were able to optimize each loop individually and observed no instabilities due to cross couplings.This stabilized laser system is an indispensable part of Advanced LIGO and fulfilled nearly all design goals concerning the maximum acceptable noise levels of the different beam pa-rameters right after installation.Furthermore all or a subset of the implemented stabilization schemes might be of interest for many other high-precision optical experiments that are limited by laser noise.Besides gravitational wave detectors,stabilized laser systems are used e.g.in the field of optical frequency standards,macroscopic quantum objects,precision spectroscopy and optical traps.In the following section the laser system,the stabilization scheme and the characterization methods are described(Section2).Then,the results of the characterization(Section3)and the conclusions(Section4)are presented.ser system and stabilizationThe PSL consists of the laser,developed and fabricated by Laser Zentrum Hannover e.V.(LZH) and neoLASE,and the stabilization,developed and integrated by AEI.The optical components of the PSL are on a commercial optical table,occupying a space of about1.5×3.5m2,in a clean,dust-free environment.At the observatory sites the optical table is located in an acoustically isolated cleanroom.Most of the required electronics,the laser diodes for pumping the laser,and water chillers for cooling components on the optical table are placed outside of this cleanroom.The laser itself consists of three stages(Fig.1).An almostfinal version of the laser,the so-called engineering prototype,is described in detail in[8].The primary focus of this article is the stabilization and characterization of the PSL.Thus only a rough overview of the laser and the minor modifications implemented between engineering prototype and reference system are given in the following.Thefirst stage,the master laser,is a commercial non-planar ring-oscillator[9,10](NPRO) manufactured by InnoLight GmbH in Hannover,Germany.This solid-state laser uses a Nd:Y AG crystal as the laser medium and resonator at the same time.The NPRO is pumped by laser diodes at808nm and delivers an output power of2W.An internal power stabilization,called the noise eater,suppresses the relaxation oscillation at around1MHz.Due to its monolithic res-onator,the laser has exceptional intrinsic frequency stability.The two subsequent laser stages, used for power scaling,inherit the frequency stability of the master laser.The second stage(medium-power amplifier)is a single-pass amplifier[11]with an output power of35W.The seed laser beam from the NPRO stage passes through four Nd:YVO4crys-#161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10619power stabilizationFig.1.Pre-stabilized laser system of Advanced LIGO.The three-staged laser(NPRO,medium power amplifier,high power oscillator)and the stabilization scheme(pre-mode-cleaner,power and frequency stabilization)are shown.The input-mode-cleaner is not partof the PSL but closely related.NPRO,non-planar ring oscillator;EOM,electro-optic mod-ulator;FI,Faraday isolator;AOM,acousto-optic modulator.tals which are longitudinally pumped byfiber-coupled laser diodes at808nm.The third stage is an injection-locked ring oscillator[8]with an output power of about220W, called the high-power oscillator(HPO).Four Nd:Y AG crystals are used as the active media. Each is longitudinally pumped by sevenfiber-coupled laser diodes at808nm.The oscillator is injection-locked[12]to the previous laser stage using a feedback control loop.A broadband EOM(electro-optic modulator)placed between the NPRO and the medium-power amplifier is used to generate the required phase modulation sidebands at35.5MHz.Thus the high output power and good beam quality of this last stage is combined with the good frequency stability of the previous stages.The reference system features some minor modifications compared to the engineering proto-type[8]concerning the optics:The external halo aperture was integrated into the laser system permanently improving the beam quality.Additionally,a few minor designflaws related to the mechanical structure and the optical layout were engineered out.This did not degrade the output performance,nor the characteristics of the locked laser.In general the PSL is designed to be operated in two different power modes.In high-power mode all three laser stages are engaged with a power of about160W at the PSL output.In low-power mode the high-power oscillator is turned off and a shutter inside the laser resonator is closed.The beam of the medium-power stage is reflected at the output coupler of the high power stage leaving a residual power of about13W at the PSL output.This low-power mode will be used in the early commissioning phase and in the low-frequency-optimized operation mode of Advanced LIGO and is not discussed further in this article.The stabilization has three sections(Fig.1:PMC,PD2,reference cavity):A passive resonator, the so called pre-mode-cleaner(PMC),is used tofilter the laser beam spatially and temporally (see subsection2.1).Two pick-off beams at the PMC are used for the active power stabilization (see subsection2.2)and the active frequency pre-stabilization,respectively(see subsection2.3).In general most stabilization feedback control loops of the PSL are implemented using analog electronics.A real-time computer system(Control and Data Acquisition Systems,CDS,[13]) which is common to many other subsystems of Advanced LIGO,is utilized to control and mon-itor important parameters of the analog electronics.The lock acquisition of various loops,a few #161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10620slow digital control loops,and the data acquisition are implemented using this computer sys-tem.Many signals are recorded at different sampling rates ranging from16Hz to33kHz for diagnostics,monitoring and vetoing of gravitational wave signals.In total four real-time pro-cesses are used to control different aspects of the laser system.The Experimental Physics and Industrial Control System(EPICS)[14]and its associated user tools are used to communicate with the real-time software modules.The PSL contains a permanent,dedicated diagnostic instrument,the so called diagnostic breadboard(DBB,not shown in Fig.1)[15].This instrument is used to analyze two different beams,pick-off beams of the medium power stage and of the HPO.Two shutters are used to multiplex these to the DBB.We are able to measurefluctuations in power,frequency and beam pointing in an automated way with this instrument.In addition the beam quality quantified by the higher order mode content of the beam was measured using a modescan technique[16].The DBB is controlled by one real-time process of the CDS.In contrast to most of the other control loops in the PSL,all DBB control loops were implemented digitally.We used this instrument during the characterization of the laser system to measure the mentioned laser beam parameters of the HPO.In addition we temporarily placed an identical copy of the DBB downstream of the PMC to characterize the output beam of the PSL reference system.2.1.Pre-mode-cleanerA key component of the stabilization scheme is the passive ring resonator,called the pre-mode-cleaner(PMC)[17,18].It functions to suppress higher-order transverse modes,to improve the beam quality and the pointing stability of the laser beam,and tofilter powerfluctuations at radio frequencies.The beam transmitted through this resonator is the output beam of the PSL, and it is delivered to the subsequent subsystems of the gravitational wave detector.We developed and used a computer program[19]to model thefilter effects of the PMC as a function of various resonator parameters in order to aid its design.This led to a resonator with a bow-tie configuration consisting of four low-loss mirrors glued to an aluminum spacer. The optical round-trip length is2m with a free spectral range(FSR)of150MHz.The inci-dence angle of the horizontally polarized laser beam is6◦.Theflat input and output coupling mirrors have a power transmission of2.4%and the two concave high reflectivity mirrors(3m radius of curvature)have a transmission of68ppm.The measured bandwidth was,as expected, 560kHz which corresponds to afinesse of133and a power build-up factor of42.The Gaussian input/output beam had a waist radius of about568µm and the measured acquired round-trip Gouy phase was about1.7rad which is equivalent to0.27FSR.One TEM00resonance frequency of the PMC is stabilized to the laser frequency.The Pound-Drever-Hall(PDH)[20,21]sensing scheme is used to generate error signals,reusing the phase modulation sidebands at35.5MHz created between NPRO and medium power amplifier for the injection locking.The signal of the photodetector PD1,placed in reflection of the PMC, is demodulated at35.5MHz.This photodetector consists of a1mm InGaAs photodiode and a transimpedance amplifier.A piezo-electric element(PZT)between one of the curved mirrors and the spacer is used as a fast actuator to control the round-trip length and thereby the reso-nance frequencies of the PMC.With a maximum voltage of382V we were able to change the round-trip length by about2.4µm.An analog feedback control loop with a bandwidth of about 7kHz is used to stabilize the PMC resonance frequency to the laser frequency.In addition,the electronics is able to automatically bring the PMC into resonance with the laser(lock acquisition).For this process a125ms period ramp signal with an amplitude cor-responding to about one FSR is applied to the PZT of the PMC.The average power on pho-todetector PD1is monitored and as soon as the power drops below a given threshold the logic considers the PMC as resonant and closes the analog control loop.This lock acquisition proce-#161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10621dure took an average of about65ms and is automatically repeated as soon as the PMC goes off resonance.One real-time process of CDS is dedicated to control the PMC electronics.This includes parameters such as the proportional gain of the loop or lock acquisition parameters.In addition to the PZT actuator,two heating foils,delivering a maximum total heating power of14W,are attached to the aluminum spacer to control its temperature and thereby the roundtrip length on timescales longer than3s.We measured a heating and cooling1/e time constant of about2h with a range of4.5K which corresponds to about197FSR.During maintenance periods we heat the spacer with7W to reach a spacer temperature of about2.3K above room temperature in order to optimize the dynamic range of this actuator.A digital control loop uses this heater as an actuator to off-load the PZT actuator allowing compensation for slow room temperature and laser frequency drifts.The PMC is placed inside a pressure-tight tank at atmospheric pressure for acoustic shield-ing,to avoid contamination of the resonator mirrors and to minimize optical path length changes induced by atmospheric pressure variations.We used only low-outgassing materials and fabri-cated the PMC in a cleanroom in order to keep the initial mirror contamination to a minimum and to sustain a high long-term throughput.The PMCfilters the laser beam and improves the beam quality of the laser by suppress-ing higher order transverse modes[17].The acquired round-trip Gouy phase of the PMC was chosen in such a way that the resonance frequencies of higher order TEM modes are clearly separated from the TEM00resonance frequency.Thus these modes are not resonant and are mainly reflected by the PMC,whereas the TEM00mode is transmitted.However,during the design phase we underestimated the thermal effects in the PMC such that at nominal circu-lating power the round-trip Gouy-phase is close to0.25FSR and the resonance of the TEM40 mode is close to that of the TEM00mode.To characterize the mode-cleaning performance we measured the beam quality upstream and downstream of the PMC with the two independent DBBs.At150W in the transmitted beam,the circulating power in the PMC is about6.4kW and the intensity at the mirror surface can be as high as1.8×1010W m−2.At these power levels even small absorptions in the mirror coatings cause thermal effects which slightly change the mirror curvature[22].To estimate these thermal effects we analyzed the transmitted beam as a function of the circulating power using the DBB.In particular we measured the mode content of the LG10and TEM40mode.Changes of the PMC eigenmode waist size showed up as variations of the LG10mode content.A power dependence of the round-trip Gouy phase caused a variation of the power within the TEM40mode since its resonance frequency is close to a TEM00mode resonance and thus the suppression of this mode depends strongly on the Gouy phase.We adjusted the input power to the PMC such that the transmitted power ranged from100W to 150W corresponding to a circulating power between4.2kW and6.4kW.We used our PMC computer simulation to deduce the power dependence of the eigenmode waist size and the round-trip Gouy phase.The results are given in section3.1.At all circulating power levels,however,the TEM10and TEM01modes are strongly sup-pressed by the PMC and thus beam pointingfluctuations are reduced.Pointingfluctuations can be expressed tofirst order as powerfluctuations of the TEM10and TEM01modes[23,24].The PMC reduces thefield amplitude of these modes and thus the pointingfluctuations by a factor of about61according to the measuredfinesse and round-trip Gouy phase.To keep beam point-ingfluctuations small is important since they couple to the gravitational wave channel by small differential misalignments of the interferometer optics.Thus stringent design requirements,at the10−6Hz−1/2level for relative pointing,were set.To verify the pointing suppression effect of the PMC we used DBBs to measure the beam pointingfluctuations upstream and downstream #161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10622Fig.2.Detailed schematic of the power noise sensor setup for thefirst power stabilizationloop.This setup corresponds to PD2in the overview in Fig.1.λ/2,waveplate;PBS,polar-izing beam splitter;BD,glassfilters used as beam dump;PD,single element photodetector;QPD,quadrant photodetector.of the PMC.The resonator design has an even number of nearly normal-incidence reflections.Thus the resonance frequencies of horizontal and vertical polarized light are almost identical and the PMC does not act as polarizer.Therefore we use a thin-film polarizer upstream of the PMC to reach the required purity of larger than100:1in horizontal polarization.Finally the PMC reduces technical powerfluctuations at radio frequencies(RF).A good power stability between9MHz and100MHz is necessary as the phase modulated light in-jected into the interferometer is used to sense several degrees of freedom of the interferometer that need to be controlled.Power noise around these phase modulation sidebands would be a noise source for the respective stabilization loop.The PMC has a bandwidth(HWHM)of about 560kHz and acts tofirst order as a low-passfilter for powerfluctuations with a-3dB corner frequency at this frequency.To verify that the suppression of RF powerfluctuations is suffi-cient to fulfill the design requirements,we measured the relative power noise up to100MHz downstream of the PMC with a dedicated experiment involving the optical ac coupling tech-nique[25].In addition the PMC serves the very important purpose of defining the spatial laser mode for the downstream subsystem,namely the input optics(IO)subsystem.The IO subsystem is responsible,among other things,to further stabilize the laser beam with the suspended input mode cleaner[26]before the beam will be injected into the interferometer.Modifications of beam alignment or beam size of the laser system,which were and might be unavoidable,e.g., due to maintenance,do not propagate downstream of the PMC tofirst order due to its mode-cleaning effect.Furthermore we benefit from a similar isolating effect for the active power and frequency stabilization by using the beams transmitted through the curved high-reflectivity mirrors of the PMC.2.2.Power stabilizationThe passivefiltering effect of the PMC reduces powerfluctuations significantly only above the PMC bandwidth.In the detection band from about10Hz to10kHz good power stability is required sincefluctuations couple via the radiation pressure imbalance and the dark-fringe offset to the gravitational wave channel.Thus two cascaded active control loops,thefirst and second power stabilization loop,are used to reduce powerfluctuations which are mainly caused by the HPO stage.Thefirst loop uses a low-noise photodetector(PD2,see Figs.1and2)at one pick-off port #161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10623of the PMC to measure the powerfluctuations downstream of the PMC.An analog electronics feedback control loop and an AOM(acousto-optic modulator)as actuator,located upstream of the PMC,are used to stabilize the power.Scattered light turned out to be a critical noise source for thisfirst loop.Thus we placed all required optical and opto-electronic components into a box to shield from scattered light(see Fig.2).The beam transmitted by the curved PMC mirror has a power of about360mW.This beam isfirst attenuated in the box using aλ/2waveplate and a thin-film polarizer,such that we are able to adjust the power on the photodetectors to the optimal operation point.Afterwards the beam is split by a50:50beam splitter.The beams are directed to two identical photode-tectors,one for the control loop(PD2a,in-loop detector)and one for independent out-of-loop measurements to verify the achieved power stability(PD2b,out-of-loop detector).These pho-todetectors consist of a2mm InGaAs photodiode(PerkinElmer C30642GH),a transimpedance amplifier and an integrated signal-conditioningfilter.At the chosen operation point a power of about4mW illuminates each photodetector generating a photocurrent of about3mA.Thus the shot noise is at a relative power noise of10−8Hz−1/2.The signal conditioningfilter has a gain of0.2at very low frequencies(<70mHz)and amplifies the photodetector signal in the im-portant frequency range between3.3Hz and120Hz by about52dB.This signal conditioning filter reduces the electronics noise requirements on all subsequent stages,but has the drawback that the range between3.3Hz and120Hz is limited to maximum peak-to-peak relative power fluctuations of5×10−3.Thus the signal-conditioned channel is in its designed operation range only when the power stabilization loop is closed and therefore it is not possible to measure the free running power noise using this channel due to saturation.The uncoated glass windows of the photodiodes were removed and the laser beam hits the photodiodes at an incidence angle of45◦.The residual reflection from the photodiode surface is dumped into a glassfilter(Schott BG39)at the Brewster angle.Beam positionfluctuations in combination with spatial inhomogeneities in the photodiode responsivity is another noise source for the power stabilization.We placed a silicon quadrant photodetector(QPD)in the box to measure the beam positionfluctuations of a low-power beam picked off the main beam in the box.The beam parameters,in particular the Gouy phase,at the QPD are the same as on the power sensing detectors.Thus the beam positionfluctuations measured with the QPD are the same as the ones on the power sensing photodetectors,assuming that the positionfluctuations are caused upstream of the QPD pick-off point.We used the QPD to measure beam positionfluctuations only for diagnostic and noise projection purposes.In a slightly modified experiment,we replaced one turning mirror in the path to the power sta-bilization box by a mirror attached to a tip/tilt PZT element.We measured the typical coupling between beam positionfluctuations generated by the PZT and the residual relative photocurrent fluctuations measured with the out-of-the-loop photodetector.This coupling was between1m−1 and10m−1which is a typical value observed in different power stabilization experiments as well.We measured this coupling factor to be able to calculate the noise contribution in the out-of-the-loop photodetector signal due to beam positionfluctuations(see Subsection3.3).Since this tip/tilt actuator was only temporarily in the setup,we are not able to measure the coupling on a regular basis.Both power sensing photodetectors are connected to analog feedback control electronics.A low-pass(100mHz corner frequency)filtered reference value is subtracted from one signal which is subsequently passed through several control loopfilter stages.With power stabilization activated,we are able to control the power on the photodetectors and thereby the PSL output power via the reference level on time scales longer than10s.The reference level and other important parameters of these electronics are controlled by one dedicated real-time process of the CDS.The actuation or control signal of the electronics is passed to an AOM driver #161737 - $15.00 USD Received 18 Jan 2012; revised 27 Feb 2012; accepted 4 Mar 2012; published 24 Apr 2012 (C) 2012 OSA7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10624。
Structural Biological Materials-Critical Mechanics-Materials Connections
DOI: 10.1126/science.1220854, 773 (2013);339 Science et al.Marc André Meyers ConnectionsStructural Biological Materials: Critical Mechanics-MaterialsThis copy is for your personal, non-commercial use only.clicking here.colleagues, clients, or customers by , you can order high-quality copies for your If you wish to distribute this article to othershere.following the guidelines can be obtained by Permission to republish or repurpose articles or portions of articles): May 1, 2013 (this information is current as of The following resources related to this article are available online at/content/339/6121/773.full.html version of this article at:including high-resolution figures, can be found in the online Updated information and services, /content/339/6121/773.full.html#ref-list-1, 12 of which can be accessed free:cites 51 articles This article/cgi/collection/mat_sci Materials Sciencesubject collections:This article appears in the following registered trademark of AAAS.is a Science 2013 by the American Association for the Advancement of Science; all rights reserved. The title Copyright American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the Science o n M a y 1, 2013w w w .s c i e n c e m a g .o r g D o w n l o a d e d f r o mStructural BiologicalMaterials:CriticalMechanics-Materials ConnectionsMarc AndréMeyers,1,2*Joanna McKittrick,1Po-Yu Chen3Spider silk is extraordinarily strong,mollusk shells and bone are tough,and porcupine quills and feathers resist buckling.How are these notable properties achieved?The building blocks of the materials listed above are primarily minerals and biopolymers,mostly in combination;the first weak in tension and the second weak in compression.The intricate and ingenious hierarchical structures are responsible for the outstanding performance of each material.Toughness is conferred by the presence of controlled interfacial features(friction,hydrogen bonds,chain straightening and stretching);buckling resistance can be achieved by filling a slender column with a lightweight foam.Here,we present and interpret selected examples of these and other biological materials.Structural bio-inspired materials design makes use of the biological structures by inserting synthetic materials and processes that augment the structures’capability while retaining their essential features.In this Review,we explain this idea through some unusual concepts.M aterials science is a vibrant field of in-tellectual endeavor and research.Thisfield applies physics and chemistry, melding them in the process,to the interrela-tionship between structure,properties,and perform-ance of complex materials with technological applications.Thus,materials science extends these rigorous scientific disciplines into complex ma-terials that have structures providing properties and synergies beyond those of pure and simple solids.Initially geared at synthetic materials,ma-terials science has recently extended its reach into biology,especially into the extracellular matrix, whose mechanical properties are of utmost im-portance in living organisms.Some of the semi-nal work and important contributions in this field are either presented or reviewed in(1–5).There are a number of interrelated features that define biological materials and distinguish them from their synthetic counterparts[inspired by Arzt(6)]: (i)Self-assembly.In contrast to many synthetic processes to produce materials,the structures are assembled from the bottom up,rather than from the top down.(ii)Multi-functionality.Many com-ponents serve more than one purpose.For exam-ple,feathers provide flight capability,camouflage, and insulation,whereas bones provide structural framework,promote the growth of red blood cells, and provide protection to the internal organs.(iii) Hierarchy.Different,organized scale levels(nano-to ultrascale)confer distinct and translatable prop-erties from one level to the next.We are starting to develop a systematic and quantitative understandingof this hierarchy by distinguishing the character-istic levels,developing constitutive descriptionsof each level,and linking them through appro-priate and physically based equations,enabling afull predictive understanding.(iv)Hydration.Theproperties are highly dependent on the level ofwater in the structure.There are some exceptions,such as enamel,but this rule applies to mostbiological materials and is of importance to me-chanical properties such as strength(which isdecreased by hydration)and toughness(which isincreased).(v)Mild synthesis conditions.Themajority of biological materials are fabricated atambient temperature and pressure as well as in anaqueous environment,a notable difference fromsynthetic materials fabrication.(vi)Evolution andenvironmental constraints.The limited availabil-ity of useful elements dictates the morphologyand resultant properties.The structures are notnecessarily optimized for all properties but arethe result of an evolutionary process leading tosatisfactory and robust solutions.(vii)Self-healingcapability.Whereas synthetic materials undergodamage and failure in an irreversible manner,biological materials often have the capability,due to the vascularity and cells embedded in thestructure,to reverse the effects of damage byhealing.The seven characteristics listed above arepresent in a vast number of structures.Nevertheless,the structures of biological materials can bedivided into two broad classes:(i)non-mineralized(“soft”)structures,which are composed of fibrousconstituents(collagen,keratin,elastin,chitin,lignin,and other biopolymers)that display widelyvarying mechanical properties and anisotropiesdepending on the function,and(ii)mineralized(“hard”)structures,consisting of hierarchicallyassembled composites of minerals(mainly,butnot solely,hydroxyapatite,calcium carbonate,and amorphous silica)and organic fibrous com-ponents(primarily collagen and chitin).The mechanical behavior of biological con-stituents and composites is quite diverse.Bio-minerals exhibit linear elastic stress-strain plots,whereas the biopolymer constituents are non-linear,demonstrating either a J shape or a curvewith an inflection point.Foams are characterizedby a compressive response containing a plastic orcrushing plateau in which the porosity is elim-inated.Many biological materials are compositeswith many components that are hierarchicallystructured and can have a broad variety of con-stitutive responses.Below,we present some of thestructures and functionalities of biological ma-terials with examples from current research.Here,we focus on three points:(i)How high tensilestrength is achieved(biopolymers),(ii)how hightoughness is attained(composite structures),and(iii)how bending resistance is achieved in light-weight structures(shells with an interior foam).Structures in Tension:Importance of BiopolymersThe ability to sustain tensile forces requires aspecific set of molecular and configurational con-formations.The initial work performed on exten-sion should be small,to reduce energy expenditure,whereas the material should stiffen close to thebreaking point,to resist failure.Thus,biopolymers,such as collagen and viscid(catching spiral)spidersilk,have a J-shaped stress-strain curve where mo-lecular uncoiling and unkinking occur with con-siderable deformation under low stress.This stiffening as the chains unfurl,straighten,stretch,and slide past each other can be repre-sented analytically in one,two,and three dimen-sions.Examples are constitutive equations initiallydeveloped for polymers by Ogden(7)and Arrudaand Boyce(8).An equation specifically proposedfor tissues is given by Fung(3).A simpler for-mulation is given here;the slope of the stress-strain(s-e)curve increases monotonically with strain.Thus,one considers two regimes:(i)unfurlingand straightening of polymer chainsd sd eºe nðn>1Þð1Þand(ii)stretching of the polymer chain backbonesd sd eºEð2Þwhere E is the elastic modulus of the chains.Thecombined equation,after integrating Eqs.1and2,iss=k1e n+1+H(e c)E(e–e c)(3)Here k1is a parameter,and H is the Heavisidefunction,which activates the second term at e=e c,where e c is a characteristic strain at whichcollagen fibers are fully extended.Subsequent straingradually becomes dominated by chain stretch-ing.The computational results by Gautieri et al.(9)on collagen fibrils corroborate Eq.3for n=1.This corresponds to a quadratic relation between1Department of Mechanical and Aerospace Engineering andMaterials Science and Engineering Program,University ofCalifornia,San Diego,La Jolla,CA92093,USA.2Department ofNanoengineering,University of California,San Diego,La Jolla,CA92093,USA.3Department of Materials Science and En-gineering,National Tsing Hua University,Hsinchu30013,Taiwan,Republic of China.*To whom correspondence should be addressed.E-mail:mameyers@ SCIENCE VOL33915FEBRUARY2013773o n M a y 1 , 2 0 1 3 w w w . s c i e n c e m a g . o r g D o w n l o a d e d f r o mstress and strain (s ºe 2),which has the char-acteristic J shape.Collagen is the most important structural bio-logical polymer,as it is the key component in many tissues (tendon,ligaments,skin,and bone),as well as in the extracellular matrix.The de-formation process is intimately connected to the different hierarchical levels,starting with the poly-peptides (0.5-nm diameter)to the tropocollagen molecules (1.5-nm diameter),then to the fibrils (~40-to 100-nm diameter),and finally to fibers (~1-to 10-m m diameter)and fascicles (>10-m m diameter).Molecular dynamics computations (9)of entire fibrils show the J -curve response;these computational predictions are well matched to atomic force microscopy (AFM)(10),small-angle x-ray scattering (SAXS)(11),and experiments by Fratzl et al .(12),as shown in Fig.1A.The effect of hydration is also seen and is of great impor-tance.The calculated density of collagen de-creases from 1.34to 1.19g/cm 3with hydration and is accompanied by a decrease in the Young ’s modulus from 3.26to 0.6GPa.The response of silk and spider thread is fascinating.As one of the toughest known ma-terials,silk also has high tensile strength and extensibility.It is composed of b sheet (10to 15volume %)nanocrystals [which consist of highly conserved poly-(Gly-Ala)and poly-Ala domains]embedded in a disordered matrix (13).Figure 1B shows the J -shape stress-strain curve and molecular configurations for the crystalline domains in silkworm (Bombyx mori )silk (14).Similar to collagen,the low-stress region corre-sponds to uncoiling and straightening of the pro-tein strands.This region is followed by entropic unfolding of the amorphous strands and then stiffening due to load transfer to the crystalline b sheets.Despite the high strength,the major mo-lecular interactions in the b sheets are weak hy-drogen bonds.Molecular dynamics simulations,Fig.1.Tensile stress-strain relationships in bio-polymers.(A )J -shaped curve for hydrated and dry collagen fibrils obtained from molecular dynamics (MD)simulations and AFM and SAXS studies.At low stress levels,considerable stretching occurs due to the uncrimping and unfolding of molecules;at higher stress levels,the polymer backbone stretches.Adapted from (9,12).(B )Stretching of dragline spider silk and molecular schematic of the protein fibroin.At low stress levels,entropic effects domi-nate (straightening of amorphous strands);at higher levels,the crystalline parts sustain the load.(C )Mo-lecular dynamics simulation of silk:(i)short stack and (ii)long stack of b -sheet crystals,showing that a higher pullout force is required in the short stack;for the long stack,bending stresses become im-portant.Hydrogen bonds connect b -sheet crystals.Adapted from (14).(D )Egg whelk case (bioelastomer)showing three regions:straightening of the a helices,the a helix –to –b sheet transformation,and b -sheet extension.A molecular schematic is shown.Adapted from (18).300.000.2Yield pointEntropic unfoldingMD simulationsStick slipStiffening β-crystal123456700012345670102030405050010001500200025050075010001250150017500.40.60.80.010.020.030.040.05MD wet (Gautieri et al)SAXS (Sasaki and Odajima)AFM (Aladin et al)MD dry (Gautieri et al)2520151050S t r e s s (M P a )S t r e s s(M P a )StrainABCDStrain (m/m)Length (nm)Length (nm)Stick-slip deformation (robust)"brittle" fracture (fragile)i iiP u l l -o u t f o r c e (p N )00.20.4Native state Unloading: reformation of α-helicesDomain 4: Extension andalignmentof β-sheets0.60.8ε=0ε4ε=01.0012345StrainS t r e s s (M P a )E n e r g y /v o l u m e (k c a l /m o l /n m 3)L e n g t hI I II II III IIIIVIVFDomain 3: Formation of β-sheetsfrom random coilsε3Domain 2: Extension of random coilsε2Domain 1: Unraveling of α-helicesinto random coilsε1Toughness (MD)Resilience (MD)T=-1°C T=20°C T=40°C T=60°C T=80°C15FEBRUARY 2013VOL 339SCIENCE 774REVIEWo n M a y 1, 2013w w w .s c i e n c e m a g .o r g D o w n l o a d e d f r o mshown in Fig.1C,illustrate an energy dissipative stick-slip shearing of the hydrogen bonds during failure of the b sheets (14).For a stack with a height L ≤3nm (left-hand side of Fig.1C),the shear stresses are more substantial than the flex-ure stresses,and the hydrogen bonds contribute to the high strength obtained (1.5GPa).How-ever,if the stack of b sheets is too high (right-hand side of Fig.1C),it undergoes bending with tensile separation between adjacent sheets.The nanoscale dimension of the b sheets allows for a ductile instead of brittle failure,resulting in high toughness values of silk.Thus,size affects the mechanical response considerably,changing the deformation characteristics of the weak hydro-gen bonds.This has also been demonstrated in bone (15–17),where sacrificial hydrogen bonds between mineralized collagen fibrils contribute to the excellent fracture resistance.Other biological soft materials have more complex responses,marked by discontinuities in d s /d e .This is the case for wool,whelk eggs,silks,and spider webs.Several mechanisms are responsible for this change in slope;for instance,the transition from a -to b -keratin,entropic changes with strain (such as those prevalent in rubber,where chain stretching and alignment decrease entropy),and others.The example of egg whelk is shown in Fig.1D (18).In this case,there is a specific stress at which a -keratin heli-ces transform to b sheets,with an associated change in length.Upon unloading,the reverse occurs,and the total reversible strain is,therefore,extensive.This stress-induced phase transforma-tion is similar to what occurs in shape-memory alloys.Thus,this material can experience sub-stantial reversible deformation (up to 80%)in a reversible fashion,when the stress is raised from 2to 5MPa,ensuring the survival of whelk eggs,which are continually swept by waves.These examples demonstrate the distinct properties of biopolymers that allow these ma-terials to be strong and highly extensible with distinctive molecular deformation characteristics.However,many interesting biological materials are composites of flexible biopolymers and stiff minerals.The combination of these two constit-uents leads to the creation of a tough material.Imparting Toughness:Importance of Interfaces One hallmark property of most biological com-posites is that they are tough.Toughness is defined as the amount of energy a material ab-sorbs before it fails,expressed asU ¼∫e fs d eð4Þwhere U is the energy per volume absorbed,s is the stress,e is the strain,and e f is the failure strain.Tough materials show considerable plastic deformation (or permanent damage)coupled with considerable strength.This maximizes the integral expression in Eq.4.Biological com-posite materials (for example,crystalline and noncrystalline components)have a plethora oftoughening mechanisms,many of which depend on the presence of interfaces.As a crack im-pinges on an interface or discontinuity in the material,the crack can be deflected around the interface (requiring more energy to propagate than a straight crack)or can drive through it.The strength of biopolymer fibers in tension im-pedes crack opening;bridges between micro-cracks are another mechanism.The toughening mechanisms have been divided into intrinsic (ex-isting in the material ahead of crack)and extrinsic (generated during the progression of failure)cat-egories (19).Thus,toughening is accomplished by a wide variety of stratagems.We illustrate this concept for four biological materials,shown in Fig.2.All inorganic materials contain flaws and cracks,which reduce the strength from the theo-retical value (~E /10to E /30).The maximum stress (s max )a material can sustain when a preexisting crack of length a is present is given by the Griffith equations max ¼ffiffiffiffiffiffiffiffiffiffi2g s E p a r ¼YK Icffiffiffiffiffip ap ð5Þwhere E is the Young ’s modulus,g s is the sur-face (or damage)energy,and Y is a geometric parameter.K Ic ¼Y −1ffiffiffiffiffiffiffiffiffiffi2g s E p is the fracture toughness,a materials property that expresses the ability to resist crack propagation.Abalone (Haliotis rufescens )nacre has a fracture tough-ness that is vastly superior to that of its major constituent,monolithic calcium carbonate,due to an ordered assembly consisting of mineral tiles with an approximate thickness of 0.5m m and a diameter of ~10m m (Fig.2A).Additionally,this material contains organic mesolayers (separated by ~300m m)that are thought to be seasonal growth bands.The tiles are connected by mineral bridges with ~50-nm diameter and are separated by organic layers,consisting of a chitin network and acidic proteins,which,when combined,have a similar thickness to the mineral bridge diame-ters.The Griffith fracture criterion (Eq.5)can be applied to predict the flaw size (a cr )at which the theoretical strength s th is achieved.With typical values for the fracture toughness (K Ic ),s th ,and E ,the critical flaw size is in the range of tens of nanometers.This led Gao et al .(20)to propose that at sufficiently small dimensions (less than the critical flaw size),materials become insensitive to flaws,and the theoretical strength (~E /30)should be achieved at the nanoscale.However,the strength of the material will be determined by fracture mechanisms operating at all hierar-chical levels.The central micrograph in Fig.2A shows how failure occurs by tile pullout.The interdigitated structure deflects cracks around the tiles instead of through them,thereby increasing the total length of the crack and the energy needed to fracture (increasing the toughness).Thus,we must de-termine how effectively the tiles resist pullout.Three contributions have been identified and are believed to operate synergistically (21).First,themineral bridges are thought to approach thetheoretical strength (10GPa),thereby strongly attaching the tiles together (22).Second,the tile surfaces have asperities that are produced during growth (23)and could produce frictional resist-ance and strain hardening (24).Third,energy is required for viscoelastic deformation (stretching and shearing)of the organic layer (25).One important aspect on the mechanical prop-erties is the effect of alignment of the mineral crystals.The oriented tiles in nacre result in an-isotropic properties with the strength and modulus higher in the longitudinal (parallel to the organic layers)than in the transverse direction.For a composite with a dispersed mineral m of volume fraction V m embedded in a biopolymer (bp)matrix that has a much lower strength and Young ’s modulus than the mineral,the ratio of the lon-gitudinal (L)and transverse (T)properties P (such as elastic modulus)can be expressed,in simpli-fied form,asP L P T ¼P mP bpV m ð1−V m Þð6ÞThus,the longitudinal properties are much higher than the transverse properties.This aniso-tropic response is also observed in other oriented mineralized materials,such as bone and teeth.Another tough biological material is the exo-skeleton of an arthropod.In the case of marine animals [for instance,lobsters (26,27)and crabs (28)],the exoskeleton structure consists of layers of mineralized chitin in a Bouligand arrange-ment (successive layers at the same angle to each other,resulting in a helicoidal stacking sequence and in-plane isotropy).These layers can be en-visaged as being stitched together with ductile tubules that also perform other functions,such as fluid transport and moisture regulation.The cross-ply Bouligand arrangement is effective in crack stopping;the crack cannot follow a straight path,thereby increasing the materials ’toughness.Upon being stressed,the mineral components frac-ture,but the chitin fibers can absorb the strain.Thus,the fractured region does not undergo physical separation with dispersal of fragments,and self-healing can take place (29).Figure 2B shows the structure of the lobster (Homarus americanus )exoskeleton with the Bouligand ar-rangement of the fibers.Bone is another example of a biological ma-terial that demonstrates high toughness.Skeletal mammalian bone is a composite of hydroxyapatite-type minerals,collagen and water.On a volu-metric basis,bone consists of ~33to 43volume %minerals,32to 44volume %organics,and 15to 25volume %water.The Young ’s modulus and strength increase,but the toughness decreases with increasing mineral volume fraction (30).Cortical (dense)mammalian bone has blood ves-sels extending along the long axis of the limbs.In animals larger than rats,the vessel is encased in a circumferentially laminated structure called the osteon.Primary osteons are surrounded by hypermineralized regions,whereas secondary SCIENCEVOL 33915FEBRUARY 2013775REVIEWo n M a y 1, 2013w w w .s c i e n c e m a g .o r g D o w n l o a d e d f r o m(remodeled)osteons are surrounded by a cement line (also of high mineral content)(31).In mam-malian cortical bone,the following intrinsic toughening mechanisms have been identified:molecular uncoiling and intermolecular sliding of collagen,fibrillar sliding of collagen bonds,and microcracking of the mineral matrix (19).Extrinsic mechanisms are collagen fibril bridging,uncracked ligament bridging,and crack deflec-tion and twisting (19).Rarely does a limb bone snap in two with smooth fracture surfaces;the crack is often deflected orthogonal to the crack front direction.In the case of (rehydrated)elk (Cervus elaphus )antler bone (shown in Fig.2C)(32),which has the highest toughness of any bone type by far (33),the hypermineralized re-gions around the primary osteons lead to crackdeflection,and the high amount of collagen (~60volume %)adds mechanisms of crack re-tardation and creates crack bridges behind the crack front.The toughening effect in antlers has been estimated as:crack deflection,60%;un-cracked ligament bridges,35%;and collagen as well as fibril bridging,5%(33).A particu-larly important feature in bone is that the fracture toughness increases as the crack propagates,as shown in the plot.This plot demonstrates the crack extension resistance curve,or R -curve,behavior,which is the rate of the total energy dissipated as a function of the crack size.This occurs by the activation of the extrinsic tough-ening mechanisms.In this manner,it becomes gradually more difficult to advance the crack.In human bone,the cracks are deflected and/ortwisted around the cement lines surrounding the secondary osteons and also demonstrate R -curve behavior (34).The final example illustrating how the presence of interfaces is used to retard crack propagation is the glass sea sponge (Euplectella aspergillum ).The entire structure of the V enus ’flower basket is shown in Fig.2D.Biological silica is amorphous and,within the spicules,consists of concentric layers,separated by an organic material,silicatein (35,36).The flexure strength of the spicule notably exceeds (by approximately fivefold)that of monolithic glass (37).The principal reason is the presence of interfaces,which can arrest and/or deflect the crack.Biological materials use ingenious meth-ods to retard the progression of cracks,therebyAbalone shell: NacreMineral bridgesLobsterDeer antlerChitin fibril networkHuman cortical boneMineral crystallitesPrimary osteonsSubvelvet/compact Subvelvet/cCompact Comp p actTransition zoneCancellousCollagen fibrilsDeep sea spongeSkeletonSpicules20 mm1 cmHuman cortical boneElk antlerTransverseIn-plane longitudinalASTM validASTM invalid Mesolayers ABCD0.1 mm500 nm500 nm ˜1 nm˜3 nm˜20 nmCrack extension, ⌬a (mm)T o u g h n e s s , J (k J m -2)50 nm200 nm 10 m500 nm2 m1 m200 m300 m˜10 m0.010.11101000.20.40.6500 00 nm50 nmFig.2.Hierarchical structures of tough biological materials demonstrating the heterogeneous interfaces that provide crack deflection.(A )Abalone nacre showing growth layers (mesolayers),mineral bridges between mineral tiles and asperities on the surface,the fibrous chitin network that forms the backbone of the inorganic layer,and an example of crack tortuosity in which the crack must travel around the tiles instead of through them [adapted from (4,21)].(B )Lobster exoskeleton showing the twisted plywood structure of the chitin (next to the shell)and the tubules that extend from the chitin layers to the animal [adapted from (27)].(C )Antler bone image showing the hard outer sheath (cortical bone)surrounding the porous bone.The collagen fibrils are highly aligned in the growth direction,with nanocrystalline minerals dispersed in and around them.The osteonal structure in a cross section of cortical bone illustrates the boundaries where cracks perpendicular to the osteons can be directed [adapted from (33)].ASTM,American Society for Testing and Mate-rials.(D )Silica sponge and the intricate scaffold of spicules.Each spicule is a circumferentially layered rod:The interfaces between the layers assist in ar-resting crack anic silicate in bridging adjacent silica layers is observed at higher magnification (red arrow)(36).15FEBRUARY 2013VOL 339SCIENCE776REVIEWo n M a y 1, 2013w w w .s c i e n c e m a g .o r g D o w n l o a d e d f r o mincreasing toughness.These methods operate at levels ranging from the nanoscale to the structur-al scale and involve interfaces to deflect cracks,bridging by ductile phases (e.g.,collagen or chitin),microcracks forming ahead of the crack,delocal-ization of damage,and others.Lightweight Structures Resistant to Bending,Torsion,and Buckling —Shells and FoamsResistance to flexural and torsional tractions with a prescribed deflection is a major attribute of many biological structures.The fundamental mechanics of elastic (recoverable)deflection,as it relates to the geometrical characteristics of beams and plates,is given by two equations:The first relates the bending moment,M ,to the curvature of the beam,d 2y /dx 2(y is the deflection)d 2y dx 2¼MEIð7Þwhere I is the area moment of inertia,which de-pends on the geometry of the cross section (I =p R 4/4,for circular sections,where R is the ra-dius).Importantly,the curvature of a solid beam,and therefore its deflection,is inversely propor-tional to the fourth power of the radius.The sec-ond equation,commonly referred to as Euler ’s buckling equation,calculates the compressive load at which global buckling of a column takes place (P cr )P cr ¼p 2EI ðkL Þ2ð8Þwhere k is a constant dependent on the column-end conditions (pinned,fixed,or free),and L is the length of the column.Resistance to buck-ing can also be accomplished by increasing I .Both Eqs.7and 8predict the principal designLongitudinal sectionToucan beak Keratin layers(i) Fibers(circumferential)Megafibrils and fibrilsBarbsBarbulesCortexCortical ridgesFoamRachisNodes(iii) Medulloidpith(ii) Fibers (longitudinal)Feather rachisPlant-Bird of ParadisePorcupine quillsNodesRebarClosed-cell foamTransverseLongitudinalCross sectionABCD5 mm 1 mm1 cm 0.1 mm5m 5 m m1c 1 c m1 mm100 m500 mFig.3.Low-density and stiff biological materials.The theme is a dense outer layer and a low-density core,which provides a high bending strength –to –weight ratio.(A )Giant bird of paradise plant stem showing the cellular core with porous walls.(B )Porcupine quill exhibiting the dense outer cortex surrounding a uniform,closed-cell foam.Taken from (42).(C )Toucan beak showing the porousinterior (bone)with a central void region [adapted from (43)].(D )Schematic view of the three major structural components of the feather rachis:(i)superficial layers of fibers,wound circumferentially around the rachis;(ii)the majority of the fibers extending parallel to the rachidial axis and through the depth of the cortex;and (iii)foam comprising gas-filled polyhedral structures.Taken from (45)SCIENCEVOL 33915FEBRUARY 2013777REVIEWo n M a y 1, 2013w w w .s c i e n c e m a g .o r g D o w n l o a d e d f r o m。