王瑞利等： 多物理耦合非线 性偏微分方程与数值解不确定度量化数学方法 身， 还是与之耦合的常微分方程、复杂函数关系及其数值求解过程中都含有未认知或随机 的变量， 使 得耦合 的微分方程组中含有大量不确定性的因素． 要通过建模与模拟 ( mo d e l i n g a n d s i mu l a t io n ， M&S ) 这条途径认证或决策大型工程设计参数或系统的可靠性， 必须对其建模与模拟中的不确定度进行量 化， 才能发展高可信度的数值模拟仿真软件． 不确定度量化 ( u n c e r t a in t y q u a n t i fi c a t i o n ， uo) 的数学方 法是评估 M& S可信度 的最佳途径， 在航空航天和国防安全等 国家众多尖端领域， 己引起学术界和工 程应用领域的高度重视． 2 多物理耦合非线性偏微分方程的数学模型 2 ． 1 非 线性 双曲型偏 微 分方栏 ( 非足 帚爆=轰流1 本刀字模 型) 这里以爆轰流体力学为例， 炸药爆轰过程所使用的基本方程是非定常可压缩流体力学方程和化学 反应动力学方程的藕合方程组， + ：0 ( 质量方程) ， ( 2 ．1 ) + p ~+vP=0 ( 动量方程) ， ( 2 -2 ) +V- p E + ． P =0 ( 能量方程) ， ( 2 ．3 ) P ： { 三 ， 、 药 ’ (状 态 方 程 )， (2 _4) I P (p ，e ， F) ， 炸药 = P( p ， e ， ) ( 炸药反应率) ， ( 2 ．5 ) 其中 E=e - 4- Ⅱ 。
， P表示密度， 表示速度， E， e和 P分别表示单位质量的总能量、内能和压力， 为 化学反应份额 ( 燃烧 函数) ， 这里采用炸药反应率的唯象模型． 2 ．2 一阶常微分方程 ( 炸药反应率唯象模型) 由于爆轰过程中化学反应过程复杂， 反应速度快， 反应过程与放能物质的特性 、状态和放能引起 的物质运动密切相关．目前， 从理论上， 严格建立反应率方程是相当困难 的， 只能采用唯象近似．国内 外研究了多种形式的唯象反应率模型， 这些唯象模型往往是一阶的常微分方程， 这里给出几种常用的 唯象模型： ( 1 ) A r r h e n i u s 反应率模型． 其形式为 一 c - ， ( 2 ．6 ) 其中 叫是未反应炸药的质量分数， ， E 和 R是常数， 是温度 ( 2 ) F o r e s t F ir e反应率模型． 其形式为 一 _e c ， 其中 W是未反应炸药的质量分数， P是压力， a t 是常数， 由实验确定 7 2 4 ( 2 ． 7 ) 中 国科学 ： 数学第 4 5卷第 6期 ( 3 ) C o c h r a n反应率模型． 其形式为 R = d F = ( 1 -F ) 1 P n - ~- W 2 B Y ] ， 其中 F是反应产物的质量分数， P是压力， W1 ， W 2和 n是常数， 由实验来定 ( 4 ) L e e 反应率模型． 其形式为 OF = ( 1 -F) 叩 +G( 1 一F) F P ， ( 2 ． 8 ) ( 2 ． 9 ) 其 中 F 是己反应炸药的分数， 叼= vv o 一1 ， V o是炸药的初始比容， V 是受冲击但未反应的炸药比容， P 是压力， I ， X ， G， Y ， Z和 均为常数． ( 5 ) Wi l k i n s反应率模型． 将时间燃烧 函数和 C — J比容燃烧 函数组合起来． F=0为凝固炸药 ( 未 反应)区； 0F1为反应区； F=1为爆炸产物区． 其形式为 F=[ma x { F 1 ， F 2 ) ] ， 其中 为 C — J比容燃烧函数， F 2为时间燃烧函数， T t b 是可调参数 ( 2 ． 1 0 ) f 0 ， 凝 固 炸 药 区 ， F1= { ， 燃区 ， (2 ．1 1) 【 1 ， 爆炸产物 区， f 0 ， ≤ ， F 2 = { ≤ + ， (2．12 ) 【 1 ，t≥ t b + △L， 这里 =藉 是c — J 比 容； t b 是 爆轰 波刚到 达计算网 格的时 刻， 即 起爆时间 ； t 是当 前计 算时刻； △ = r 6 AR ，AR是网格宽度， Dj是爆轰速度， 是可调参数． 2 ． 3 函数关系式 ( 物态方程唯象模型) 在炸药爆轰驱动装置中， 涉及炸药和非炸药 ( 一般为金属)两类材料， 这些材料随着状态的变化， 其行为相当复杂．目前， 大部分采用平衡态下系统的温度与状态参量之间的函数关系式描述， 即物态 方程．针对某种材料， 可能有多种形式的函数关系式描述， 这里给 出炸药爆炸产物和金属采用的物态 方程唯象模型． ( 1 ) 炸药． 爆炸产物常采用 J wL( J o n e s — Wi l k i n s ． L e e ) 形式的状态方程， 其形式为 P = ( 一 ) e- R ~V + B ( 一 ) e- R：V + w E ，(2．13 ) 其中 V： ， E为 比热容力学能， ， B， R1 ， R 2和 W 是可调参数， u和 是 比容． 或爆轰产物的状态 方程采用理想气体： P=( —1 ) p e ， 用燃烧函数 F把凝固炸药和爆轰产物的状态方程连接起来， 即用 P=f 一1 ) p e F． 7 2 5 王瑞利等：多物理耦合非线性偏微分方程与数值解不确定度量化数学方法 在 J wL状态方程中， 其可调参数 A， B， Rl ， R2和 埘之问又有相关性．由爆轰波阵面上的守恒关 系和 C ． J条件， 可推出 A ， B 和 C与 R 1和 R2 相关性的线性方程组， 即给定一组 R 1和 R2 ， 可以求出 ，B 和 ． A e - R 1 +B e - R 2 +c v 2 + =P j ， A Rl e - R 1 +BB ：2 e - R ~ + ( +1 ) V j - ( + ) =p o D ， A e - R 1Vj + ~ e- R ~V j + C v - w = 加+ 去 (1 一 ) ( 2 ． 1 4 ) ( 2 ．1 5 ) ( 2 ．1 6 ) 其 中y j = 岳 ， D是 爆 速 由 实 验 确 定 ， P j 和y j 分 别 是 炸 药c — J 状 态 下 的 爆 压 和 比 容 ， Q 是 爆 热 ， A ， B 和 C三个参数依赖于不确定性参数 R1和 R2 ． ( 2 ) 金属． 金属材料可采用理想气体： P=( 7—1 ) p e ， 也可采用更接近实际的复杂函数关系式， 如 G r u n n e i s e n形式的状态方程 其形式为 P = ( 1 一 r ) + F p (e - eo )， P _ 1 j (2 ．17 ) ： ， 若 ， ( 2 ． 1 8 ) ：{ [1 一 ( 一 1 ) ] ’ ～’ ( ) I P o ， 若 0 ， 其中 P 0 ， e 0和 分别为材料的物性参数， 和 r为常数． 从爆轰流体力学数学模型可看 出， 模型中涉及初边值条件、 方程形式和参数等多种不确定性因素． 3 多物理耦合数学模型的数值求解方法 多物理耦合数学模型的数值求解方法， 一般涉及计算空间离散 ( 即网格) 、方程时空离散和离散方 程的求解算法． 3 ． 1 非定常爆轰流体力学计算方法 爆轰波 由冲击波和波后的化学反应 区组成， 是一个非常强的冲击波． 对于强间断 ( 即冲击波) 的计 算有两种方法 [1 o ] ．一种是激波装配法， 将间断面看成是块块连续解的边界， 在连续解区域内用差分方 法求解， 而在间断面上用 Ra n k i n e — Hu g o n io t条件作为 内边界条件． 另一种是激波捕捉法， 在差分格式 中引入人为黏性将间断解光滑化． 激波捕捉法的思想和做法可归纳为三点： f 1 ) 在动量方程和能量方 程 的压力项 中加上人为黏性， 这就是流体运动中引进 了某种人为耗散机制， 使得冲击波间断解变成一 个在相当狭窄的过渡区域内急剧变化 的 但却是连续的解； ( 2 1 要求所加的人为黏性项只起到使冲击波 间断光滑化 的作用， 而基本上不会影响冲击波过渡区域以外的连续解 的计算结果， 可以近似地满足冲 击波间断条件； ( 3 ) 冲击波过渡区的范围应限制在几个空间步长 以内， 并且这个过渡区在计算过程中不 会扩大， 而过渡 区的速度应逼近真实的冲击波速度． v o n N e u ma n n和 Ric h t m y e r引入人为黏性 q ， 他们 在动量方程 ( 2 ． 2 ) 和能量方程 ( 2 ． 3 )中将压力 P换成 P=P+q ． 人为黏性 q的经典形式为 q =q N R+钆， ( 3 ．1 ) 其中 q N R为 v o n Ne u ma n n — Ric h t my e r 人为黏性 ( 二次黏性) ， q L为 L a n d s h o ff人为黏性 ( 一次黏性或线 性黏性1 ． 7 2 6 中国科学 ： 数学第 4 5卷第 6期 ( 1 ) y o n N e u ma n n — R i c h t my e r 人为黏性 ( 二次黏性) ： q ⅣR ： z R ( ) ， 1 0 ) 和非线性问题． 对高维、小样本和强非线性等复杂 问题， 经过不断努 力发展了很多非参数型的拟合方法， 如 Ga u s s响应面方法、多元自适应回归样条方法 ( MA RS ) 、人工 神经网络 f A NN ) 和支持向量机 f S V M) ， 尤其是 G a u s s 过程响应面还能够考虑响应面模型本身的不确 定性， 因此， 在复杂工程建模与模拟 ( MgS ) ， 尤其是结构力学 M＆s的不确定性量化研究中受到广泛的 关注 ．响应面方法存在的问题是， 逼近 的响应面是否是最佳逼近以及以什么速度收敛到最佳逼近， 这 是具有挑战性的问题． 5 ． 4 随机微分方程 随机微分方程是 2 0世纪中叶发展起来的一个不确定性量化方法．广义上的随机微分方程 [4 0 , 4 1 ] f s t o c h a s t i c d i ff e r e n t i a l e q u a t i o n s 1 包括 由随机过程驱动的微分系统和系数为随机量的微分方程， 狭义 上的随机微分方程常常专指前者． 对于由随机过程驱动的微分系统， 方程 的解一般是一个随机过程函 数， 一般求解微分方程的方法不适用于求解随机微分方程． 然而， 有些复杂过程可以直接用随机过程建 模， 模型的不确定性直接表现在模型本身上而不是输入参数的不确定性， 随机微分方程是一种研究不 确定性量化的很好途径．目前随机微分方程 的分析及数值求解是应用数学界的热点和前沿方向， 特别 是近三十年来随着随机微分方程越来越广泛地应用于系统科学、 工程控制、 生态学和复杂工程等各个 方面， 使得随机微分方程的理论和应用有了迅速的发展， 内容十分丰富． 因此， 结合实际应用中的典型 微分方程， 发展相应的随机微分方程模型及高效求解方法是当前和今后的重要研究方向． 系数为随机 量的微分方程可以直接处理为随机参数的微分方程， 常见的分析方法是随机谱方法 ( s p e c t r a l me t h o d ) ． 随机谱方法 【4 2 — 4 7 ] 是使用谱逼近随机微分方程或参量， 然后将其分解为独立的确定性分量和随机分量 来量化不确定性． 其数学思想是把数学模型中的每个参量利用正交多项式 f 如 H e r mit e多项式) 展开 成无穷级数项， 在实际应用中取有限项． 无穷级数第一项表示参量的平均值， 第二项表示 Ga u s s随机 波动， 第三项及高阶项表示非 Ga u s s随机波动． 近几十年来， 随机谱方法在建模与模拟不确定性量化 中得到 了广泛应用．常见的两类方法：K L展开 f K a r 。
