CN114117877B - 一种基于等几何粒子描述的拓扑优化方法 - Google Patents

Abstract

一种基于等几何粒子描述的拓扑优化方法，采用等几何粒子流体动力学方法进行物理场仿真计算，结合基于变密度法SIMP的拓扑表征方法完成结构优化设计，根据拉格朗日观点思想和非定常流场的流场仿真的结果，完成目标函数关于设计变量的灵敏度分析，迭代优化设计变量，从而推动目标函数按照设计意图变化，促使结构向预定性能逼近，最终完成给定约束下，材料布局和结构设计的最优组合；本发明在设计之初能保证结构设计的正确性，可显著缩短结构设计周期。

Description

一种基于等几何粒子描述的拓扑优化方法

技术领域

本发明涉及复杂物理场拓扑优化技术领域，具体涉及一种基于等几何粒子描述的拓扑优化方法。

背景技术

电力电子设备的集成设计、航空航天装备设计、核能设备设计、微纳流控结构、热控结构设计等复杂高端装备技术领域经常涉及复杂耦合物理场的求解及结构优化设计问题。对于复杂物理场的求解，现有的数值分析手段主要以有限元或有限体等欧拉网格法和光滑粒子流体动力学无网格法为主，此类数值分析手段尽管已经被广泛使用，并有大量的商业软件，但对于复杂耦合物理场的计算，欧拉网格法耗费时间长，且计算结果的准确程度严重依赖于选用的数值分析模型以及网格的数量和质量，且由于方法原理上的限制，欧拉网格法无法准确描述物理场的发展过程；而光滑粒子流体动力学方法计算精度受离散程度限制，难以同时保持计算精度和计算效率。

电子设备的集成化、航空航天装备的轻量化，核能设备的可靠性，微纳流控、热控结构的微小化都对结构的性能提出了更高的要求，要求在较为恶劣的工况下和极为狭小的空间内，最大程度的发挥结构的性能这对相关领域的结构设计工作带来了更加困难的挑战。以往设计者根据工程经验的指导设计相应的结构再进行实验测试，验证其是否满足使用需求，然后给出改进措施。这样的设计方法虽然有工程实际作为指导，且可以较为简便的对设计进行优化。但是在结构设计功构一体化的发展需求下，上述设计模式无法在设计之初保证结构设计的正确性，难以快速高效的挖掘出结构设计潜能，不能发挥材料布局的最大能力，缺乏科学完善的设计理论依据。

发明内容

为了克服上述现有技术的缺点，本发明的目的在于提供一种基于等几何粒子描述的拓扑优化方法，在设计之初能保证结构设计的正确性，可显著缩短结构设计周期。

为达到上述目标，本发明采取的技术方案为：

一种基于等几何粒子描述的拓扑优化方法，采用等几何粒子流体动力学方法进行物理场仿真计算，结合基于变密度法SIMP的拓扑表征方法完成结构优化设计，根据拉格朗日观点思想和非定常流场的流场仿真的结果，完成目标函数关于设计变量的灵敏度分析，迭代优化设计变量，从而推动目标函数按照设计意图变化，促使结构向预定性能逼近，最终完成给定约束下，材料布局和结构设计的最优组合。

一种基于等几何粒子描述的拓扑优化方法，包括以下步骤：

1)建立基于SIMP法的拓扑优化模型：

1.1)建立SIMP法结构表征模型：

根据设计域的真实CAD模型，建立设计域的欧拉离散网格；采用SIMP拓扑优化方法，以各向同性的伪密度单元为结构优化的最小单位对设计域进行离散；SIMP法以单元的伪密度值γ作为设计变量，其反映了材料密度与材料性质之间对应关系；伪密度值的1和0分别代表了该位置结构的有无，设计变量场γ＝{γ₁,γ₂,...,γ_i,...}^T表征了设计域中的结构分布；

1.2)不同相/场的统一描述：

设计域内包含不同伪密度值的伪密度单元，其所代表的相/场对于物理场的作用程度与伪密度值相关，并作为额外项添加到物理场控制方程中；伪密度值为0和1的单元对于物理场的作用效果容易计算，但是中间值的伪密度单元的作用则需要通过设计相关量反映到物理场控制方程中；

