CN108345741B - 基于无网格rkpm各向异性材料二维热变形和热应力分析方法 - Google Patents

Abstract

本发明公开了基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法，主要包括以下步骤：(1)对计算模型进行RKPM节点离散前处理；(2)求解RKPM热变形位移值：组装RKPM整体力刚度矩阵和整体热载荷列向量；施加边界，采用罚函数法处理第一类边界条件；建立各向异性材料结构无网格RKPM热应力离散控制方程，并求解节点的RKPM热变形位移参数值；(3)求解RKPM热应力值，采用重构核近似对所得热变形位移参数值进行逼近，可计算出高斯点的热应力，进而得到节点的RKPM热应力值；(4)对计算结果进行后处理。本发明基于无网格RKPM进行各向异性材料结构的二维热变形和热应力分析，数值方法稳定、精度高。

Description

基于无网格RKPM各向异性材料二维热变形和热应力分析方法

技术领域

本发明属于计算机辅助工程中的数值传热学领域，具体涉及一种基于无网格重构核粒子法(Reproducing Kernel Particle Method,RKPM)的各向异性材料结构二维热变形和热应力分析方法。

背景技术

热应力广泛存在于机械工程和动力工程实际问题中，它能造成设备损坏，甚至引起严重的生产事故和巨大的经济损失。当物体温度改变时，由于外在约束以及内部各部分之间的相互约束，使其不能完全自由胀缩而产生的应力，称为热应力，又称为变温应力。发动机活塞顶部区域、汽轮机转子、化工塔设备的裙座支撑区和重轨等，这些地方都是热应力产生的高发区，容易产生热应力失效，破坏设备的正常功能。因此，对这些热应力高发区域进行热变形和热应力分析，并采取相应的控制措施，才能最大限度地减少热应力造成的不良后果。

复合材料是一种混合物，是由两种或两种以上不同性质的材料，通过物理或化学的方法，在宏观(微观)上组成具有新性能的材料。各种材料在性能上互相取长补短，产生协同效应，使复合材料的综合性能优于原组成材料而满足各种不同的要求。与传统材料相比，复合材料具有重量轻、强度高、加工成型方便、弹性优良、耐化学腐蚀和耐候性好等特点，已逐步取代金属合金等传统材料，广泛应用于机械工程、航空航天、汽车工程和电子电气等领域。除此之外，复合材料的一个显著特点是各向异性，在不同方向具有不同的热、力性能，这种传热性能和力学性能的差异不仅与材料本身的结构形式有关，还与材料布局有关。因此，各向异性材料的导热和力学性能具有明显地独立性，即各向异性导致同一部位不同方向的传热能力和受力情况不同。

目前，热应力的数值计算方法主要有有限元法、有限体积法和有限差分法等，但它们都是基于网格的数值计算方法，网格划分过程繁琐，耗费较大，而且需要高质量的网格，不利于分析复杂几何形状的零部件。另外，它们在处理大变形、裂纹扩展、相变和动边界等问题时，容易因网格畸变或单元奇异等问题而造成较大计算误差，甚至求解失败。为了克服有限元法等依赖网格划分技术的传统数值计算方法的不足，学者们提出了无网格法。无网格法作为一种较新的数值计算方法，在构造形函数时仅依靠一系列可灵活布置的节点，摆脱了网格束缚，简化了前处理，并且无单元连续性要求，目前已有很多种不同形态的无网格方法，包括：无网格Galerkin法(Element-free Galerkin Method，EFGM)、重构核粒子法(Reproducing Kernel Particle Method，RKPM)、无网格局部Petrov-Galerkin法(Meshless Local Petrov-Galerkin Method，MLPG)等。无网格重构核粒子法(无网格RKPM)是众多无网格法中的一种，它具有其它无网格方法所不具备的多分辨率和变时频特性的优势，其应用领域也更加广阔。

目前，无网格法在工程热应力问题的数值计算中已得到了一些应用。苑素玲等采用EFGM求解了汽轮机转子的热传导和热应力问题；洪新兰等利用EFGM研究了正方形板结构的热变形和热应力问题；王峰等利用改进的EFG-SBM法(即无网格伽辽金法和边界有限元法相结合的新方法)研究了连续体的热应力问题，但上述研究大多局限于各向同性材料或具有简单几何形状的结构热应力问题。张建平和周国强等虽然利用EFGM和无网格RKPM分别研究了二维各向异性材料结构的温度场分布和二维各向异性材料热结构拓扑优化设计，但并没有考虑更深层次的热变形和热应力问题，即没有研究热边界和位移边界共同约束、热载荷和机械载荷共同作用下的热变形和热应力问题，只是单纯地考虑了材料热导率正交各向异性因子和各向异性材料方向角对二维各向异性材料结构的传热性能的影响，即对温度场分布规律的影响，而没有进一步开展基于无网格RKPM的二维各向异性材料复杂结构的热变形和热应力研究。另外，无网格法目前尚无通用的商业软件，因此学者们都是首先通过算法公式推导，然后自行编程进行无网格法的理论和应用研究。

发明内容

目前，能源动力、航空航天及微机械系统(MEMS)等领域的快速发展对材料和结构的综合性能提出了越来越高的要求，除了需要具有很好的力学性能和轻量化要求外，还需要具备优良的散热和热强度性能等。各向异性的复合材料因具有重量轻、强度高、加工成型方便、弹性优良、耐化学腐蚀和耐候性好等特点，在能源动力、航空航天和机械工程等众多工程领域得到了广泛应用，并且在很多领域代替了传统材料。为了提高求解各向异性材料复杂几何结构的热变形和热应力问题的计算精度，本发明基于无网格RKPM理论，结合热-力学的热变形和热应力知识，采用MATLAB语言进行编程，编写了基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法的程序，该方法摆脱了网格束缚，简化了前处理过程，避免了传统的基于有限元法的热变形和热应力分析可能出现的网格畸变、单元奇异和计算精度较低的问题，适用于求解工程中各向异性材料复杂几何结构的热变形和热应力问题。

