CN112231946B - 一种基于最优权重因子的激光烧蚀高精度数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于最优权重因子的激光烧蚀高精度数值模拟方法，适用于激光与物质相互作用的数值模拟技术领域。本方法旨在利用有限差分法数值求解热传导偏微分方程，从而模拟靶材物质在激光烧蚀作用下的温度场分布演化。有限差分加权隐格式中一般常用的权重因子取值为0、1或0.5，而本方法则根据当前问题中具体的方程参数和离散时空网格参数，结合线性拟合迭代与非线性拟合方法来搜寻专门针对当前问题的最优权重因子值。本发明搜寻最优权重因子的方法简捷高效、准确度高，而基于搜寻得到的最优权重因子来求解偏微分方程，其数值计算结果的精确度能够显著高于基于一般常用权重因子的结果。本发明适用于激光脉冲烧蚀固体靶材的数值模拟。

Description

一种基于最优权重因子的激光烧蚀高精度数值模拟方法

技术领域

本发明涉及激光与物质相互作用的数值模拟技术领域，具体而言是一种基于最优权重因子的激光烧蚀高精度数值模拟方法，该方法最主要的特征是计算结果的精确度能够显著高于基于一般常用权重因子的结果。

背景技术

激光与物质相互作用(Laser-Matter Interaction)可以带来各种物理、化学、生物等效应，因此有着广泛的应用领域，例如：材料加工(表面热处理、切割焊接等)、材料制备(沉积薄膜、微纳材料等)、光谱分析(激光诱导击穿光谱、激光吸收光谱等)、生物医学(诊断、手术)等等。激光烧蚀是激光与固体物质相互作用最常见的物理效应之一，因此靶材物质在激光加热作用下的温度场分布演化是人们研究的热点。数值模拟是真实实验、解析理论之外的第三类重要科研方法，近几十年来随着计算机技术的不断进步而蓬勃发展，也在激光与物质相互作用领域得到了广泛而深入的应用。对于研究靶材温度场分布演化而言，其核心是求解热传导方程。热传导方程的数学形式通常是二阶偏微分方程，而求解此类方程最常见的数值模拟方法之一就是有限差分法。

有限差分法的基本原理是用有限的离散化时空网格来近似表示连续的时间和空间，从而将连续分布的核心物理量进行离散化。对于热传导偏微分方程来说，核心物理量就是温度，有限差分法利用离散网格点的差分来近似替代数学上的微分，从而将热传导偏微分方程转化为以各网格点温度值为未知量的线性代数方程组，并得到各网格点温度的数值解。用不同方式定义差分可以得到不同的差分格式，一般常用的有向前欧拉格式、向后欧拉格式、Crank-Nicolson格式等。在这三种格式中，向前欧拉格式是显格式，而后两者都是隐格式。事实上，这三种格式都可以认为是广义上的加权隐格式：当权重因子θ分别取θ＝0、θ＝1和θ＝0.5时，即对应为向前欧拉格式、向后欧拉格式和Crank-Nicolson格式(向前欧拉格式由于权重因子为0，因此是显格式)。文献1(L.Lapidus,G.F.Pinder,NumericalSolution of Partial Differential Equations in Science and Engineering.JohnWiley,1982.)通过一个极简单的例子从理论上论证了，通过选取一个适当的权重因子值θ，可以使数值解的截断误差达到指定的量级，从而使得结果具有较高的精确度。文献2(冯新龙,王焕焕,冯新龙,阿布都热西提,求解一维热传导方程数值解的高精度方法.新疆大学学报,第17卷第3期,2000.)和文献3(王焕焕,冯新龙,阿布都热西提,最优加权因子初步探讨.新疆大学学报,第17卷第4期,2000.)分别通过类似的简单例子指出，在给定的网格参数条件下，理论上存在一个最优权重因子值θ_opt，使得数值解的精确度能高于一般常用权重因子(即0、1和0.5)的结果。

上述现有技术的不足之处主要体现在以下几点：

一、对于向前欧拉格式：精确度不够高；结果不是无条件稳定，而是网格参数必须满足一定的条件才行。

二、对于向后欧拉格式：虽然结果无条件稳定，但是精确度低，通常不如向前欧拉格式。