1.3)建立SIMP拓扑优化数学模型：

考虑约束优化问题的形式：

min J(u,γ)

s.t.e(u,γ)＝0,inΩ

c(u,γ)≤0

0≤γ≤1

其中u为状态变量，γ为设计变量，e和c分别为物理控制方程和约束条件，至此，建立了基于SIMP法的拓扑优化模型；

2)等几何粒子流体动力学物理场仿真计算：

2.1)确定待求解的计算域和设计域：

根据真实几何模型，简化提取出计算域和设计域；引入非均匀有理B样条曲线NURBS基函数，构建二阶B样条曲线，对计算域进行划分；在参数空间内定义控制计算域几何模型的曲线组及相应的控制节点向量组；根据边界曲线组定义二维NURBS曲面，从而完整表征出计算域和设计域；

2.2)建立待求解物理问题的分析模型：

在步骤1.1)的基础上，对物理问题的计算域进行双线性离散，并对计算域内的控制粒子赋予信息；根据NURBS基函数的特点，对同一曲线上的控制粒子信息的线性组合，构建待求解物理问题的控制粒子层离散模型；控制粒子层的参数坐标与其物理空间坐标的转换关系采用基于NURBS基函数的方法插值：

其中R_i,j为双线性NURBS基函数，c_i,j为控制粒子坐标向量，η,ξ,ζ为参数坐标系下的坐标，F为基于NURBS基函数的坐标转换；

物理粒子层建立在控制粒子层之上，物理粒子所承载的物理信息由控制粒子物理信息和NURBS基函数插值得到：

其中，p_i,j为物理粒子承载的待求解物理场信息，d_i,j为控制粒子承载的物理信息；

2.3)计算域边界条件施加：

提前标记边界控制粒子和内部控制粒子；求解时也将解向量分解为边界控制解和内部控制解；

边界条件通过NURBS基函数，基于边界力方法施加在边界层控制粒子上，边界条件包括壁面碰撞边界条件或热力学条件；

2.4)物理粒子相互作用：

对物理粒子施加相应的物性参数，控制粒子不承载物理属性；结合有相位移场思想，通过设置粒子间干涉半径，粒子间相互作用方程的形式，计算并施加物理粒子间的相互作用；

2.5)求解非定常物理场信息：

求解非定常偏微分控制方程获得待求解问题的物理场分布，对于方程中的高阶偏微分项，基于NURBS基函数结合雅可比矩阵对其处理如下：

采用瞬态分析方法，计算从初始时刻到某一时刻的待求解物理场的变化过程；以给定的时间步长为单位，逐步求解未知物理场控制粒子信息的发展情况，再根据步骤2)的映射插值关系计算物理粒子信息，得到的非定常物理场分布逐步计算形式如下所示：

其中为第l步物理粒子所承载的待求解物理场信息，/>为第l-1步物理粒子所承载的待求解物理场信息，/>为待求解物理场信息关于时间的变化发展速率，Δt为给定的时间步长；

3)灵敏度分析：

3.1)设计域的欧拉离散：

离散的网格数量、形状、大小与建立的基于SIMP法的拓扑优化模型对应；

3.2)非定常场的物理粒子信息输出：

非定常场的物理粒子信息输出分为两级输出：分别是粒子物理信息的欧拉投影输出和非定常场的物理信息输出；

所述的粒子物理信息的欧拉投影输出含义为粒子向欧拉离散网格内投影，根据粒子的位置坐标将其分散到各个欧拉离散网格中，粒子依据位置分组，从属于所在的网格；同一网格内粒子上携带的物理场信息数据取平均值作为该网格的物理场数据值，计算方法如下列公式所示：

其中，A^e是欧拉离散网格的物理场数值，A^p是流场粒子的物理信息，i为欧拉离散网格编号，n_i为第i个网格内包含的粒子总数；

3.3)敏度公式推导：

采用非定常思想，逐帧求解物理场敏度，拓扑优化的一般数学模型为：

min J(u,γ)