本发明解决其技术问题所采用的技术方案是：采用MATLAB语言编程实现基于无网格RKPM的各向异性材料结构二维热变形和热应力分析，其通过热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt和材料方向角θ来控制正交各向异性材料的热变形能力和热应力性能，以灵活地实施不同各向异性材料结构的无网格RKPM热变形和热应力分析；采用平面转轴矩阵、正交各向异性总弹性矩阵、热膨胀弹性矩阵将各向异性材料坐标系1—2的热应力与应变关系转换为与几何坐标系X—Y相一致的热应力与应变关系，具体变换为：

上式中

是坐标变换矩阵，α₁，α₂分别为材料坐标系主轴1，2方向的热膨胀系数，定义热膨胀正交各向异性因子Pt＝α₁/α₂，

是正交各向异性弹性矩阵，

定义主次泊松比因子Bt＝ν₁₂/ν₂₁，λ₁和λ₂分别代表材料坐标系主轴1，2方向的热导率，定义材料的热导率正交各向异性因子Ht＝λ₁/λ₂，通过修改热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt和材料方向角θ就可进行不同的各向异性材料结构的无网格RKPM热变形和热应力分析；然后对计算模型进行无网格RKPM前处理，再对各向异性材料结构进行无网格RKPM热变形分析，将无网格RKPM计算模型的实际温度与参考温度的温差引起的温差力载荷施加到模型上，并施加位移边界条件以及力载荷，获得计算模型的无网格RKPM热变形位移场，将获得的无网格RKPM热变形位移代入应力计算公式得到高斯点的无网格RKPM应力值，利用得到的高斯点无网格RKPM应力求得节点无网格RKPM应力，最后将计算结果进行后处理。

本发明所述技术方案的具体实施步骤如下：

(1)对计算模型进行无网格RKPM前处理：根据实际工程中所选计算模型的工况，确定无网格RKPM各向异性材料结构的主热导率λ₁、主弹性模量E₁、主泊松比ν₁₂、主热膨胀系数α₁、剪切模量G₁₂、热导率正交各向异性因子Ht、主次泊松比因子Bt、热膨胀正交各向异性因子Pt和材料方向角θ等材料属性；根据给定的材料属性，计算热导率坐标变换矩阵、正交各向异性弹性矩阵、平面转轴矩阵、正交各向异性总弹性矩阵、热膨胀弹性矩阵和热膨胀系数列向量；导入计算模型的几何模型、无网格RKPM离散节点信息和边界条件；生成无网格RKPM积分背景网格；确定高斯点的坐标、权重系数和雅可比；

(2)基于无网格RKPM分析各向异性材料结构的二维热变形位移值，无网格RKPM区别于其他无网格法最大的不同在于其形函数采用重构核近似进行构造，在任意节点x处的热变形位移场u^h(x)可由其影响域内的节点值u_I拟合出来，

式中，C(x；x_I-x)为校正函数，ω(x_I-x)为核函数，Δx_I为第I个节点x_I所对应的面积，φ_I(x)为对应第I个节点x_I的无网格RKPM形函数，写成矩阵形式为：

Φ(x)＝[φ₁(x),φ₂(x)···φ_n(x)]＝p^T(0)M^-1(x)H(x) (3)

其中，

p^T(0)＝[1,0,0,0,0,0] (4)

H(x)＝[p(x₁-x)ω(x₁-x)Δx₁,p(x₂-x)ω(x₂-x)Δx₂,...,p(x_I-x)ω(x_I-x)Δx_I] (6)

p(x_I-x)＝[1,x_I-x,y_I-y,(x_I-x)²,(x_I-x)(y_I-y),(y_I-y)²]^T (7)

基于无网格RKPM的热变形位移场分析的具体步骤为：(A)计算单个高斯点影响域内的RKPM力刚度矩阵并组装整体RKPM力刚度矩阵；(B)计算单个高斯点影响域内由于温差引起的RKPM温差载荷向量并组装整体RKPM温差载荷列向量，温差为该高斯点的RKPM实际温度与RKPM参考温度的差值，高斯点的RKPM实际温度求解具体步骤为：(a)首先根据节点和高斯点坐标求各高斯点与节点之间的距离并从小到大排序，并规定影响域内的节点个数，一般取8～10个为宜，按照距离先后顺序采用可视性准则判断RKPM温度场的非连续性，即节点是否跨越边界，如果跨界则删掉该节点，并补充距离高斯点最近的新的节点进入，再次进行可视性准则判断，直到满足条件的节点达到规定的节点个数，取符合条件的节点中距离高斯点最远的节点和该高斯点之间的距离为s[k]，该高斯点的动态影响域半径d_mI＝scale×s[k]，scale是大于1的乘子，取值为1.2～1.5，采用圆形影响域，核函数ω_I(x)＝ω(x-x_I)，在第I个节点x_I处的值最大，且具有紧支性，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mI，将r定义为

其中x_I，y_I为第I个节点x_I在x方向和y方向的节点坐标；其次逐个搜索高斯点影响域内的节点并计算其重构核近似形函数及其导数，建立单个高斯点影响域内的RKPM热传导矩阵并组装整体RKPM热传导刚度矩阵；(b)施加混合传热边界条件，采用罚函数法处理第一类传热边界条件，罚因子取值范围为2×10⁶～2×10⁸；(c)形成系统总的RKPM热刚度矩阵和RKPM温差载荷列向量；(d)建立各向异性材料结构无网格RKPM稳态传热离散控制方程，并求解计算模型离散节点的RKPM温度参数值；(e)采用重构核近似对得到的节点RKPM温度参数值进行插值，就可获得的计算模型高斯点的RKPM实际温度值；(C)施加位移边界条件和力载荷，采用罚函数法处理位移边界条件，罚因子的经验取值为10³×E，其中E为主弹性模量；(D)建立各向异性材料结构的无网格RKPM热应力离散控制方程，并求解计算模型离散节点的RKPM热变形位移参数值，具体方程为：

