热辐射交换系数计算方法

1.本发明涉及传热学和计算力学设计领域，具体是一种基于nurbs和蒙特卡洛射线追踪的热辐射交换系数计算方法。


背景技术：

2.有限元法(fem)和有限差分法(fdm)是两种重要的数值分析方法。一般来说，fem因其计算准确性和通用性而被广泛使用。
3.但是对于传热分析，特别是在同时考虑热传导和热辐射的航空航天领域，fdm是一种主要方法。
4.在进行热辐射分析时，有限元法采用有限元网格将几何体离散为众多的单元，计算大量的单元之间的热辐射交换系数(ref)时非常耗时。相比而言，fdm处理辐射传热时的计算工作量更小。但是其相比fem而言计算精度有待提高。
5.高精度的热辐射交换系数的计算需要准确且连续变化的表面法向量。fem虽然计算精度较高，但是当考虑镜面反射时，有限元网格生成过程产生的小平面使得表面法线方向不连续，可能会影响精度。
6.对于上述问题，相关研究提出以下几种改进方案：1.p.vueghs等人提出了一种将原始几何描述与有限元网格相结合的方法。(vueghs,p.,de koning,h.p.,pin,o.,&beckers,p.(2010).use of geometry in finite element thermal radiation combined with ray tracing.journal of computational and applied mathematics,234(7),2319-2326.)该方法使用球体、圆柱体来描述几何形状并进行蒙特卡洛射线追踪，使用有限元网格进行热分析。该方法存在一定局限性：需要将有限元网格投影到原始几何体上，以确定射线的来源单元；该方法几何形状描述能力有限，对于更高阶表面，通常没有解析解来计算射线与表面单元的交点。2.t.walker等人基于上述方法通过引入仿射变换可以描述更多的形状。(walker t,xue s c,barton g w.numerical determination of radiative view factors using ray tracing[j].journal of heat transfer,2010,132(7).)。3.l.jacques等人提出将有限元网格的顶点拟合到二次曲面中，后将单元上生成的采样点投影到二次曲面上，从而用近似曲面法向量来替代精确的法向量。(jacques,l.,bechet,e.,&kerschen,g.(2017).finite element model reduction for space thermal analysis.finite elements in analysis and design,127,6-15.)
[0007]
上述方法提出了通过fem提高热辐射交换系数的计算方法，但3种方法中，无论是使用一致的有限元网格来进行射线追踪和热分析，还是使用两套不同的网格来执行这两个过程，其计算结果都不是理想的。
[0008]
此外，fem方法由于必须将几何体离散成单元进行分析，不可避免地导致离散后的几何表面不光滑连续或难以做到高阶连续。而精准的热辐射交换系数的计算需要准确且连续变化的表面法向量。上述3种基于有限元离散网格的分析方法从本质上并没有解决这一问题。


技术实现要素：