s.t.e(u,γ)＝0,inΩ

c(u,γ)≤0

0≤γ≤1

所述的敏度分析指求解目标函数J对设计变量γ的导数，物理场的非定常求解要求敏度分析同样为非定常形式，即逐帧求解：

取流动稳定后的一段时间的逐帧敏度平均值作为驱动设计变量变化的最终敏度，即：

其中，f为计算敏度所用的帧数；

4)设计变量迭代优化：

使用移动渐近线算法对数学模型进行优化，直至目标函数收敛，获得优化结构；

5)结构后处理：

结合加工工艺要求或者装配要求对得到的优化结构进行后续的圆整处理，得到最终的结构优化设计。

所述的步骤1.2)中设计相关量有材料的渗透率、有效导热系数、人工材料模量。

所述的步骤2.5)中求解非定常偏微分控制方程获得待求解问题的物理场分布，NURBS插值基函数支持控制场量信息关于位置的高阶偏导项，处理过程涉及到参数空间和物理空间之间的空间转换，首先完成基函数在参数空间内的位置偏导项计算，再使用雅克比矩阵实现参数空间位置偏导项向物理空间位置偏导项的映射方程中的高阶偏微分项。

所述的步骤3.2)中所述的非定常场的物理信息输出是指由于等几何粒子流体动力学方法是一种拉格朗日观点下的非定常物理场求解方法，为得到用于拓扑优化的稳态物理场数据，输出物理状态稳定后一定帧数的整合后的欧拉投影物理场数据，物理场数据以帧数为单位逐帧输出。

本发明的有益效果是：本发明以一种拉格朗日观点下的等几何粒子流体动力学方法，实现了复杂耦合物理场的高效求解；结合基于SIMP法的拓扑表征方法，提出了一种基于等几何粒子描述的拓扑优化方法，将物理场分析与结构设计结合起来，摆脱了设计人员的经验依赖，为设计过程提供了充足的理论依据，使设计结果更加合理可靠，性能更为优异，充分发挥材料潜能。本发明应用范围广泛，可有效解决电力电子设备的集成设计、航空航天装备设计、核能设备设计、微纳流控结构、热控结构设计等复杂高精尖技术领域的结构设计优化难题，缩短了装备升级周期。

附图说明

图1是本发明的实施例风道散热结构模型图。

图2是本发明的流程图。

图3是本发明实施例的初始拓扑模型图。

图4是本发明实施例的等几何粒子流体动力学计算模式图。

图5是本发明实施例的计算域投影网格示意图。

图6是本发明实施例的粒子投影算法示意图。

图7是本发明实施例的迭代优化后的风道结构示意图。

具体实施方式

下面结合实施例和附图对本发明做进一步说明。

如图1所示，本发明的实施例为某型号电子设备的风冷流道设计，设计区域长160mm，宽120mm；结构设计的条件为两个方形热源11的热流密度为2.75W/cm²，两个入口12流速5m/s，两个出口13的压强为大气表压，流体工质为空气，热源材料为铝，设计目标为尽量降低热源的温度。

如图2所示，一种基于等几何粒子描述的拓扑优化方法，包括以下步骤：

1)建立基于SIMP法的拓扑优化模型：

1.1)建立SIMP法结构模型：

如图3所示，根据风冷流道设计域的形状，建立设计域的欧拉离散网格，得到了初始的SIMP结构模型；采用SIMP拓扑优化方法，以各向同性的伪密度单元为结构优化的最小单位对设计域进行离散；SIMP法以单元的伪密度值γ作为设计变量，其反映了单元的相对密度与材料性质之间对应关系；伪密度值的1和0分别代表了该位置材料的有无，设计变量场γ＝{γ₁,γ₂,...,γ_i,...}^T表征了设计域中的材料分布；

1.2)不同相/场的统一描述：

设计域内包含不同伪密度值的伪密度单元，其所代表的固体材料对于物理场的作用程度与伪密度值相关，并作为额外项添加到控制方程中；伪密度值为0和1的单元对于物理场的作用效果可以较为容易的计算出来，但是中间值的伪密度单元的作用则需要通过设计相关量反映到物理场控制方程中；在本实施例中，选择材料的渗透率作为设计相关量，渗透率是多孔介质材料的属性，在此处表示固体材料对气体流动的阻碍作用，渗透率与单元伪密度的关系如下所示：