Ku＝f (8)

其中，

其中，D₁＝[T]^-1[Q]([T]^-1)^T，D₂＝[T]^-1[Q] (11)

式中：

其中，D₁为正交各向异性弹性矩阵，D₂为热膨胀弹性矩阵，[T]^-1为平面转轴矩阵，[Q]为正交各向异性弹性矩阵，φ_i,x和φ_i,y分别为无网格RKPM形函数关于坐标x和y的导数，∫_ΩΦ_if_idΩ为RKPM温差载荷列向量，

为处理位移边界条件引起的RKPM力载荷列向量修正项，α为罚因子，

为外力引起的RKPM力载荷列向量，当x或y方向有位移约束时，相应的s_x(或s_y)等于1；否则为0；θ为材料方向角；E₁和ν₁₂分别为材料沿x方向的弹性模量和泊松比，E₂和ν₂₁分别为材料沿y方向的弹性模量和泊松比，要保证E₁ν₂₁＝E₂ν₁₂；(E)采用重构核近似对得到的RKPM热变形位移参数值进行插值，就可获得计算模型离散节点的RKPM实际热变形位移值；

(3)求解无网格RKPM热应力值：利用步骤(2)得到的RKPM热变形位移参数值，计算出高斯点的RKPM应力值，进而得到计算模型离散节点的RKPM应力值，具体方程为

式中，ε₀为温度应变，对于平面问题温度应变为ε₀＝[α₁(θ_T-θ₀) α₂(θ_T-θ₀)0]^T，α₁，α₂分别为材料坐标系主轴1，2方向的热膨胀系数，θ₀为参考温度，θ_T为无网格RKPM计算模型高斯点的RKPM实际温度值，即基于无网格RKPM各向异性材料结构二维稳态热传导问题的解；求解无网格RKPM热应力值的具体步骤为：(A)首先根据输入高斯点和节点坐标求各节点与高斯点之间的距离并从小到大排序，并规定影响域内的高斯点个数，一般取9～12个为宜，按照距离先后顺序采用可视性准则判断RKPM热变形位移场的非连续性，即高斯点是否跨越边界，如果跨界则删掉该高斯点，并补充距离节点最近的新的高斯点进入，再次进行可视性准则判断，直到满足条件的高斯点达到规定的高斯点个数，取符合条件的高斯点中距离节点最远的高斯点和该节点之间的距离为为g[k]，该节点的动态影响域半径d_mx＝scale×g[k]，scale是大于1的乘子，取值为1.2～1.5，核函数ω_I(x)＝ω(x-x_I)，在第I个节点x_I处的值最大，且具有紧支性，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mx，将r定义为

其中gx，gy为高斯点坐标；(B)逐个搜索节点影响域内的高斯点，采用重构核近似对节点影响域内的高斯点的RKPM应力进行插值就可得到节点的RKPM应力值；

(4)对无网格RKPM热变形和热应力计算结果进行后处理，并将无网格RKPM的计算结果与有限元的计算结果以及参考解进行比较：(A)读入有限元的计算结果和参考解；(B)无网格RKPM得到的计算模型温度结果值输出及温度云图显示，无网格RKPM得到的计算模型热变形位移结果值的输出及位移云图显示，无网格RKPM得到的计算模型热应力结果值的输出及热应力云图显示；(C)利用输出的结果绘制无网格RKPM、有限元法和参考解的沿某一指定路径解的温度曲线图、热变形位移曲线图和热应力曲线图，并对三者进行直观的精度比较。

进一步地，所述基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法，步骤(2)和步骤(3)中采用圆形影响域会导致计算过程中产生分母为零的无穷大现象，为了防止除数为零，采用eps函数克服这一现象，具体为dis＝dis+(dis＝＝0)×eps，其中核函数采用三次样条函数，具体为：

进一步地，所述基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法，根据后处理结果可得到关于各向异性材料结构传热以及热变形和热应力的结论：无网格RKPM法的计算结果和有限元法的计算结果吻合得很好，且更逼近参考解，说明本方法是有效可行的；不同的热导率正交各向异性因子、热膨胀正交各向异性因子、主次泊松比因子和材料方向角对计算模型的RKPM热变形位移场和Mises应力场的分布产生了明显影响，热导率正交各向异性因子、热膨胀正交各向异性因子、主次泊松比因子主要影响RKPM热变形场和Mises应力场的强度，材料方向角主要影响RKPM热变形场和Mises应力场的分布方向。

本发明的有益效果是：本发明通过修改热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt和材料方向角θ就可灵活地进行不同的各向异性材料结构的无网格RKPM热变形和热应力分析；避免了基于有限元法的热变形和热应力分析可能出现的网格畸变、单元奇异或计算精度较低的问题，能更高效灵活地进行高斯点的布置，方法收敛速度快，数值稳定性好和计算精度高；本发明可求解不同的各向异性材料复杂结构的热变形和热应力问题，能与工程实践紧密结合，具有较好的理论研究和工程应用价值。

附图说明

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

图1是本发明各向异性材料结构的几何坐标系X—Y和材料坐标系1—2

图2是本发明无网格RKPM热变形和热应力分析流程框图

图3是本发明计算过程中无网格RKPM热变形和热应力分析程序调用关系图

图4是本发明算例重轨的几何模型

图5是本发明算例重轨的无网格RKPM离散节点示意图

图6是发明算例重轨的无网格RKPM积分背景网格示意图

图7是本发明算例有限元法计算得到的重轨材料方向角为0°时的热变形总位移场云图

图8是本发明算例无网格RKPM法计算得到的重轨材料方向角为0°时热变形总位移场云图