三、对于Crank-Nicolson格式：虽然结果无条件稳定，通常精确度也高于向前欧拉格式和向后欧拉格式，但是与前两种格式一样，权重因子也是固定值，缺乏基于具体问题和具体网格参数的针对性，因此对于大多数问题，其精确度和理想水平相比仍有较大差距。

四、对于上述文献论述的基于最优权重因子的加权隐格式：这些文献实际上仅仅是通过简单的例子指出了最优权重因子值在理论上的存在性，但并未给出如何在实际问题中搜寻出其值；另外，文献中所举的例子较为简单，仅限于没有外源项的情形，而在激光烧蚀过程中，激光脉冲持续加热靶材物质，热传导方程有外部热源项，远比文献所述情形复杂。因此，上述现有文献中论述的内容无法满足求解实际激光烧蚀热传导方程的应用需求。

发明内容

针对现有技术的上述不足，本发明提供了一种基于最优权重因子的激光烧蚀高精度数值模拟方法，适用于激光与物质相互作用的数值模拟技术领域。本方法基于有限差分法来求解一维热传导方程，将其应用于数值模拟激光脉冲烧蚀靶材的温度场分布演化。本发明的核心工作是设计了一套算法，该算法可以准确高效地搜寻出针对当前具体问题和具体网格参数的有限差分最优权重因子值。本发明的技术方案如下：

本技术方案的总体流程可分为五大步骤，如说明书附图1所示。具体阐述如下：步骤1.预备工作：

(1)确定激光参数，包括激光脉冲的峰值光强I₀、脉冲宽度τ、脉冲中心时刻t_c；确定靶材物质的物性参数，包括导热系数κ、比热容C_p、质量密度ρ、光吸收系数a、表面反射率R。

(2)确定描述激光烧蚀靶材温度场的热传导方程，假定热传导仅沿着垂直于物质表面的方向进行，靶材温度场T(x,t)满足以下一维热传导偏微分方程

其中T表示温度，x表示距离物质表面的深度，t表示时间，I(x,t)表示靶材物质受到的激光光强，其表达式为

其中关于时间t的e指数项表示激光脉冲的时间包络为高斯型包络。

(3)将热传导方程表示成一种标准形式，其表达式为

其中

β为热传导方程二阶项系数，f(x,t)为外部热源项，并满足β>0，f(x,t)>0。步骤2.确定数值计算的离散化时空网格：

(1)确定计算网格的空间范围[x_min,x_max]，其中x_min为x的最小值，x_min＝0表示物质的表面位置，x_max为x的最大值，表示距离物质表面的最深的位置，确定网格的空间步长dx，表示空间维度上相邻网格格点的间距。

(2)确定计算网格的时间范围[t_min,t_max]，其中t_min为t的最小值，t_min＝0表示温度场演化的起始时刻，t_max为t的最大值，表示结束时刻，确定网格的时间步长dt，表示时间维度上相邻网格格点的间距。

(3)确定计算网格的总行数和总列数，总行数表示空间维度的总格点数，为Nx+1，总列数表示时间维度的总格点数，为Nt+1，其中Nx＝(x_max-x_min)/dx，Nt＝(t_max-t_min)/dt；用一个(Nx+1)×(Nt+1)大小的矩阵来表示整个温度场T(x,t)，将第j行第k列的矩阵元记为(x_j,t_k)，其中1≤j≤Nx+1，1≤k≤Nt+1。

(4)确定计算网格的步长比参数r，其定义式为

并满足r>0。

步骤3.对原激光烧蚀实际问题进行特定改造，构建一个具有解析解的定解问题：

(1)对于f(x,t)中关于时间t的e指数项f(t)

选定一个截止时间t_z，计算出[0,t_z]时间范围的f(t)并对其进行指数函数拟合，拟合函数F(t)的形式为

F(t)＝b₁·exp(b₂t) (8)

得到拟合系数b₁和b₂，其中b₁>0，b₂>0。

(2)构建以下定解问题，该定解问题包括偏微分方程、初值条件和边界条件，其中偏微分方程为

其中

为定解问题的外源项系数，并满足f_s>0；

初值条件为

T(x,0)＝u_s·exp(-ax) (10)

边界条件为

对于由方程(9)、方程(10)和方程(11)构成的定解问题，存在一个精确的解析解T_ex(x,t)，其形式为

T_ex(x,t)＝u_s·exp(-ax)·exp(b₂t) (12)

其中解析解的系数u_s＝f_s/(b₂-βa²)。