其中α_max是渗透率的最大值，取为720，α_min为渗透率的最小值，取为0，q为一正实数，取为100；

1.3)建立SIMP拓扑优化数学模型：

本实施例的风冷流道结构优化目标为尽量降低热源温度。约束为材料总用量不高于设计域体积的40％。根据该目标和约束，建立如下的拓扑优化数学模型：

0≤γ_j≤1

其中T_hi为热源处的网格温度数值，n_e为设计变量的数量，α为单元的渗透率，n_ph为热源处的粒子数量；至此，建立了基于SIMP法的拓扑优化模型；

2)等几何粒子流体动力学物理场仿真计算：

如图4所示，本发明的一种等几何粒子描述的拓扑优化方法首先要进行等几何粒子流体动力学物理场仿真计算，其核心思想是通过等几何NURBS基函数实现控制粒子与物理粒子的引导式离散，提高降低参与计算的粒子数量，提高计算效率且不损失精度；

2.1)确定待求解的计算域和设计域：

根据风冷流道计算域的几何形状，简化提取出计算域和设计域；引入非均匀有理B样条曲线(NURBS)基函数，构建二阶B样条曲线，对计算域进行划分；B样条基函数的递推公式如下：

其中N为B样条曲线基函数，p为基函数阶次，ξ为参数坐标系下的节点；

在参数空间内定义控制风冷流道计算域几何模型的曲线组及相应的控制节点向量组；根据边界曲线组定义二维NURBS曲面，从而完整表征出计算域和设计域；NURBS基函数表达形式如下：

其中R为双线性NURBS基函数，N为B样条曲线基函数，p为基函数阶次，ω为投射投影权因子；

2.2)建立待求解物理问题的分析模型：

在步骤1.1的基础上，对风冷流道的计算域进行双线性离散，并对计算域内的控制粒子赋予信息；根据NURBS基函数的特点，对同一曲线上的控制粒子信息的线性组合，构建风冷流道中冷却工质的控制粒子层离散模型；控制粒子层的参数坐标与其物理空间坐标的转换关系采用基于NURBS基函数的方法插值：

其中R_i,j为双线性NURBS基函数，c_i,j为控制粒子坐标向量，η,ξ,ζ为参数坐标系下的坐标，F为基于NURBS基函数的坐标转换；

物理粒子层建立在控制粒子层之上，物理粒子所承载的物理信息由控制粒子物理信息和NURBS基函数插值得到：

其中，p_i,j为物理粒子承载的待求解的流速、温度、压强等信息，d_i,j为控制粒子承载的信息；

2.3)计算域边界条件施加：

对于风冷流道这样的耦合复杂物理场而言，边界层物理控制方程与内部物理控制方程有所区别，故需要提前标记边界控制粒子和内部控制粒子；求解时也将解向量分解为边界控制解和内部控制解；

边界条件通过NURBS基函数，基于边界力方法施加在边界层控制粒子上，边界条件包括壁面碰撞边界条件和热力学条件；

2.4)物理粒子相互作用：

风冷流道物理场的求解由气体属性粒子和热源处的固体粒子共同参与；对气体粒子施加气体物性参数，本实施例中的气体为空气，密度设置为1.205kg/m³，导热系数0.026W/m·k，恒压热容为1.013kJ/kg·k；对固体粒子施加固体物性参数，本实施例中的固体为铝，密度为2.719×10³kg/m³，导热系数为202.4W/m·k，恒压热容为871J/kg·k；控制粒子不承载物理属性；结合有相位移场思想，通过设置粒子间干涉半径，粒子间相互作用方程的形式，计算并施加物理粒子间的相互作用；

2.5)求解非定常物理场信息：

求解流体N-S方程和导热方程从而得到流道的物理场分布，控制方程形式为：

粒子离散形式的控制方程为：