图9是本发明算例有限元方法计算得到的重轨材料方向角为0°时的Mises应力场云图

图10是本发明算例无网格RKPM法计算得到的重轨材料方向角为0°时的Mises应力场云图

图11是本发明算例有限元方法计算得到的重轨材料方向角为30°时的热变形总位移场云图

图12是本发明算例无网格RKPM法计算得到的重轨材料方向角为30°时的热变形总位移场云图

图13是本发明算例有限元方法计算得到的重轨材料方向角为30°时的Mises应力场云图

图14是本发明算例无网格RKPM法计算得到的重轨材料方向角为30°时的Mises应力场云图

图15是本发明算例有限元方法计算得到的重轨材料方向角为60°时的热变形总位移场云图

图16是本发明算例无网格RKPM法计算得到的重轨材料方向角为60°时的热变形总位移场云图

图17是本发明算例有限元方法计算得到的重轨材料方向角为60°时的Mises应力场云图

图18是本发明算例无网格RKPM法计算得到的重轨材料方向角为60°时的Mises应力场云图

图19是本发明算例重轨材料方向角取为0°、30°、60°时沿路径AB热变形总位移的无网格RKPM解、有限元解和参考解对比曲线图

图20是本发明算例重轨材料方向角取为0°、30°、60°时沿路径AB的Mises应力的无网格RKPM解、有限元解和参考解对比曲线图

图21是本发明算例重轨材料热膨胀正交各向异性因子为0.1、1/3、3、10时沿路径AB的热变形总位移的无网格RKPM解、有限元解和参考解对比曲线图

图22是本发明算例重轨材料热膨胀正交各向异性因子为0.1、1/3、3、10时沿路径AB的Mises应力的无网格RKPM解、有限元解和参考解对比曲线图

具体实施方式

参见图1和图2，基于无网格RKPM的各向异性材料结构二维热变形和热应力的分析方法主要包括如下步骤：

首先，建立热导率、应力分量和应变分量在几何坐标系X—Y和材料坐标系1—2之间的转换关系，沿着几何坐标轴方向的热流密度为：

上式中，λ_ij(i,j＝1,2)是随几何坐标变化的热导率系数，

是变换矩阵，λ₁和λ₂是材料坐标系主轴1，2方向的热导率。定义材料的热导率正交各向异性因子Ht＝λ₁/λ₂。

建立热应力分量与应变分量之间的关系为：

上式中

是坐标变换矩阵，热膨胀系数列向量的具体公式为alpxvec＝[α₁ α₂ 0]^T，α₁，α₂分别为材料坐标系主轴1，2方向的热膨胀系数，定义热膨胀正交各向异性因子Pt＝α₁/α₂，

是正交各向异性弹性矩阵，

定义主次泊松比因子Bt＝ν₁₂/ν₂₁，通过修改热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt和材料方向角θ就可进行不同的各向异性材料结构的无网格RKPM热变形和热应力分析。

其次，参见图2，基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法的具体步骤如下：

(1)对计算模型进行无网格RKPM前处理：根据实际工程中所选计算模型的工况，确定基于无网格RKPM的各向异性材料结构的主热导率λ₁、主弹性模量E₁、主泊松比ν₁₂、主热膨胀系数α₁、剪切模量G₁₂、热导率正交各向异性因子Ht、主次泊松比因子Bt、热膨胀正交各向异性因子Pt和材料方向角θ等材料属性；根据给定的材料属性，计算热导率坐标变换矩阵、正交各向异性弹性矩阵、平面转轴矩阵、正交各向异性总弹性矩阵、热膨胀弹性矩阵和热膨胀系数列向量；导入计算模型的几何模型、无网格RKPM离散节点信息和边界条件；生成无网格RKPM积分背景网格；确定高斯点的坐标、权重系数和雅可比；

(2)基于无网格RKPM分析各向异性材料结构的热变形位移场：任意节点处x的位移场u^h(x)可由其影响域内的节点值

拟合出来

式中，

由此得到域内任一点x处的应变为：

式中，

域内任一点x处的应力为：

式中，

D₁＝[T]^-1[Q]([T]^-1)^T，D₂＝[T]^-1[Q]，

ε₀为温度应变，对于平面问题温度应变为：

ε₀＝[α₁(θ_T-θ₀) α₂(θ_T-θ₀) 0]^T (24)

α₁，α₂分别为材料坐标系主轴1，2方向的热膨胀系数，θ₀为参考温度，θ_T为平面稳态热传导边值问题的解。

对位移边界条件引入罚函数

从而有变分式：

将式(20)～(24)代入式(25)得到无网格RKPM热变形位移场的离散控制方程为：

Ku＝f (26)

当x或y方向有位移约束时，相应的s_x(或s_y)等于1；否则为0。

基于无网格RKPM分析各向异性材料结构的热变形位移场的详细步骤如下：

(2.1)计算单个高斯点影响域内的RKPM力刚度矩阵并组装整体RKPM力刚度矩阵；

(2.2)计算单个高斯点影响域内由于温差引起的RKPM热载荷向量并组装整体RKPM温差载荷列向量，温差为该高斯点的RKPM实际温度与参考温度的差值，高斯点的RKPM实际温度求解具体步骤为：(2.2.1)首先根据输入节点和高斯点坐标求各高斯点与节点之间的距离并从小到大排序，取排序为8～10之间的距离为s[k]，该高斯点的动态影响域半径d_mI＝scale×s[k]，scale是大于1的乘子，取值为1.2～1.5，采用圆形影响域，核函数ω_I(x)＝ω(x-x_I)，在第I个节点x_I处的值最大，且具有紧支性，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mI，将r定义为