步骤4.结合线性拟合迭代与非线性拟合，搜寻针对该定解问题的有限差分最优权重因子：

(1)基于步骤2第3点中已确定的时间格点参数Nt，选定M个试验时间格点参数Nt₁,Nt₂,……Nt_M，其中M≥2，并满足Nt_m<Nt，其中m＝1,2,……M；然后根据每一个Nt_m值，计算出对应的时间步长dt_m，以及步长比参数r_m，并确保max{r₁,r₂,……r_M}≤1/2。

(2)针对每一个Nt_m值，根据有限差分加权隐格式的标准求解形式，求解步骤3第2点中所述的定解问题，该标准求解形式为

其中为θ的加权隐格式的权重因子，0≤θ≤1；

表示网格点(x_j,t_k)的温度值T(x_j,t_k)；

表示外源项f_s·exp(-ax)·exp(b₂t)在虚拟网格点(x_j,t_k+θ)的函数值，该值是f_s·exp(-ax)·exp(b₂t)在网格点(x_j,t_k+1)和(x_j,t_k)的函数值的加权组合，可表示为

(3)分别取θ＝1和θ＝0，按照通用的有限差分方法，求解上述定解问题，分别求得温度场的数值解T_θ(x,t)，即T₁(x,t)和T₀(x,t)。

(4)对于每一个θ值求得的数值解T_θ(x,t)，求出其在各个网格点上的相对误差值ε_jk,θ

并求出其平均值<ε_jk>_θ：

相对误差平均值的绝对值|<ε_jk>_θ|作为评价数值解T_θ(x,t)精确度的标准，|<ε_jk>_θ|值越小，表明结果的精确度越高；

(5)对于θ＝1和θ＝0，分别求出<ε_jk>₁和<ε_jk>₀；如果<ε_jk>₁·<ε_jk>₀>0，则说明当前定解问题的最优权重因子θ_opt不在(0,1)区间内，此时直接比较|<ε_jk>₁|和|<ε_jk>₀|的大小，若|<ε_jk>₁|<|<ε_jk>₀|则判定θ_opt＝1，若|<ε_jk>₁|>|<ε_jk>₀|则判定θ_opt＝0，搜寻完成；如果<ε_jk>₁·<ε_jk>₀<0，则说明当前定解问题的最优权重因子θ_opt在(0,1)区间内，继续进行下一步搜寻。

(6)对于<ε_jk>₁·<ε_jk>₀<0的情形，二者必然一正一负，将正项记为<ε>₊，将负项记为<ε>_-，将正项对应的θ值记为θ₊，将负项对应的θ值记为θ_-。

(7)设定最优权重因子值θ_opt的准确度要求，以可容忍误差w来表示。

(8)取最优权重因子初始试探值

(9)取权重因子

求出数值解T_θ(x,t)，以及相应的相对误差平均值<ε_jk>_θ，若<ε_jk>_θ>0，则将<ε_jk>_θ赋值给<ε>₊，将θ值赋值给θ₊，若<ε_jk>_θ<0，则将<ε_jk>_θ赋值给<ε>_-，将θ值赋值给θ_-。

(10)对(θ₊,<ε>₊)和(θ_-,<ε>_-)这两组数据进行线性拟合，拟合函数形式为

<ε>＝p₁·θ+p₂ (16)

求得拟合系数p₁和p₂，并计算更新一步的最优权重因子值

其计算式为

(11)判断

是否成立：如果成立，则判定

搜寻完成，如果不成立，则令

(12)重复进行步骤4第9点至第11点的操作，直至

成立，迭代搜寻完成，得到最优权重因子值θ_opt。

(13)对于每一个Nt_m值及其对应的步长比参数r_m，求得一个相应的最优权重因子值θ_opt,m，总共得到M组(r_m,θ_opt,m)数据后，对其进行非线性拟合，拟合函数形式为

求得拟合系数p₃和p₄(通常p₃>0，p₄<0)。

(14)对于方程(6)中已确定的网格步长比参数r，由方程(17)计算出最优权重因子值θ_opt。若θ_opt≥0，则这个值就作为步骤3所述定解问题的最终最优权重因子值；若θ_opt<0，则说明步长比参数r过小，应适当修改网格参数以获得一个更大的r值，然后重新计算θ_opt，直至满足θ_opt≥0。

步骤5.基于最优权重因子求解原问题：