其中ρ为气体密度，μ为运动粘度，k为空气的导热系数，C_p为恒压热容，v_i为粒子速度，T_i为粒子温度；

NURBS插值基函数支持控制场量信息关于位置的高阶偏导项，处理过程涉及到参数空间和物理空间之间的空间转换，首先完成基函数在参数空间内的位置偏导项计算，再使用雅克比矩阵实现参数空间位置偏导项向物理空间位置偏导项的映射方程中的高阶偏微分项；计算规则如下

依据上述规则，对于控制方程中的高阶偏微分项基于NURBS基函数结合雅可比矩阵对其处理如下：

采用瞬态分析方法，计算从初始时刻到某一时刻的速度场和温度场的变化过程；本实施例的每帧时间为1/120秒，逐帧求解速度场和温度场控制粒子数值的发展情况，再根据上述的映射插值关系计算物理粒子信息；非定常速度场和温度场分布逐步计算形式如下所示：

其中T_i ^l为第l步物理粒子所承载的速度和温度信息，/>为第l-1步物理粒子所承载的速度和温度信息，/>为速度和温度信息关于时间的变化发展速率，Δt为给定的每帧时间；

3)灵敏度分析：获取物理场信息后进行敏度分析，求解目标函数约束函数对于目标的敏度，从而使结构趋优；

3.1)设计域的欧拉离散：

结合了拉格朗日和欧拉两种观点，设计域的离散形式是欧拉观点下的网格法，因此在计算敏度从而驱动伪密度单元变化时，也必须以单元作为基本单位，故必须对设计域进行欧拉离散，离散的网格数量、形状、大小与建立的SIMP拓扑优化模型对应；

3.2)非定常速度场和温度场的物理粒子信息输出：

对于拓扑优化而言，敏度计算离不开物理场状态信息；非定常场的物理粒子信息输出分为两级输出：分别是粒子物理信息的欧拉投影输出和非定常场的物理信息输出；

所述的粒子物理信息的欧拉投影输出含义为粒子向欧拉离散网格内投影，根据粒子的位置坐标将其分散到各个欧拉离散网格中，粒子依据位置分组，从属于所在的网格，本实施例的风冷流道的投影网格如图5所示；

如图6所示，本发明的粒子投影算法指同一网格内粒子上携带的物理场信息数据取平均值作为该网格的物理场数据值，本实施例的速度场和温度场网格数值计算方法如下列公式所示：

其中，v^e是欧拉离散网格的速度场数值，T^e是欧拉离散网格的温度场数值，j为欧拉离散网格编号，为第j个网格内包含的粒子总数；

所述的非定常场的物理信息输出是指由于等几何粒子流体动力学方法是一种拉格朗日观点下的非定常物理场求解方法，为得到用于拓扑优化的稳态物理场数据，输出物理状态稳定后一定帧数的整合后的欧拉投影物理场数据。在本实施例中，以帧数为单位逐帧输出物理粒子的速度场投影数据和温度场投影数据。

3.3)敏度公式推导：

等几何粒子流体动力学是典型的非定常物理场求解模式，在敏度求解时也采用非定常思想，逐帧求解物理场敏度。风冷流道的数学模型为：

0≤γ_j≤1

其中l为粒子的位置，n为设计变量的个数。

所述的敏度分析指求解目标函数J对设计变量γ的导数，物理场的非定常求解要求敏度分析同样为非定常形式，即逐帧求解；本实施例是典型的强迫对流换热场景，热源温度与受流场影响，设计变量以渗透率的形式带入控制方程，且前述的两者之间表达式十分清晰，故求解目标函数热源温度对设计变量γ的敏度关键在于求解温度场关于渗透率的敏度/>本实施例温度场关于渗透率的敏度具体求解方式如下：

首先根据控制方程求解粒子温度对于设计变量的敏度：

其中，T_i1是粒子当前的速度，T_i0是粒子前一帧的速度，Δt为两帧之间的时间间隔，i为粒子编号，由链式法则可知：

同理有：

综合以上各式，将目标函数温度的敏度通过链式法则转换为速度场的敏度，速度场的敏度按照如下方式求解：

