[发明专利]梯度折射率介质内瞬态辐射传递求解方法

权利要求书

1.梯度折射率介质内瞬态辐射传递求解方法，其特征在于，所述方法具体包括以下步骤：

步骤一、根据预设的计算域进行网格划分，依据每一个单元几何信息和多项式阶数生成单元节点，并将所有单元节点信息存储；具体过程为：

步骤一一、根据计算域的几何信息生成有限元网格单元，根据每一个单元几何信息和多项式阶数生成该单元节点；

步骤一二、采用结构类型数组，统一存储节点坐标及其对应单元号及局部索引信息；

步骤一三、采用切比雪夫多项式作为基函数，计算标准单元[-1,1]上的切比雪夫多项式节点坐标，计算公式为：

其中，N为每个单元中积分节点数目；

求出在单元[x_k+1,x_k]上的切比雪夫求解节点坐标：

其中，Nelem是单元数目；ξ是标准单元[-1,1]上的积分节点坐标；

步骤一四、计算任意单元上的积分节点坐标与积分权值；

利用公式(3)计算得到标准单元[-1,1]上的积分节点坐标ξ：

其中，ξ_i为高斯点坐标；H_i为相应高斯点的求积系数；

采用如下坐标变换，得到单元[x_k+1,x_k]上的积分节点坐标：

相对应的单元[x_k+1,x_k]上的积分权值为：

ω_k＝ω×yocabi (5)

其中，ω表示标准单元上的积分权值，k＝1,2,…,Nelem；

步骤一五、将每一个节点存入无重复的边界节点结构型数组，若有坐标重复则不存储，同时将重复节点的单元信息累加到已装入的节点单元集，局部索引添加到重复节点的局部索引；循环结束后，建立一个无重复边界节点结构型数组；

步骤二、根据步骤一单元信息计算标准单元上切比雪夫节点基函数及其导数在给定节点集上的函数值和权值，并生成系数矩阵；具体过程为：

依据步骤一中的公式(4)和计算出的单元上的求解节点坐标和积分节点坐标，得到标准单元上的切比雪夫节点基函数及其导数在给定节点集上的函数值和权值，并利用得到的切比雪夫节点基函数的函数值及其导数在给定节点集上的函数值和权值生成系数矩阵：

其中，eBSF(x_k)表示基函数在x_k处的矩阵值，eDBSF(x_k)表示基函数导数在x_k处的矩阵值，Wd为权值矩阵，上标“T”表示矩阵转置符号；

步骤三、利用平行脉冲照射下的梯度折射率介质内的扩散近似方程与步骤二获得的系数矩阵，计算刚度矩阵，求出辐射强度；具体过程为：

步骤三一、初始化辐射强度，设置时间步循环t＝1,…,P；

步骤三二、设介质左边界受如下激光平行脉冲q_p(t)照射：

式中，t_p为脉冲宽度，q₀为初始辐射热流；H(t)为赫维赛德阶梯函数；

步骤三三、介质内辐射分为两部分处理，即平行光分量和漫射光分量，其中平行光分量满足贝尔定律，计算式为：

式中q_c(x,t)表示在t时刻介质内部x处的平行光辐射热流；β为衰减系数，n表示折射率分布；c₀表示真空中的光速；

步骤三四、通过求解瞬态扩散近似方程获得漫射光分量，具体如下：

含激光入射的扩散近似控制方程的离散形式为：

将式(12)整理得到：

将式(13)代入式(11)得：

令

C＝1+Δt^*L_Rk_a，

D＝Δt^*L_R，

S＝Δt^*L_Rk_sG_c，

则式(14)简化为：

式中，Δt^*为无量纲时间步长；β为衰减系数；A₁为各相异性散射系数；L_R为衰减系数的倒数；κ_s为散射系数，κ_a表示吸收系数；G和G_c分别代表投射辐射和平行光投射辐射；q和q_c分别代表辐射热流和平行光辐射热流；上标“Δ”表示上一时刻场分布，为哈密顿算子；

离散后边界条件为：

其中，ε为壁面发射率，n_w为壁面法向量；

令

式(16)得到：

依据式(15)和(17)最终离散式写为矩阵形式为：KG＝H

其中：

K＝C*M+D*B/A*Mxx

H＝M*S-D*Mx^T*q_s+D/E*N*R

R＝E×q_s

漫射光分量的矩阵形式为：

其中，inv(M)表示矩阵M的逆矩阵；

步骤三五、依据时间步，循环计算，并将运行的结果统一存储；

步骤四、调用步骤三辐射强度值计算半球穿透率与半球反射率，得到半球穿透率与时间的曲线、半球反射率与时间的曲线。