[0009]
本发明为了解决现有技术的问题，提供了一种热辐射交换系数计算方法，通过等几何方法(iga)使得用于计算ref的蒙特卡洛射线追踪和热分析都可以在同一个nurbs网格上进行，保证了几何体表面几何精确、高阶连续等特性，提升了ref计算效率和结果精度，具有较强的通用性，因而易于在热辐射和热传导等分析中推广使用。
[0010]
iga是一种新型的有限元方法，其将cad中描述几何形状的非均匀有理b样条(nurbs)引入到等参有限元中，消除了计算机辅助几何设计(cad)与计算机辅助分析(cae)之间反复的数据转换过程，节省了大量的前处理时间，而且其具有几何精确、高阶连续等特性。
[0011]
为实现以上目的，本发明包括以下技术步骤：
[0012]
步骤1：根据研究的热辐射问题建立几何体的nurbs模型及其控制网格；
[0013]
步骤2：基于完全漫反射表面假设和琅勃特余弦定律，在表面半球空间内对射线方向进行随机均匀采样；
[0014]
步骤3：假设几何体表面光学特性处处相同，采用接受-拒绝采样方法在nurbs模型表面进行均匀采样；
[0015]
步骤4：利用贝塞尔曲面的凸包特性，通过贝塞尔剪裁法获取射线和nurbs模型曲面的交叉点；
[0016]
步骤5：在步骤2、3、4的采样完成的基础上，使用蒙特卡洛数值积分将热辐射交换系数(ref)的计算转化为样本期望值的计算，从而获得该热辐射问题的ref。至此整个分析过程结束。
[0017]
以上所述步骤中：
[0018]
步骤1中在参数空间中定义nurbs曲面，并可以通过改变控制点获得几何精确且高阶连续的曲面。步骤2-5中的蒙特卡洛射线追踪和热分析过程均可以在步骤1建立的同一套nurbs网格上进行；
[0019]
步骤2中基于完全漫反射表面假设和琅勃特余弦定律，射线在表面半球空间内均匀分布，采样的具体公式如下：
[0020][0021][0022]
式中：θ和分别为球坐标下的极角和方位角；r
θ
和r
φ
为[0,1]内的均匀采样随机数。
[0023]
步骤3中使用了接受-拒绝采样方法，该采样的具体方式为：
[0024]
步骤3.1：在参数空间中进行随机采样生成一组随机数(ru，rv)和随机点(u,v)；
[0025]
步骤3.2：剔除不满足约束条件的随机点。约束条件具体为：
[0026][0027][0028]
式中：为nurbs曲面的雅可比矩阵的行列式；
[0029]
对于连续的nurbs曲面可写成如下参数形式：
[0030][0031]
式中x为连续曲面的参数形式，(u,v)是参数空间中的坐标，ω是参数域，x,y,z是(u,v)的可微函数。则e、f、g具体表达式为：
[0032][0033][0034][0035]
步骤4中，由于没有针对射线和nurbs曲面交叉点的解析解，此步骤采用数值解法—贝塞尔剪裁法。其利用贝塞尔曲面的凸包特性，无需初始猜测便可以获得射线和nurbs曲面交叉点。
[0036]
步骤5中基于步骤2-4的采样方案，计算ref的表达式如下：
[0037][0038][0039]
式中：dai和daj是两几何体表面的面积微元；s
ij
是两微元之间连线的长度；θi和θj分别是连线与表面积微元法向量ni和nj的夹角；dω
θ
为微分立体角；χj(x,θ)为预测值仅可取0或1。
[0040]
采用概率密度函数的蒙特卡洛数值积分，将ref计算转化为样本的期望值计算。具体公式为：
[0041][0042][0043][0044]
至此整个分析过程结束。
[0045]
本发明有益效果在于：
[0046]
通过等几何方法(iga)将非均匀有理b样条(nurbs)引入到等参有限元中，对于热辐射问题的几何体建立nurbs模型及其控制网格，使得用于计算ref的蒙特卡洛射线追踪和
热分析两个过程都可以在同一个nurbs网格上进行，节省了大量前处理时间，提高了ref计算效率。同时nurbs模型保证了几何体表面几何精确、高阶连续等特性，提高了ref的结果精度。相比现有计算ref的fdm和fem方法，本发明具有较强的通用性，易于在热辐射和热传导等分析中推广使用。
附图说明
[0047]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
[0048]
图1为球面-圆柱面热辐射模型示意图；
[0049]
图2为球体nurbs模型及接受-拒绝采样结果示意图；
[0050]
图3为球面-圆柱面的热辐射交换系数曲线图；
[0051]
图4为热辐射交换系数ref相对误差-对数采样射线数关系图。
具体实施方式
[0052]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
[0053]
本发明的具体步骤如下：
[0054]
步骤1：本发明所示实施例为球体表面对圆柱体表面的热辐射模型(如图1所示)，其中r和l分别为圆柱体的半径和高度；s是球心与圆柱体底部之间的距离；圆柱体的中轴线穿过球心。
[0055]
如图2所示，在参数空间[0,1]
×
[0,1]中定义nurbs曲面，并通过改变控制点获得几何精确且高阶连续的曲面，从而建立球面-圆柱面的nurbs模型及其控制网格，使得后续的蒙特卡洛射线追踪和热分析过程均可以在此步骤建立的同一套nurbs网格上进行。
[0056]
步骤2：基于完全漫反射表面假设和琅勃特余弦定律，射线在球-圆柱表面半球空间内均匀分布，对射线方向进行随机均匀采样，采样的具体公式如下：
[0057][0058][0059]
式中：θ和分别为球坐标下的极角和方位角；r
θ
和r
φ
为[0,1]内的均匀采样随机数。
[0060]
步骤3：假设几何体表面光学特性处处相同，采用接受-拒绝采样方法在nurbs模型表面进行均匀采样，具体方式为：
[0061]
(1)在参数空间[0,1]
×
[0,1]中进行随机采样生成一组随机数(ru，rv)和随机点(u,v)；
[0062]
(2)剔除不满足约束条件的随机点。约束条件具体为：
[0063][0064][0065]
式中：为球的nurbs曲面雅可比矩阵的行列式；
[0066][0067]
式中(u,v)是参数空间中的坐标，x,y,z是(u,v)的可微函数。则e、f、g具体表达式为：
[0068][0069][0070][0071]
则
[0072][0073]
如图2所示，对球体的nurbs模型表面进行采样。执行(1)后所得随机采样点并未均匀分布在表面，而在球的两极点较为集中。执行(2)后对部分采样点进行剔除进而获得均匀分布的表面采样。
[0074]
一般来说，此采样方案适用于nurbs描述的任意自由曲面的表面采样。
[0075]
步骤4：此步骤求解步骤2所得采样射线与球-圆柱nurbs曲面的交点。由于没有针对射线和nurbs曲面交点的解析解，此步骤采用数值解法—贝塞尔剪裁法。其利用贝塞尔曲面的凸包特性，无需初始猜测便可以获得步骤2中所得的采样射线和球-圆柱nurbs曲面的交点。
[0076]
步骤5：基于步骤2-4的采样方案，计算ref的表达式如下：
[0077][0078][0079]
式中：dai和daj是两几何体表面的面积微元；s
ij
是两微元之间连线的长度；θi和θj分别是连线与表面积微元法向量ni和nj的夹角；dω
θ
为微分立体角；χj(x,θ)为预测值仅可取0或1。
[0080]
对上述重积分采用数值积分求解。此处引入概率密度函数的蒙特卡洛数值积分，将ref计算转化为样本的期望值计算。具体公式为：
[0081][0082][0083][0084]
从而获得该热辐射问题的ref。至此整个计算分析过程结束。
[0085]
对于球面-圆柱面热辐射这一经典问题，其ref存在解析解，具体表达式为：
[0086][0087]
式中r,z取决于圆柱的几何外形及球体-圆柱的相对位置，r＝r/s，z＝l/r。
[0088]
如图3所示，对于不同的r,z取值，使用本发明计算方法获得的ref与相应解析解吻合较好。
[0089]
如图4所示为ref相对误差-对数采样射线数关系图。随着采样射线数的增加，相对误差快速下降。当采样射线数n
ray
达到105时，相对误差将低于0.1％。相比传统的fem和fdm，本发明实施的计算结果精度更高。
[0090]
综上所述本发明提出了一种通用的热辐射交换系数计算方法，该方法通过建立nurbs模型及其控制网格，使得用于计算ref的蒙特卡洛射线追踪和热分析可以在同一个nurbs网格上进行，同时nurbs模型保证了几何体表面几何精确、高阶连续等特性，提高了ref计算效率和结果精度。
[0091]
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于设备实施例而言，以上所述仅是本发明的优选实施方式，由于其基本相似于方法实施例，所以描述得比较简单，相关之处参见方法实施例的部分说明即可。以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，对于本技术领域的普通技术人员来说，可轻易想到的变化或替换，在不脱离本发明原理的前提下，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。