其中x_I，y_I为第I个节点x_I在x方向和y方向的节点坐标；其次逐个搜索高斯点影响域内的节点，建立单个高斯点影响域内的RKPM热传导矩阵并组装整体RKPM热传导刚度矩阵；(2.2.2)施加混合传热边界条件，采用罚函数法处理第一类传热边界条件，罚因子取值范围为2×10⁶～2×10⁸；(2.2.3)形成系统总的RKPM热刚度矩阵和RKPM温差载荷列向量；(2.2.4)建立各向异性材料无网格RKPM传热离散控制方程，并求解重轨离散节点的RKPM温度参数值；(2.2.5)采用重构核近似对得到的RKPM温度参数值进行插值，就可获得高斯点的RKPM实际温度值；

(2.3)施加位移边界条件和力载荷，采用罚函数法处理位移边界条件；

(2.4)建立各向异性材料无网格RKPM热应力离散控制方程，并求解计算模型的RKPM热变形位移参数值；

(2.5)采用重构核近似对得到的位移参数值进行插值，就可获得计算模型离散节点的RKPM实际热变形位移值；

(3)求解RKPM热应力值：利用得到的RKPM热变形位移参数值，计算出高斯点的RKPM应力值，进而得到节点的RKPM应力值，其详细步骤如下：

(3.1)根据输入节点和高斯点坐标求各节点与高斯点之间的距离并从小到大排序，取排序为9～15的距离为g[k]，该节点的动态影响域半径d_mx＝scale×g[k]，scale是大于1的乘子，取值为1.2～1.5，核函数ω_I(x)＝ω(x-x_I)，在第I个节点x_I处的值最大，且具有紧支性，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mx，将r定义为

其中gx，gy为高斯点坐标；

(3.2)逐个搜索节点影响域内的高斯点并采用重构核近似对节点影响域内的高斯点的RKPM应力进行插值就可得到节点的RKPM应力值；

(4)对计算结果进行后处理，并将无网格RKPM的计算结果与有限元的计算结果以及参考解比较，其详细步骤如下：

(4.1)读入有限元的计算结果和参考解；

(4.2)无网格RKPM计算域内温度结果值的输出及温度云图显示，无网格RKPM计算域内热变形位移结果值的输出及位移云图显示，无网格RKPM计算域内热应力结果值的输出及热应力云图显示；

(4.3)利用输出的结果绘制无网格RKPM、有限元法和参考解的沿某一指定路径的温度曲线图、热变形位移曲线图和热应力曲线图，对三者进行直观地精度比较，并得出关于基于无网格RKPM的各向异性材料结构热变形和热应力分析的结论。

参见图3，本发明基于无网格RKPM各向异性材料结构二维热变形和热应力分析方法的各个子程序的主要功能介绍如下：

HeavyRail_Main：该程序是本分析方法的主程序；

INdata：该程序的主要功能是完成各种初始条件的读入，包括主热导率λ₁、主弹性模量E₁、主泊松比ν₁₂、主热膨胀系数α₁、剪切模量G₁₂、热导率正交各向异性因子Ht、主次泊松比因子Bt、热膨胀正交各向异性因子Pt和材料方向角θ等材料属性，并导入计算模型的几何模型、无网格RKPM离散节点信息和边界条件；

Gauss_Weight：该程序的主要功能是完成高斯积分点的权重系数的指定；

Gauss_Point:该程序的主要功能是确定无网格RKPM背景积分网格以及高斯积分点的坐标、权重系数和雅可比；

IngluenceDomain：该程序的主要功能是确定高斯点影响域内的无网格RKPM节点；

JudgeDiscontinuous：该程序的主要功能是对高斯点影响域内的无网格RKPM节点进行非连续性的可视性准则判断；

JudgeDomain：该程序的主要功能是定义无网格RKPM计算模型的非连续区域；

ShapeFunction：该程序的主要功能是计算重构核近似的形函数；

DomainWeight：该程序的主要功能是计算圆形影响域内的权重和权重的方导数；

Boundry_Conditon：该程序的主要功能是施加三类传热边界条件；

Result_Output：该程序的主要功能是输出无网格RKPM计算结果的数据；

Post_Processing：该程序的主要功能是对无网格RKPM计算结果进行后处理，绘制热变形以及热应力的曲线图等。

下面是本发明方法应用于工程实践的一个实例：

工程用重轨在制造过程中因材料热处理会产生热变形和热应力，图4所示为重轨的截面示意图，其中管重轨上下表面的温度为900℃，其他位置的对流换热系数为100W/(m²·℃)，环境温度为20℃，参考温度定义为20℃且无内部热源，下表面全约束，材料主弹性模量为70MPa，剪切模量为75MPa，为主泊松比为0.1，主次泊松比因子为1/3，主热膨胀系数为1.15×10^-5，热膨胀正交各向异性因子为1/3，选定材料的导热系数λ₁＝35W/(m·℃)，热导率正交各向异性因子为0.1。材料方向角θ分别取0°、30°、60°，热膨胀正交各向异性因子Pt分别取0.1、1/3、3、10，在图4中，所取对比路径AB两个端点的坐标分别为点A(0mm，38.7mm)，点B(0mm，28.8mm)。本发明针对该实例的具体实施步骤如下：

(A)根据重轨的实际工况，确定无网格RKPM各向异性材料结构的主热导率λ₁、主弹性模量E₁、主泊松比ν₁₂、主热膨胀系数α₁、剪切模量G₁₂、热导率正交各向异性因子Ht、主次泊松比因子Bt、热膨胀正交各向异性因子Pt和材料方向角θ等材料属性；根据给定的材料属性，计算热导率坐标变换矩阵、正交各向异性弹性矩阵、平面转轴矩阵、正交各向异性总弹性矩阵、热膨胀弹性矩阵和热膨胀系数列向量；参见图5和图6，导入计算模型的几何模型、离散节点信息和边界条件；生成无网格RKPM积分背景网格；确定高斯点的坐标、权重系数和雅可比；