其中，v_i1是粒子当前的速度，v_i0是粒子前一帧的速度，Δt为两帧之间的时间间隔，α_i是该粒子所在SIMP单元的渗透率；其中的微分项可由上述的基于NURBS基函数结合雅可比矩阵进行处理，至此，求解出了粒子温度关于渗透率的敏度渗透率关于设计变量的敏度为：

由链式法则有：

再按照粒子投影形式将每个粒子的敏度投影到各个网格中即求得了每个单元的温度场敏度；

取流动稳定后的一段时间的逐帧敏度平均值作为驱动设计变量变化的最终敏度，即：

其中，f为计算敏度所用的帧数；

4)设计变量迭代优化：

使用移动渐近线算法对数学模型进行优化，直至目标函数收敛，收敛条件为相邻两次迭代的目标函数差值小于0.001，获得优化结构；

5)结构后处理：

结合加工工艺要求或者装配要求对得到的优化结构进行后续的圆整处理，得到最终的结构优化设计，如图7所示。

Claims (4)

1.一种基于等几何粒子描述的拓扑优化方法，其特征在于：采用等几何粒子流体动力学方法进行物理场仿真计算，结合基于变密度法SIMP的拓扑表征方法完成结构优化设计，根据拉格朗日观点思想和非定常流场的流场仿真的结果，完成目标函数关于设计变量的灵敏度分析，迭代优化设计变量，从而推动目标函数按照设计意图变化，促使结构向预定性能逼近，最终完成给定约束下，材料布局和结构设计的最优组合；

所述的一种基于等几何粒子描述的拓扑优化方法，包括以下步骤：

1)建立基于SIMP法的拓扑优化模型：

1.1)建立SIMP法结构表征模型：

根据设计域的真实CAD模型，建立设计域的欧拉离散网格；采用SIMP拓扑优化方法，以各向同性的伪密度单元为结构优化的最小单位对设计域进行离散；SIMP法以单元的伪密度值γ作为设计变量，其反映了材料密度与材料性质之间对应关系；伪密度值的1和0分别代表了该位置结构的有无，设计变量场γ＝{γ₁,γ₂,...,γ_i,...}^T表征了设计域中的结构分布；

1.2)不同相/场的统一描述：

设计域内包含不同伪密度值的伪密度单元，其所代表的相/场对于物理场的作用程度与伪密度值相关，并作为额外项添加到物理场控制方程中；伪密度值为0和1的单元对于物理场的作用效果容易计算，但是中间值的伪密度单元的作用则需要通过设计相关量反映到物理场控制方程中；

1.3)建立SIMP拓扑优化数学模型：

考虑约束优化问题的形式：

min J(u,γ)

s.t.e(u,γ)＝0,inΩ

c(u,γ)≤0

0≤γ≤1

其中u为状态变量，γ为设计变量，e和c分别为物理控制方程和约束条件，至此，建立了基于SIMP法的拓扑优化模型；

2)等几何粒子流体动力学物理场仿真计算：

2.1)确定待求解的计算域和设计域：

根据真实几何模型，简化提取出计算域和设计域；引入非均匀有理B样条曲线NURBS基函数，构建二阶B样条曲线，对计算域进行划分；在参数空间内定义控制计算域几何模型的曲线组及相应的控制节点向量组；根据边界曲线组定义二维NURBS曲面，从而完整表征出计算域和设计域；

2.2)建立待求解物理问题的分析模型：

在步骤1.1)的基础上，对物理问题的计算域进行双线性离散，并对计算域内的控制粒子赋予信息；根据NURBS基函数的特点，对同一曲线上的控制粒子信息的线性组合，构建待求解物理问题的控制粒子层离散模型；控制粒子层的参数坐标与其物理空间坐标的转换关系采用基于NURBS基函数的方法插值：

其中R_i,j为双线性NURBS基函数，c_i,j为控制粒子坐标向量，η,ξ,ζ为参数坐标系下的坐标，F为基于NURBS基函数的坐标转换；

物理粒子层建立在控制粒子层之上，物理粒子所承载的物理信息由控制粒子物理信息和NURBS基函数插值得到：