(1)确定原激光烧蚀靶材问题的初值条件和边界条件。

(2)基于步骤4求得的最优权重因子值θ_opt，用有限差分法求解原问题。

本发明的作用原理如下：

首先，对于有限差分加权隐格式来说，假定构建的定解问题的精确解T_ex(x,t)取值恒为正数，那么在权重因子取值θ＝1时，数值解T_θ(x,t)是位于精确解T_ex(x,t)的正方向，即对于任意网格点(x_j,t_k)，均有相对误差值ε_jk,θ>0，因而平均误差值<ε_jk>_θ>0，而在权重因子取值θ＝0时，数值解T_θ(x,t)是位于精确解T_ex(x,t)的负方向，即对于任意网格点(x_j,t_k)，均有相对误差值ε_jk,θ<0，因而平均误差值<ε_jk>_θ<0。当θ的取值从1逐渐变到0时，<ε_jk>_θ是从正值逐渐变为负值，因此，理论上必然存在一个介于0到1之间的权重因子值θ，使得<ε_jk>_θ＝0，即数值解完全等于精确解，该θ值即为理论上的最优权重因子值θ_opt。虽然在网格参数非常特殊的极端条件下，可能会出现θ_opt不在(0,1)区间内的情况，但这种情况不在本发明的讨论范围以内，已在相关步骤中予以排除。

其次，要搜寻出最优权重因子值θ_opt，需要有一个可供参考的精确解析解，而实际的激光烧蚀问题由于偏微分方程、初值条件和边界条件的复杂性，通常并无解析解，因此，本发明通过构建一个具有精确解析解的定解问题，可以获得一个精确解供搜寻时使用。

第三，对于一个具体的定解问题而言，最优权重因子θ_opt值和偏微分方程参数以及数值计算的网格参数有关，而与初值条件、边界条件无关。在本发明构建的定解问题中，偏微分方程与原问题的方程是高度近似的(内部热传导部分完全一致，外部热源部分在截取的脉冲前沿时间范围[0,t_z]内高度近似)，最后求解最优权重因子时所用的网格参数与原问题的网格参数是完全一致的，因此定解问题的最优权重因子可以用于求解原问题。

第四，在固定的Nx参数下，Nt越大，步长比参数r越小，意味着时间维度划分越细，结果也越精确，而计算耗时也越多。由于最优权重因子θ_opt和步长比参数r之间满足方程(17)所示的非线性关系，因此本发明采用的方法是先在几个较大的r值下通过线性拟合迭代来搜寻各自的最优权重因子，而后通过非线性拟合来求出在较小的r值下的最优权重因子，和直接在较小的r值下迭代搜寻最优权重因子相比，本发明的方法可以有效地节省计算耗时，使整个算法更高效。

最后，由迭代算法搜寻得到的最优权重因子值θ_opt有可能小于1/2，如果确实θ_opt<1/2，那么取θ＝θ_opt进行有限差分求解时，其结果将不是无条件稳定的，必须要满足r≤1/2(1-2θ)才能确保结果是稳定的。在步骤4第1点中，已经限制了max{r₁,r₂,……r_M}≤1/2，而最后的实际目标网格步长比参数r又比所有的r_m更小，因此一定满足r<1/2，又因为0<θ<1/2，因此1/2(1-2θ)一定是大于1/2的值，所以稳定性条件r≤1/2(1-2θ)必然得到满足。由此可知，基于本发明求得的θ_opt值在相应的r参数条件下来进行有限差分求解，可以确保其结果是稳定的。

有益效果

与现有技术相比，一种基于最优权重因子的激光烧蚀高精度数值模拟方法具有以下优点：

一、与常规的有限差分加权隐格式方法相比，本方法专门根据当前具体问题和具体网格参数来搜寻最优权重因子，更有针对性，因此，基于最优权重因子来求解热传导方程，其数值解的精确度能够显著高于基于一般常用权重因子的结果(相对误差平均值通常可以小若干个数量级)，并且可以确保数值结果是稳定的。

二、与现有相关文献论述的基于最优权重因子的加权隐格式相比，本发明不是停留在讨论最优权重因子理论存在性的层面，而是给出了切实可行的在实际问题中搜寻最优权重因子值的方法。另外，本发明的方法不仅适用于热传导方程不含外源项的情形，还适用于热传导方程包含外源项的情形，因此能够满足激光烧蚀数值模拟的实际应用需求。