(B)逐个搜索各高斯点影响域内的RKPM节点并计算其重构核近似的形函数，建立单个高斯点影响域内由于温差引起的RKPM力载荷向量并组装整体RKPM温差力载荷列向量ff，建立单个高斯点影响域内的RKPM力刚度矩阵并组装整体RKPM力刚度矩阵kf；温差为该高斯点的RKPM实际温度与RKPM参考温度的差值，高斯点的RKPM实际温度求解具体步骤为：(a)根据节点和高斯点坐标求各高斯点与节点之间的距离并从小到大排序，并规定影响域内的节点个数，个数取9，按照距离先后顺序采用可视性准则判断RKPM温度场的非连续性，即节点是否跨越边界，如果跨界则删掉该节点，并补充距离高斯点最近的新的节点进入，再次进行可视性准则判断，直到满足条件的节点达到规定的节点个数，取符合条件的节点中距离高斯点最远的节点和该高斯点之间的距离为s[k]，该高斯点的动态影响域半径大小为d_mI＝scale×s[k]，scale取值为1.2，采用圆形影响域，核函数为ω_I(x)＝ω(x-x_I)，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mI，将r定义为

其中x_I，y_I为第I个节点x_I在x方向和y方向的节点坐标；(b)逐个搜索各高斯点影响域内的节点，建立单个高斯点影响域内的RKPM热传导矩阵并组装整体RKPM热传导刚度矩阵k；(c)通过输入的边界节点信息，求各边界上高斯点的信息，包括边界高斯点的坐标、权重系数和雅可比，遍寻重轨上下边界的高斯点并求其重构核近似形函数，并将各节点形函数与温度T、罚因子、权重系数、雅可比的乘积作为该节点的RKPM热载荷施加附加项，并组装成处理第一类边界条件产生的RKPM热载荷列向量修正项ft1，同时将各个节点之间的形函数乘积与罚因子、权重系数、雅可比相乘以作为处理第一类边界条件产生的RKPM热刚度矩阵修正项kt1，罚因子取为2×10⁸；遍寻重轨剩余边界上的高斯点并求其重构核近似形函数，并将各节点形函数与对流系数h、环境介质温度T_f、权重系数、雅可比的乘积作为该节点的RKPM热载荷施加附加项，并组装成处理第三类边界条件产生的RKPM热载荷列向量修正项fh3，同时将各个节点之间的形函数乘积与对流系数h、权重系数、雅可比相乘以作为处理第三类边界条件产生的RKPM热刚度矩阵修正项kh3；整体RKPM热刚度矩阵为k＝k+kt1+kh3，整体RKPM热载荷列向量为f＝ft1+fh3；(d)建立各向异性材料无网格RKPM传热离散控制方程，并求解重轨离散节点的RKPM温度参数值；(e)逐步搜索计算域内各节点影响域内的节点并求对应节点的重构核近似的形函数，利用各个节点处的RKPM温度参数值求高斯点的RKPM实际温度值；

(C)利用重轨边界上高斯点的信息，遍寻重轨下边界的高斯点并求其重构核近似的形函数，并将各节点形函数与位移向量、罚因子、权重系数、雅可比的乘积作为该节点的RKPM力载荷施加附加项，并组装成处理第一类边界条件产生的RKPM力载荷列向量修正项fp1，同时将各个节点之间的形函数乘积与罚因子、权重系数、雅可比相乘以作为处理第一类边界条件产生的RKPM力刚度矩阵修正项kp1，罚因子取为2×10¹⁴；整体RKPM力刚度矩阵为k2＝kf+kp1，整体RKPM力载荷列向量为f2＝ff+fp1；

(D)建立各向异性材料无网格RKPM热应力离散控制方程，并求解重轨离散节点的RKPM热变形位移参数值；

(E)逐步搜索计算域内各节点影响域内的节点并求对应节点的重构核近似形函数，结合各个节点处的RKPM热变形位移参数值求节点的RKPM实际热变形位移值；

(F)逐个搜索各高斯点影响域内的节点，建立单个高斯点影响域内弹性矩阵和几何矩阵，并利用得到的位移参数值代入应力公式(23)进行计算，就可得到高斯点的RKPM应力值；

(G)根据输入高斯点和节点坐标求各节点与高斯点之间的距离并从小到大排序，并规定影响域内的高斯点个数，个数取为9，按照距离先后顺序采用可视性准则判断RKPM热变形位移场的非连续性，即高斯点是否跨越边界，如果跨界则删掉该高斯点，并补充距离节点最近的新的高斯点进入，再次进行可视性准则判断，直到满足条件的高斯点达到规定的高斯点个数，取符合条件的高斯点中距离节点最远的高斯点和该节点之间的距离为为g[k]，该节点的动态影响域半径d_mx＝scale×g[k]，scale是大于1的乘子，取值为1.5，核函数为ω_I(x)＝ω(x-x_I)，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mx，将r定义为

其中gx，gy为高斯点坐标；

(H)逐个搜索节点影响域内的高斯点并计算其重构核近似的形函数，用影响域内高斯点的形函数乘对应高斯点的RKPM应力值就可得到节点的RKPM应力值；

(I)读入有限元的计算结果和参考解；

(J)无网格RKPM计算域内温度结果值的输出及温度云图显示，无网格RKPM计算域内热变形位移结果值的输出及位移云图显示，无网格RKPM计算域内热应力结果值的输出及热应力云图显示；

(K)利用输出的结果绘制无网格RKPM、有限元法和参考解的沿某一指定路径解的温度曲线图、热变形位移曲线图和热应力曲线图，对三者进行直观的精度比较并得出关于各向异性材料结构热变形和热应力的结论。