其中，p_i,j为物理粒子承载的待求解物理场信息，d_i,j为控制粒子承载的物理信息；

2.3)计算域边界条件施加：

提前标记边界控制粒子和内部控制粒子；求解时也将解向量分解为边界控制解和内部控制解；

边界条件通过NURBS基函数，基于边界力方法施加在边界层控制粒子上，边界条件包括壁面碰撞边界条件或热力学条件；

2.4)物理粒子相互作用：

对物理粒子施加相应的物性参数，控制粒子不承载物理属性；结合有相位移场思想，通过设置粒子间干涉半径，粒子间相互作用方程的形式，计算并施加物理粒子间的相互作用；

2.5)求解非定常物理场信息：

求解非定常偏微分控制方程获得待求解问题的物理场分布，对于方程中的高阶偏微分项，基于NURBS基函数结合雅可比矩阵对其处理如下：

采用瞬态分析方法，计算从初始时刻到某一时刻的待求解物理场的变化过程；以给定的时间步长为单位，逐步求解未知物理场控制粒子信息的发展情况，再根据步骤2)的映射插值关系计算物理粒子信息，得到的非定常物理场分布逐步计算形式如下所示：

其中为第l步物理粒子所承载的待求解物理场信息，/>为第l-1步物理粒子所承载的待求解物理场信息，/>为待求解物理场信息关于时间的变化发展速率，Δt为给定的时间步长；

3)灵敏度分析：

3.1)设计域的欧拉离散：

离散的网格数量、形状、大小与建立的基于SIMP法的拓扑优化模型对应；

3.2)非定常场的物理粒子信息输出：

非定常场的物理粒子信息输出分为两级输出：分别是粒子物理信息的欧拉投影输出和非定常场的物理信息输出；

所述的粒子物理信息的欧拉投影输出含义为粒子向欧拉离散网格内投影，根据粒子的位置坐标将其分散到各个欧拉离散网格中，粒子依据位置分组，从属于所在的网格；同一网格内粒子上携带的物理场信息数据取平均值作为该网格的物理场数据值，计算方法如下列公式所示：

其中，A^e是欧拉离散网格的物理场数值，A^p是流场粒子的物理信息，i为欧拉离散网格编号，n_i为第i个网格内包含的粒子总数；

3.3)敏度公式推导：

采用非定常思想，逐帧求解物理场敏度，拓扑优化的一般数学模型为：

min J(u,γ)

s.t.e(u,γ)＝0,inΩ

c(u,γ)≤0

0≤γ≤1

所述的敏度分析指求解目标函数J对设计变量γ的导数，物理场的非定常求解要求敏度分析同样为非定常形式，即逐帧求解：

取流动稳定后的一段时间的逐帧敏度平均值作为驱动设计变量变化的最终敏度，即：

其中，f为计算敏度所用的帧数；

4)设计变量迭代优化：

使用移动渐近线算法对数学模型进行优化，直至目标函数收敛，获得优化结构；

5)结构后处理：

结合加工工艺要求或者装配要求对得到的优化结构进行后续的圆整处理，得到最终的结构优化设计。

2.根据权利要求1所述的方法，其特征在于：所述的步骤1.2)中设计相关量有材料的渗透率、有效导热系数、人工材料模量。

3.根据权利要求1所述的方法，其特征在于：所述的步骤2.5)中求解非定常偏微分控制方程获得待求解问题的物理场分布，NURBS插值基函数支持控制场量信息关于位置的高阶偏导项，处理过程涉及到参数空间和物理空间之间的空间转换，首先完成基函数在参数空间内的位置偏导项计算，再使用雅克比矩阵实现参数空间位置偏导项向物理空间位置偏导项的映射方程中的高阶偏微分项。

4.根据权利要求1所述的方法，其特征在于：所述的步骤3.2)中所述的非定常场的物理信息输出是指由于等几何粒子流体动力学方法是一种拉格朗日观点下的非定常物理场求解方法，为得到用于拓扑优化的稳态物理场数据，输出物理状态稳定后一定帧数的整合后的欧拉投影物理场数据，物理场数据以帧数为单位逐帧输出。