三、本发明搜寻最优权重因子的方法结合了线性拟合迭代与非线性拟合，该方法能够简捷高效地搜寻出针对当前问题和网格参数的最优权重因子，并且求得的最优权重因子能够具有很高的准确度。

综上，本发明能够基于有限差分最优权重因子对激光烧蚀进行高精度数值模拟，并具有搜寻最优权重因子方法准确高效，基于最优权重因子数值求解结果精确稳定的优点，本发明适用于激光脉冲烧蚀靶材过程的数值模拟，尤其适用于着重关注脉冲前沿对靶材加热过程的数值模拟。

附图说明

图1为技术方案总体流程示意图。

图2为线性拟合迭代搜寻最优权重因子方法流程示意图。

图3为具体实施案例中非线性拟合搜寻最优权重因子效果示意图。

具体实施方式

下面结合一个激光烧蚀数值模拟具体实例，阐明在发明内容部分所述方法的应用过程：

1.在本实例中，物理参数均使用国际单位制。设定激光参数如下：峰值光强I₀＝1.0×10¹³W/m²，脉冲宽度τ＝4ns＝4×10^-9s，脉冲中心时刻t_c＝6ns＝6×10^-9s；设定靶材物性参数如下：导热系数k＝14.63W/(m·K)，热容C_p＝540J/(kg·K)，质量密度ρ＝4510kg/m³，吸收系数a＝1.0×10⁷/m，表面反射率R＝0.6。由于下面将开始数值计算，为简明起见，表达式一般不再包含物理单位。

由方程(1)和方程(2)可知，靶材温度场T(x,t)满足偏微分方程

其中

简单计算可知，热传导方程可以表示成如下的一种标准形式：

其中

由以上表达式易知满足β>0，f(x,t)>0。

2.确定计算网格的空间范围[x_min,x_max]，设定考察的起始位置为物质表面，终止位置为距离表面10.0μm处，即x_min＝0，x_max＝10.0×10^-6，设定空间步长为10.0nm，即dx＝10.0×10^-9。确定计算网格的时间范围[t_min,t_max]，设定考察的起始时刻为时间零点，结束时刻为15.0ns，即t_min＝0，t_max＝15.0×10^-9，设定时间步长为0.00075ns，即dt＝7.5×10^-13。易知Nx＝1000，Nt＝20000，整个温度场T(x,t)由一个1001×20001的矩阵来表示。网格的步长比参数r为

满足r>0。

3.对于f(x,t)中关于时间t的e指数项f(t)

选定截止时间为0.45ns，即t_z＝4.5×10^-10，对[0,4.5×10^-10]时间范围的f(t)按方程(8)所示的指数函数形式进行拟合，拟合结果为b₁＝0.002041，b₂＝1.8876×10⁹，即F(t)＝0.002041·exp(1.8876×10⁹t)，满足b₁>0，b₂>0。

构建以下定解问题，该定解问题的偏微分方程为

其中

满足f_s>0。

初值条件为

T(x,0)＝u_s·exp(-1.0×10⁷x)

边界条件为

其中

u_s＝3.34724×10¹⁰/[1.8876×10⁹-6.0×10^-6×(1.0×10⁷)²]＝25.99596对于这个定解问题，存在一个精确的解析解T_ex(x,t)

T_ex(x,t)＝25.99596·exp(-1.0×10⁷x)·exp(1.8876×10⁹t)

易知，对于任意x和t，均有T_ex(x,t)>0。

4.基于已确定的时间格点参数Nt＝20000，选定M个试验时间格点参数，这里取M＝6，并设定Nt₁＝2000，Nt₂＝2500，Nt₃＝3000，Nt₄＝3600，Nt₅＝4000，Nt₆＝4500；其对应的时间步长分别为dt₁＝0.7500×10^-11，dt₂＝0.6000×10^-11，dt₃＝0.5000×10^-11，dt₄＝0.4166×10^-11，dt₅＝0.3750×10^-11，dt₆＝0.3333×10^-11；其对应的步长比参数分别为r₁＝0.45，r₂＝0.36，r₃＝0.30，r₄＝0.25，r₅＝0.225，r₆＝0.20，易知max{r₁,r₂,……r₆}＝r₁，满足r₁≤1/2。

以Nt₁＝2000的情况为例，有限差分加权隐格式的标准求解形式为