图7-图22是本发明算例的后处理结果图，图7是本发明算例有限元法计算得到的重轨材料方向角为0°时的热变形总位移场云图，图8是本发明算例无网格RKPM法计算得到的重轨材料方向角为0°时热变形总位移场云图，图9是本发明算例有限元方法计算得到的重轨材料方向角为0°时的Mises应力场云图，图10是本发明算例无网格RKPM法计算得到的重轨材料方向角为0°时的Mises应力场云图，图11是本发明算例有限元方法计算得到的重轨材料方向角为30°时的热变形总位移场云图，图12是本发明算例无网格RKPM法计算得到的重轨材料方向角为30°时的热变形总位移场云图，图13是本发明算例有限元方法计算得到的重轨材料方向角为30°时的Mises应力场云图，图14是本发明算例无网格RKPM法计算得到的重轨材料方向角为30°时的Mises应力场云图，图15是本发明算例有限元方法计算得到的重轨材料方向角为60°时的热变形总位移场云图，图16是本发明算例无网格RKPM法计算得到的重轨材料方向角为60°时的热变形总位移场云图，图17是本发明算例有限元方法计算得到的重轨材料方向角为60°时的Mises应力场云图，图18是本发明算例无网格RKPM法计算得到的重轨材料方向角为60°时的Mises应力场云图，图19本发明算例重轨材料方向角取为0°、30°、60°时沿路径AB热变形总位移的无网格RKPM解、有限元解和参考解对比曲线图，图20本发明算例重轨材料方向角取为0°、30°、60°时沿路径AB的Mises应力的无网格RKPM解、有限元解和参考解对比曲线图，图21本发明算例重轨材料热膨胀正交各向异性因子为0.1、1/3、3、10时沿路径AB的热变形总位移的无网格RKPM解、有限元解和参考解对比曲线图，图22本发明算例重轨材料热膨胀正交各向异性因子为0.1、1/3、3、10时沿路径AB的Mises应力的无网格RKPM解、有限元解和参考解对比曲线图。从图7-图22中的对比分析可以看出，无网格RKPM法的计算结果与有限元法的计算结果吻合得很好，而且更逼近于参考解，说明本发明所提供的方法计算精度更高，本发明所提供的模型、算法和计算程序完全是有效可行的，本发明可进一步推广到求解其他工程应用中的不同各向异性材料复杂结构的热变形和热应力问题，具有较好的理论研究和工程应用价值。

虽然参考本实施例对本发明进行了详细描述，但以上所述并不限定本发明的保护范围，任何依据本发明思路下的修改和改进等，均视为本发明保护范围。

Claims (3)

1.基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法，其特征在于包括以下步骤：

(1)对计算模型进行无网格RKPM前处理：根据实际工程中所选计算模型的工况，确定无网格RKPM各向异性材料结构的主热导率λ₁、主弹性模量E₁、主泊松比ν₁₂、主热膨胀系数α₁、剪切模量G₁₂、热导率正交各向异性因子Ht、主次泊松比因子Bt、热膨胀正交各向异性因子Pt和材料方向角θ等材料属性；根据给定的材料属性，计算热导率坐标变换矩阵、正交各向异性弹性矩阵、平面转轴矩阵、正交各向异性总弹性矩阵、热膨胀弹性矩阵和热膨胀系数列向量；导入计算模型的几何模型、无网格RKPM离散节点信息和边界条件；生成无网格RKPM积分背景网格；确定高斯点的坐标、权重系数和雅可比；

(2)基于无网格RKPM分析各向异性材料结构的二维热变形位移场，无网格RKPM区别于其他无网格法最大的不同在于其形函数采用重构核近似进行构造，在任意节点x处的热变形位移场u^h(x)可由其影响域内的节点值u_I拟合出来，

式中，C(x；x_I-x)为校正函数，ω(x_I-x)为核函数，Δx_I为第I个节点x_I所对应的面积，φ_I(x)为对应第I个节点x_I的无网格RKPM形函数，写成矩阵形式为：

Φ(x)＝[φ₁(x),φ₂(x)···φ_n(x)]＝p^T(0)M^-1(x)H(x)

其中，

p^T(0)＝[1,0,0,0,0,0]

H(x)＝[p(x₁-x)ω(x₁-x)Δx₁,p(x₂-x)ω(x₂-x)Δx₂,...,p(x_I-x)ω(x_I-x)Δx_I]

p(x_I-x)＝[1,x_I-x,y_I-y,(x_I-x)²,(x_I-x)(y_I-y),(y_I-y)²]^T

基于无网格RKPM的热变形位移场分析的具体步骤为：(A)计算单个高斯点影响域内的RKPM力刚度矩阵并组装整体RKPM力刚度矩阵；(B)计算单个高斯点影响域内由于温差引起的RKPM温差载荷向量并组装整体RKPM温差载荷列向量，温差为该高斯点的RKPM实际温度与RKPM参考温度的差值，高斯点的RKPM实际温度求解的具体步骤为：(a)首先根据节点和高斯点坐标求各高斯点与节点之间的距离并从小到大排序，并规定影响域内的节点个数，一般取8～10个为宜，按照距离先后顺序采用可视性准则判断RKPM温度场的非连续性，即节点是否跨越边界，如果跨界则删掉该节点，并补充距离高斯点最近的新的节点进入，再次进行可视性准则判断，直到满足条件的节点达到规定的节点个数，取符合条件的节点中距离高斯点最远的节点和该高斯点之间的距离为s[k]，该高斯点的动态影响域半径d_mI＝scale×s[k]，scale是大于1的乘子，取值为1.2～1.5，采用圆形影响域，核函数ω_I(x)＝ω(x-x_I)，在第I个节点x_I处的值最大，且具有紧支性，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mI，将r定义为