取θ＝1，按照通用的有限差分求解方法，求得数值解T₁(x,t)，由方程(15)可计算出<ε_jk>₁＝0.01009305，类似地取θ＝0，求得数值解T₀(x,t)，可计算出<ε_jk>₀＝-0.00926238。显然<ε_jk>₁·<ε_jk>₀<0，说明当前定解问题的最优权重因子θ_opt在(0,1)区间内，继续进行下一步搜寻。将<ε_jk>₁记为<ε>₊，将<ε_jk>₀记为<ε>_-，令θ₊＝1，θ_-＝0。设定最优权重因子值θ_opt的准确度要求，设定可容忍误差w＝0.00001。取最优权重因子初始试探值

取权重因子θ＝0.5，求出数值解T_0.5(x,t)，并可计算出<ε_jk>_0.5＝0.00038501。由于<ε_jk>_0.5>0，因此相关值更新为<ε>₊＝<ε_jk>_0.5，以及θ₊＝0.5。接下去，对(0.5,0.00038501)和(0,-0.00926238)这两组数据按方程(16)进行线性拟合，得到拟合系数p₁＝0.01929479，p₂＝-0.00926238，并求出

显然

因此需要进行下一步迭代计算。将

值更新为

取权重因子θ＝0.48004，求出数值解T_0.48004(x,t)，并可计算出<ε_jk>_0.48004＝-1.2691516×10^-6。由于<ε_jk>_0.48004<0，因此相关值更新为<ε>_-＝<ε_jk>_0.48004，以及θ_-＝0.48004。接下去，对(0.5,0.00038501)和(0.48004,-1.2691516×10^-6)这两组数据进行线性拟合，得到拟合系数p₁＝0.01935282，p₂＝-0.00929140，并求出

显然

因此仍需要进行下一步迭代计算。将

值更新为

取权重因子θ＝0.48010，求出数值解T_0.48010(x,t)，并可计算出<ε_jk>_0.48010＝-1.0812718×10^-7。类似地，相关值更新为<ε>-＝<ε_jk>_0.48010，以及θ-＝0.48010，并对(0.5,0.00038501)和(0.48010,-1.0812718×10^-7)这两组数据进行线性拟合，可以求出

这次，

满足

因此线性拟合迭代搜寻完成，得到Nt₁＝2000情况下的最优权重因子值θ_opt,1＝0.48010，也即(r₁,θ_opt,1)＝(0.45,0.48010)。

完全类似地，分别求出其它5个Nt参数下的最优权重因子值，结果分别为：θ_opt,2＝0.47566，θ_opt,3＝0.47114，θ_opt,4＝0.46566，θ_opt,5＝0.46198，θ_opt,6＝0.45737。对于这6组(r_m,θ_opt,m)已知点数据，即(0.45,0.48010)，(0.36,0.47566)，(0.30,0.47114)，(0.25,0.46566)，(0.225,0.46198)，(0.20,0.45737)，按照方程(17)的形式进行非线性拟合，得到p₃＝0.49838733，p₄＝-0.00819219。得到拟合曲线后，即可预测实际目标步长比参数r对应的θ_opt值。对于目标网格步长比参数r＝0.045，其最优权重因子值为θ_opt＝0.31634，由于满足θ_opt>0，因此可以直接认定为当前定解问题的最终最优权重因子值。

为展现最优权重因子的优点，对于当前定解问题，在Nt＝20000参数条件下，基于各个不同权重因子值θ进行有限差分数值求解，所得结果的相对误差平均值<ε_jk>θ如表1所示。由表1可以看出，基于最优权重因子的计算结果的误差要比一般常用权重因子的结果小得多，即使是通常认为“很好”的0.5，其误差还是要比0.31634的误差大2个数量级。由此可见，基于最优权重因子计算结果的精确度能够显著高于基于一般常用权重因子的结果。

权重因子值θ	1	0	0.5	0.31634
					平均相对误差&lt;ε<sub>jk</sub>&gt;<sub>θ</sub>	1.33×10<sup>-3</sup>	-6.05×10<sup>-4</sup>	3.62×10<sup>-4</sup>	6.98×10<sup>-6</sup>

表1

另外，如果直接用本发明的线性拟合迭代方法来精确搜寻Nt＝20000参数条件下的最优权重因子，可得其值为

对比θ_opt和

可知，通过本发明的非线性拟合预测得到的结果和直接精确搜寻得到的结果之间的相对误差仅为1％左右，因此具有很高的准确度。而相比于直接精确搜寻，用非线性拟合的方法仅需在若干个较小的Nt参数下进行计算，而不用直接在Nt＝20000这样的大参数下进行计算，可以有效地节省计算耗时，使整个方法更为简捷高效。

5.获得最优权重因子值θ_opt之后，只需要根据实际激光烧蚀靶材问题的初值条件和边界条件，将θ_opt作为权重因子，用有限差分法求解该问题即可，在此不作赘述。

备注：

1.以上具体实施方案体现了本发明的基本原理、主要特征和优点。相关技术人员应该了解，本发明不受上述具体实施案例的限制。在不脱离本发明精神和范围的前提下，本发明可能有各种变化和改进。这些变化和改进都落入要求保护的本发明的范围内。本发明要求的保护范围由所附的权利要求书及其等同物界定。

Claims (1)

1.一种基于最优权重因子的激光烧蚀高精度数值模拟方法，其特征在于包含以下步骤：

1)预备工作；

1-1)确定激光参数，包括激光脉冲的峰值光强I₀、脉冲宽度τ、脉冲中心时刻t_c；确定靶材物质的物性参数，包括导热系数κ、比热容C_p、质量密度ρ、光吸收系数a、表面反射率R；

1-2)确定描述激光烧蚀靶材温度场的热传导方程，假定热传导仅沿着垂直于物质表面的方向进行，靶材温度场T(x，t)满足以下一维热传导偏微分方程

其中T表示温度，x表示距离物质表面的深度，t表示时间，I(x，t)表示靶材物质受到的激光光强，其表达式为

其中关于时间t的e指数项表示激光脉冲的时间包络为高斯型包络；

1-3)将热传导方程表示成一种标准形式，其表达式为

其中

β为热传导方程二阶项系数，f(x，t)为外部热源项，并满足β＞0，f(x，t)＞0；

2)确定数值计算的离散化时空网格；

2-1)确定计算网格的空间范围[x_min，x_max]，其中x_min为x的最小值，x_min＝0表示物质的表面位置，x_max为x的最大值，表示距离物质表面的最深的位置，确定网格的空间步长dx，表示空间维度上相邻网格格点的间距；

2-2)确定计算网格的时间范围[t_min，t_max]，其中t_min为t的最小值，t_min＝0表示温度场演化的起始时刻，t_max为t的最大值，表示结束时刻，确定网格的时间步长dt，表示时间维度上相邻网格格点的间距；

2-3)确定计算网格的总行数和总列数，总行数表示空间维度的总格点数，为Nx+1，总列数表示时间维度的总格点数，为Nt+1，其中Nx＝(x_max-x_min)/dx，Nt＝(t_max-t_min)/dt；用一个(Nx+1)×(Nt+1)大小的矩阵来表示整个温度场T(x，t)，将第j行第k列的矩阵元记为(x_j，t_k)，其中1≤j≤Nx+1，1≤k≤Nt+1；

2-4)确定计算网格的步长比参数r，其定义式为：

并满足r＞0；

3)对原激光烧蚀实际问题进行特定改造，构建一个具有解析解的定解问题；

3-1)对于f(x，t)中关于时间t的e指数项f(t)

选定一个截止时间t_z，计算出[0，t_z]时间范围的f(t)并对其进行指数函数拟合，拟合函数F(t)的形式为：

F(t)＝b₁·exp(b₂t)

得到拟合系数b₁和b₂，其中b₁＞0，b₂＞0；

3-2)构建以下定解问题，该定解问题包括偏微分方程、初值条件和边界条件，其中偏微分方程为：

其中

为定解问题的外源项系数，并满足f_s＞0；

初值条件为：

T(x，0)＝u_s·exp(-ax)

边界条件为

对于由以上偏微分方程、初值条件和边界条件构成的定解问题，存在一个精确的解析解T_ex(x，t)，其形式为：

T_ex(x，t)＝u_s·exp(-ax)·exp(b₂t)

其中解析解的系数u_s＝f_s/(b₂-βa²)；

4)结合线性拟合迭代与非线性拟合，搜寻针对该定解问题的有限差分最优权重因子；