其中x_I，y_I为第I个节点x_I在x方向和y方向的节点坐标；其次逐个搜索高斯点影响域内的节点并计算其重构核近似形函数及其导数，建立单个高斯点影响域内的RKPM热传导矩阵并组装整体RKPM热传导刚度矩阵；(b)施加混合传热边界条件，采用罚函数法处理第一类传热边界条件，罚因子α取值范围为2×10⁶～2×10⁸；(c)形成系统总的RKPM热刚度矩阵和RKPM温差载荷列向量；(d)建立各向异性材料结构无网格RKPM稳态传热离散控制方程，并求解计算模型离散节点的RKPM温度参数值；(e)采用重构核近似对得到的RKPM温度参数值进行插值，就可获得计算模型高斯点的RKPM实际温度值；(C)施加位移边界条件和力载荷，采用罚函数法处理位移边界条件，罚因子的经验取值为10³×E，其中E为主弹性模量；(D)建立各向异性材料结构无网格RKPM热应力离散控制方程，并求解计算模型离散节点的RKPM热变形位移参数值，具体方程为Ku＝f，其中

其中，D₁＝[T]^-1[Q]([T]^-1)^T，D₁为正交各向异性弹性矩阵，D₂＝[T]^-1[Q]，D₂为热膨胀弹性矩阵，

[T]^-1为平面转轴矩阵，

[Q]为正交各向异性弹性矩阵，其中

Q₃₃＝G₁₂，

φ_i,x和φ_i,y分别为无网格RKPM形函数关于坐标x和y的导数，∫_ΩΦ_if_idΩ为RKPM温差载荷列向量，

为处理位移边界条件引起的RKPM力载荷列向量修正项，α为罚因子，

为外力引起的RKPM力载荷列向量，当x或y方向有位移约束时，相应的s_x(或s_y)等于1；否则为0；θ为材料方向角；E₁和ν₁₂分别为材料沿x方向的弹性模量和泊松比，E₂和ν₂₁分别为材料沿y方向的弹性模量和泊松比，要保证E₁ν₂₁＝E₂ν₁₂；(E)采用重构核近似对得到的RKPM热变形位移参数值进行插值，就可获得计算模型离散节点的RKPM实际热变形位移值；

(3)求解RKPM热应力值：利用步骤(2)得到的RKPM热变形位移参数值，计算出高斯点的RKPM应力值，进而得到计算模型离散节点的RKPM应力值，具体方程为

式中，ε₀为温度应变，对于平面问题温度应变为ε₀＝[α₁(θ_T-θ₀) α₂(θ_T-θ₀) 0]^T，α₁，α₂分别为材料坐标系主轴1，2方向的热膨胀系数，θ₀为参考温度，θ_T为步骤(2)所获得的计算模型高斯点的RKPM实际温度值，即基于无网格RKPM各向异性材料结构二维稳态热传导问题的解；具体步骤为：(A)首先根据输入高斯点和节点坐标求各节点与高斯点之间的距离并从小到大排序，并规定影响域内的高斯点个数，一般取9～12个为宜，按照距离先后顺序采用可视性准则判断RKPM热变形位移场的非连续性，即高斯点是否跨越边界，如果跨界则删掉该高斯点，并补充距离节点最近的新的高斯点进入，再次进行可视性准则判断，直到满足条件的高斯点达到规定的高斯点个数，取符合条件的高斯点中距离节点最远的高斯点和该节点之间的距离为g[k]，该节点的动态影响域半径d_mx＝scale×g[k]，scale是大于1的乘子，取值为1.2～1.5，核函数ω_I(x)＝ω(x-x_I)，在第I个节点x_I处的值最大，且具有紧支性，将高斯点与节点的距离正则化处理，r＝||gx-x_I||/d_mx，将r定义为

其中gx，gy为高斯点坐标；(B)逐个搜索节点影响域内的高斯点，采用重构核近似对节点影响域内的高斯点的RKPM应力进行插值就可得到节点的RKPM应力值；

(4)对计算结果进行后处理，并将无网格RKPM的计算结果与有限元的计算结果以及参考解比较：(A)读入有限元的计算结果和参考解；(B)无网格RKPM计算模型温度结果值的输出及温度云图显示，无网格RKPM计算模型热变形位移结果值的输出及位移云图显示，无网格RKPM计算模型热应力结果值的输出及热应力云图显示；(C)利用输出的结果绘制无网格RKPM、有限元法和参考解的沿某一指定路径解的温度曲线图、热变形位移曲线图和热应力曲线图，并对三者进行直观的精度比较。

2.根据权利要求1所述基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法，其特征在于通过修改热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt和材料方向角θ来控制正交各向异性材料结构的RKPM热变形和RKPM热应力性能，以实施不同各向异性材料结构的热变形和热应力分析，其采用平面转轴矩阵、正交各向异性总弹性矩阵、热膨胀弹性矩阵将各向异性材料坐标系1—2的热应力与应变关系转换为与几何坐标系X—Y相一致的热应力与应变关系，其具体变换为：

上式中，

是坐标变换矩阵，热膨胀系数列向量的具体公式为alpxvec＝[α₁ α₂ 0]^T，α₁，α₂分别为材料坐标系主轴1，2方向的热膨胀系数，定义热膨胀正交各向异性因子Pt＝α₁/α₂，

是正交各向异性弹性矩阵，

定义主次泊松比因子Bt＝ν₁₂/ν₂₁，通过修改热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt和材料方向角θ就可进行基于无网格RKPM的各种不同各向异性材料结构的热变形和热应力分析。

3.根据权利要求1所述基于无网格RKPM的各向异性材料结构二维热变形和热应力分析方法，其特征在于根据后处理结果可得到关于各向异性材料热变形和热应力的结论：无网格RKPM法的计算结果与有限元法的计算结果吻合得很好，且更逼近于参考解，其计算精度更高；不同的热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt和材料方向角θ对无网格RKPM计算模型的热变形位移场和Mises应力场的分布产生了明显影响，热导率正交各向异性因子Ht、热膨胀正交各向异性因子Pt、主次泊松比因子Bt主要影响RKPM热变形位移场和Mises应力场的强度，材料方向角θ主要影响RKPM热变形位移场和Mises应力场的分布方向。

Images