4-1)基于步骤2-3中已确定的时间格点参数Nt，选定M个试验时间格点参数Nt₁，Nt₂，……Nt_M，其中M≥2，并满足Nt_m＜Nt，其中m＝1，2，……M；然后根据每一个Nt_m值，计算出对应的时间步长dt_m，以及步长比参数r_m，并确保max{r₁，r₂，……r_M}≤1/2；

4-2)针对每一个Nt_m值，根据有限差分加权隐格式的标准求解形式，求解步骤3-2中的定解问题，该标准求解形式为：

其中θ为加权隐格式的权重因子，0≤θ≤1；

表示网格点(x_j，t_k)的温度值T(x_j，t_k)；

表示外源项f_s·exp(-ax)·exp(b₂t)在虚拟网格点(x_j，t_k+θ)的函数值，该值是f_s·exp(-ax)·exp(b₂t)在网格点(x_j，t_k+1)和(x_j，t_k)的函数值的加权组合，可表示为

4-3)分别取θ＝1和θ＝0，按照通用的有限差分方法，求解上述定解问题，分别求得温度场的数值解T_θ(x，t)，即T₁(x，t)和T₀(x，t)；

4-4)对于每一个θ值求得的数值解T_θ(x，t)，求出其在各个网格点上的相对误差值ε_jk，θ：

并求出其平均值<ε_jk>_θ：

相对误差平均值的绝对值|<ε_jk>_θ|作为评价数值解T_θ(x，t)精确度的标准，|<ε_jk>θ|值越小，表明结果的精确度越高；

4-5)对于θ＝1和θ＝0，分别求出<ε_jk>₁和<ε_jk>₀；如果<ε_jk>₁·<ε_jk>₀＞0，则说明当前定解问题的最优权重因子θ_opt不在(0，1)区间内，此时直接比较|<ε_jk>₁|和|<ε_jk>₀|的大小，若|<ε_jk>₁|＜|<ε_jk>₀|则判定θ_opt＝1，若|<ε_jk>₁|＞|<ε_jk>₀|则判定θ_opt＝0，搜寻完成；如果<ε_jk>₁·<ε_jk>₀＜0，则说明当前定解问题的最优权重因子θ_opt在(0，1)区间内，继续进行下一步搜寻；

4-6)对于<ε_jk>₁·<ε_jk>₀＜0的情形，二者必然一正一负，将正项记为<ε>₊，将负项记为<ε>_-，将正项对应的θ值记为θ₊，将负项对应的θ值记为θ_-；

4-7)设定最优权重因子值θ_opt的准确度要求，以可容忍误差w来表示；

4-8)取最优权重因子初始试探值

4-9)取权重因子

求出数值解T_θ(x，t)，以及相应的相对误差平均值<ε_jk>_θ，若<ε_jk>_θ＞0，则将<ε_jk>_θ赋值给<ε>₊，将θ值赋值给θ₊，若<ε_jk>_θ＜0，则将<ε_jk>_θ赋值给<ε>_-，将θ值赋值给θ_-；

4-10)对(θ₊，<ε>₊)和(θ_-，<ε>_-)这两组数据进行线性拟合，拟合函数形式为：

<ε>＝p₁·θ+p₂

求得拟合系数p₁和p₂，并计算更新一步的最优权重因子值

其计算式为

4-11)判断

是否成立：如果成立，则判定

搜寻完成，如果不成立，则令

4-12)重复进行步骤4-9至4-11的操作，直至

成立，线性拟合迭代搜寻完成，得到最优权重因子值θ_opt；

4-13)对于每一个Nt_m值及其对应的步长比参数r_m，求得一个相应的最优权重因子值θ_opt，m，总共得到M组(r_m，θ_opt，m)数据后，对其进行非线性拟合，拟合函数形式为：

求得拟合系数p₃和p₄；

4-14)对于步骤2-4中已确定的网格步长比参数r，由步骤4-13中求得的拟合公式，计算出最优权重因子值θ_opt；若θ_opt≥0，则这个值就作为步骤3所述定解问题的最终最优权重因子值；若θ_opt＜0，则说明步长比参数r过小，应适当修改网格参数以获得一个更大的r值，然后重新计算θ_opt，直至满足θ_opt≥0；

5)基于最优权重因子求解原问题；

5-1)确定原激光烧蚀靶材问题的初值条件和边界条件；

5-2)基于步骤4)求得的最优权重因子值θ_opt，用有限差分法求解原问题。

Images