一种以定转子最小对称单元为几何模型的电机有限元热分析方法

摘要

本申请公开了的一种以定转子最小对称单元为几何模型的电机有限元热分析方法，本方法通过在定、转子最小对称单元气隙内新增气隙节点，并基于能量守恒原理，通过非对称的气隙边界对流刚度矩阵将定、转子最小对称单元气隙内边界与新增的气隙节点联系起来进行温度场有限元计算；使得定、转子均可取最小对称单元进行计算，既实现了定转子真实几何模型的建模，又最大程度的降低了剖分网格数，简化了计算，提高了计算速度。例如对于12/10极永磁磁通切换电机采用普通的有限元方法进行热分析时至少要取1/2的转子和1/2的定子进行建模分析，而采用本方法则只需要取1/12的定子和1/10的转子进行分析就可以了。

说明书

技术领域

本发明涉及有限元热分析领域，尤其涉及一种应用于电机有限元热分析中处理定转 子对称性不同的方法。

背景技术

随着对于能源危机的加剧，社会对高效率，低成本的电机需求持续增加。而准确的 电机设计不仅需要进行电机电磁性能的计算，还需要准确计算电机工作时的温升以保证 其安全运行。传统的电机温升计算主要是基于热路法，这种方法存在通常只能计算绕组 的平均温升，且热路模型的建立和参数的校准需要结合实验。而随着有限元方法的出现 和计算机技术的飞速发展，许多成熟的商用有限元软件平台不断发展壮大。CAE (computeraidedengineering，计算机辅助工程)技术已经成为不可或缺的工具。

但采用现有有限元热分析模型进行温度场建模时不能很好的解决定转子对称性不 同的问题，在进行有限元建模时要求定转子必须取相同的比例的定转子单元进行分析， 由于温度场的3维特性往往需要较大规模的剖分网格数来进行计算。

发明内容

发明目的：针对上述现有技术，提出一种以定转子最小对称单元为几何模型的电机 有限元热分析方法，最大程度上减少网格剖分数量，提高了计算速度。

技术方案：一种以定转子最小对称单元为几何模型的电机有限元热分析方法，包括 如下具体步骤：

1)，以定、转子最小对称单元为几何模型，通过网格剖分程序对几何模型进行网格 剖分得到有限元计算所需的网格和节点信息；

2)，在节点信息中新增一个表示定、转子之间气隙温度的气隙节点，并将该节点与 定、转子气隙边界面上的网格单元关联起来，得到稳态导热微分方程；其中，所述稳态 导热微分方程包含定、转子之间气隙边界与气隙节点温度之间满足的对流边界条件；

3)，根据所述稳态导热微分方程以及定、转子气隙边界流出热量守恒定律得到非对 称的定、转子气隙边界对流单元刚度矩阵；其中，所述非对称的定、转子气隙边界对流 单元刚度矩阵包括定子气隙边界的对流刚度矩阵和转子气隙边界的对流刚度矩阵；

4)，将所述非对称的定、转子气隙边界对流单元刚度矩阵叠加到整体刚度矩阵中， 利用非对称求解器求解整个几何模型区域上的温度分布。

进一步的，所述步骤2)中，定、转子之间气隙边界与气隙节点温度之间满足的对 流边界条件为：

$-k\frac{\partial T}{\partial n}{|}_{{\Gamma }_{gap}}(T-T_{gap})$

其中，k为导热系数，h为定、转子与气隙接触的边界面对流散热系数，T为温度， n为外边界的单位法向量，Γ_gap为定、转子与气隙接触的面边界，T_gap为气隙节点温度；

得到所述稳态导热微分方程的弱解积分形式为：

$\left(\begin{array}{c}{\int }_{\Omega }(k_x·\frac{\partial T}{\partial x}\frac{\partial \delta T}{\partial x}+k_y·\frac{\partial T}{\partial y}\frac{\partial \delta T}{\partial y}+k_z·\frac{\partial T}{\partial z}\frac{\partial \delta T}{\partial z})·d\Omega +{\int }_{{\Gamma }_a}h·T·\delta T·d\Gamma +{\int }_{{\Gamma }_{gap}}h(T-T_{gap})·\delta T·d\Gamma \\ ={\int }_{{\Gamma }_a}h·T_a·\delta T·d\Gamma +{\int }_{\Omega }q·\delta T·d\Omega \end{array}\right)$

其中，Ω为求解区域，对应为进行网格剖分的几何模型；(x,y,z)为直角坐标系下 坐标，k_x、k_y、k_z分别为导热系数k在x、y、z轴方向上的分量，δT为虚位移，Γ_a为几 何模型除Γ_gap以外的对流边界，Γ为边界变量，T_a为对流边界Γ_a对应的环境温度，q为 热源密度。

进一步的，所述步骤3)中，所述定、转子气隙边界流出热量守恒定律为：从定子 侧气隙边界流出的热量与从转子侧气隙边界流出的热量之和为0的能量守恒关系，即：

$\frac{360}{{\theta }_2}·{\int }_{{\Gamma }_{stator}}h(T-T_{gap})·d\Gamma +\frac{360}{{\theta }_1}·{\int }_{{\Gamma }_{rotor}}h(T-T_{gap})·d\Gamma =0$

其中，θ₁为转子的旋转对称角度，θ₂为定子的旋转对称角度，Γ_stator为定子侧气隙 边界，Γ_rotor为转子侧气隙边界；

采用有限元法对求解区域进行离散剖分时，若定、转子与气隙接触的边界上的一个 单元e上的温度T^(e)表示为：

$T^{\left(e\right)}=\underset{j=1}{\overset{m}{\Sigma }}{T}_j·N_j$

其中，m为单元e包含的节点总数，T_j为节点j的温度，N_j为节点j对应的单元插 值基函数；则所述定子气隙边界的对流刚度矩阵[N]_s为：

${\left[N\right]}_s=\left(\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& ...& a_{1,m}& b_{1,(m+1)}\\ a_{2,1}& a_{2,2}& ...& a_{2,m}& b_{2,(m+1)}\\ ...& ...& ...& ...& ...\\ a_{m,1}& a_{m,2}& ...& a_{m,m}& b_{m,(m+1)}\\ b_{(m+1),1}& b_{(m+1),2}& ...& b_{(m+1),m}& c_{(m+1),(m+1)}\end{array}\right)$

所述转子气隙边界的对流刚度矩阵[N]_r为：

${\left[N\right]}_r=\left(\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& ...& a_{1,m}& b_{1,(m+1)}\\ a_{2,1}& a_{2,2}& ...& a_{2,m}& b_{2,(m+1)}\\ ...& ...& ...& ...& ...\\ a_{m,1}& a_{m,2}& ...& a_{m,m}& b_{m,(m+1)}\\ \frac{{\theta }_2}{{\theta }_1}{b}_{(m+1),1}& \frac{{\theta }_2}{{\theta }_1}{b}_{(m+1),2}& ...& \frac{{\theta }_2}{{\theta }_1}{b}_{(m+1),m}& \frac{{\theta }_2}{{\theta }_1}{c}_{(m+1),(m+1)}\end{array}\right)$

其中，m+1对应于新增气隙节点，矩阵元素a_i,j、b_i,j和c_i,j为：

$a_{i,j}={\int }_{{\Gamma }_e}h·N_i·N_j·d\Gamma ,i,j=\mathrm{1.2.3}...m$

$b_{i,(m+1)}=b_{(m+1),i}=-\underset{{\Gamma }_e}{\int }h·N_i·d\Gamma ,i=1,2,3...m$

c_(m+1),(m+1)＝h·S_e

其中，S_e为单元e的面积，N_i为节点i对应的单元插值基函数，N_j为节点j对应的 单元插值基函数，Γ_e为单元e所在区域。。

有益效果：现存的一些电机，如永磁磁通切换电机，其定转子对称性不同。如附图 1所示的永磁磁通切换电机定子具有旋转θ₂的中心对称性，即最小可取定子的θ₂/360° 部分进行分析；而转子则具有旋转θ₁的中心对称性，即最小对称单元是转子的θ₁/360°。 在有限元热分析模型中，为了减少网格的数目，希望尽可能的取最小单元进行分析。但 在现有的有限元方法中要求定转子必须取相同比例的定转子单元进行分析，这意味着至 少要取θ₁、θ₂的最小公倍数除以360°的定、转子部分进行分析。本申请的方法通过在定、 转子最小对称单元气隙内新增气隙节点，并基于能量守恒原理，通过非对称的气隙边界 对流刚度矩阵将定、转子最小对称单元气隙内边界与新增的气隙节点联系起来进行温度 场有限元计算；使得定、转子均可取最小对称单元进行计算，既实现了定转子真实几何 模型的建模，又最大程度的降低了剖分网格数，简化了计算，提高了计算速度。例如对 于12/10极永磁磁通切换电机采用普通的有限元方法进行热分析时至少要取1/2的转子 和1/2的定子进行建模分析，而采用本文所述的方法则只需要取1/12的定子和1/10的 转子进行分析就可以了。

附图说明

图1为12/10极永磁磁通切换电机的截面示意图；

图2是生成有限元程序所需的pde文件，文件名为ell.pde；

图3是生成有限元程序所需的fbc文件，文件名为ell.fbc；

图4是生成有限元程序所需的mdi文件，文件名为le.mdi；

图5是生成有限元程序所需的gcn文件，文件名为le.gcn；

图6是修改网格几何信息的程序流程图；

图7是alq4.for单元子程序的程序流程图；

图8是采用本发明提出的方法计算得到的一个12/10极永磁磁通切换电机温度分。

具体实施方式

下面结合附图对本发明做更进一步的解释。

本实施例以图1所示的12/10极永磁磁通切换电机为例，对本方法进行进一步说明。 其中θ₁等于36°，θ₂等于30°，转子由硅钢片叠压而成，其在平行于截面方向的导热率 为46W/(m·℃)，垂直于截面方向的导热率为1.0W/(m·℃)。定子上硅钢片的导热率与转 子上的相同，定子侧的永磁体导热率为9W/(m·℃)，绕组在平行于截面方向的导热率为 0.5W/(m·℃)垂直于截面方向的导热率为200W/(m·℃)。在边界11、13、14-17上均设置 为对称边界条件以反应电机的周期对称性。在边界12、18-20上均设置环境温度为7℃、 对流散热系数为15W/(m²·℃)的自然对流散热边界。定子气隙内边界面边界(6-10)和转子 气隙内边界面边界(1-5)与新增气隙节点关联起来，定、转子与气隙接触的内边界面对流 散热系数为113W/(m²·℃)。另外，在三维模型轴向的两个端面均设置为对称边界条件。

以定转子最小对称单元为几何模型的电机有限元热分析方法，包括如下具体步骤：

步骤1)，以定、转子最小对称单元为几何模型，通过网格剖分程序对几何模型进行 网格剖分得到有限元计算所需的网格和节点信息。本实施例中，有限元体单元的网格剖分 以最常用的四面体单元来剖分，边界单元为三角形；但不限于该剖分方方法，有限元体单元的 网格剖分还可采用六面体，三棱柱；边界单元也可为四边形或其它单元。

步骤2)，在节点信息中新增一个表示定、转子之间气隙温度的气隙节点，并将该节 点与定、转子气隙边界面上的网格单元关联起来，得到稳态导热微分方程。稳态导热微 分方程可表示为：

$\left(\begin{array}{cc}▿·(k▿T)+q=0& \left(1\right)\\ T{|}_{{\Gamma }_0}=T_0& \left(2\right)\\ -k·\frac{\partial T}{\partial n}{|}_{{\Gamma }_0}=h·(T-T_a)& \left(3\right)\end{array}\right)$

其中，为哈密顿算子，k为导热系数，q为热源密度，T为温度，n为外边界的单 位法向量，h为对流散热系数，T_a为对流边界对应的环境温度。式(2)表示在Г₀边界上 的温度为恒定值T₀为。式(3)表示在Г_a边界上的环境温度为T_a，对流散热系数为h；由 于新增了气隙节点，式(3)的对流边界条件在表示定、转子之间气隙内边界与气隙节点温 度之间时，具体可表示为：

$-k\frac{\partial T}{\partial n}{|}_{{\Gamma }_{gap}}=h(T-T_{gap})---\left(4\right)$

其中，Γ_gap为定、转子与气隙接触的面边界，T_gap为气隙节点温度，此处的h即为 定、转子与气隙接触的边界面对流散热系数。

将式(3)两端同时乘以未知函数的虚位移形式，并在求解域上积分得到控制方程 的弱解积分形式为：

$\begin{array}{c}{\int }_{\Omega }(k_x·\frac{\partial T}{\partial x}\frac{\partial \delta T}{\partial x}+k_y·\frac{\partial T}{\partial y}\frac{\partial \delta T}{\partial y}+k_z·\frac{\partial T}{\partial z}\frac{\partial \delta T}{\partial z})·d\Omega +{\int }_{{\Gamma }_a}h·T·\delta T·d\Gamma +{\int }_{{\Gamma }_{gap}}h(T-T_{gap})·\delta T·d\Gamma \\ ={\int }_{{\Gamma }_a}h·T_a·\delta T·d\Gamma +{\int }_{\Omega }q·\delta T·d\Omega \end{array}---\left(5\right)$

其中，Ω为求解区域，对应为进行网格剖分的几何模型；(x,y,z)为直角坐标系下 坐标，k_x、k_y、k_z分别为导热系数k在x、y、z轴方向上的分量，δT为虚位移，Γ为边 界变量；式(5)中Γ_a为几何模型除Γ_gap以外的对流边界，T_a为对流边界Γ_a对应的环境 温度。

3)，传统电机有限元热分析方法中未增加气隙节点的情况下，其稳态导热微分方程 弱解积分形式为：

${\int }_{\Omega }(k_x·\frac{\partial T}{\partial x}\frac{\partial \delta T}{\partial x}+k_y·\frac{\partial T}{\partial y}\frac{\partial \delta T}{\partial y}+k_z·\frac{\partial T}{\partial z}\frac{\partial \delta T}{\partial z})·d\Omega +{\int }_{{\Gamma }_a}h·T·\delta T·d\Gamma ={\int }_{{\Gamma }_a}h·T_a·\delta T·d\Gamma +{\int }_{\Omega }q·\delta T·d\Omega ---\left(6\right)$

由于(5)相对于(6)对比，由于新增了一个气隙节点温度作为未知量，需要新增一个 方程来联立求解所有未知节点的温度。新增方程根据定、转子气隙边界流出热量守恒定 律得出，即从定子侧气隙边界流出的热量与从转子侧气隙边界流出的热量之和为0的能 量守恒关系，方程具体为：

$\frac{360}{{\theta }_2}·{\int }_{{\Gamma }_{stator}}h(T-T_{gap})·d\Gamma +\frac{360}{{\theta }_1}·{\int }_{{\Gamma }_{rotor}}h(T-T_{gap})·d\Gamma =0---\left(7\right)$

其中，θ₁为转子的旋转对称角度，θ₂为定子的旋转对称角度，Γ_stator为定子侧气隙 边界，Γ_rotor为转子侧气隙边界。

本实施例中，有限元体单元的网格剖分以最常用的四面体单元来剖分，边界单元为三角 形，那么得到(7)具体可表示为：

$\frac{360}{{\theta }_2}·\underset{{\Gamma }_{stator}}{\Sigma }(-\frac{h}{3}·(T_e^1+T_e^2+T_e^3)·s_e+h·s_e·T_{gap})+\frac{360}{{\theta }_1}·\underset{{\Gamma }_{rotor}}{\Sigma }(-\frac{h}{3}·(T_e^1+T_e^2+T_e^3)·s_e+h·s_e·T_{gap})=0---\left(8\right)$

其中，S_e为边界单元e的面积，T_e¹、T_e²、T_e³为单元e的三个节点的温度。

此式(8)可以从物理意义上进行解释：第一项表示从气隙流入定子侧的热量，第二项 表示从从气隙流入转子侧的热量，两者之和为0刚好符合能量守恒定律。另外，由于上 式中含有θ₂和θ₁，我们可以分别取转子和定子的最小对称单元进行建模，减少了网格 剖分的数量。

根据稳态导热微分方程以及定、转子气隙边界流出热量守恒定律得到非对称的定、 转子气隙边界对流单元刚度矩阵。采用有限元法对求解区域进行离散剖分时，若定、转 子与气隙接触的边界上的一个单元e上的温度T^(e)表示为：

$T^{\left(e\right)}=\underset{j=1}{\overset{m}{\Sigma }}{T}_j·N_j---\left(9\right)$

其中，m为单元e包含的节点总数，T_j为节点j的温度，N_j为节点j对应的单元插 值基函数；则所述定子气隙边界的对流刚度矩阵[N]_s为：

${\left[N\right]}_s=\left(\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& ...& a_{1,m}& b_{1,(m+1)}\\ a_{2,1}& a_{2,2}& ...& a_{2,m}& b_{2,(m+1)}\\ ...& ...& ...& ...& ...\\ a_{m,1}& a_{m,2}& ...& a_{m,m}& b_{m,(m+1)}\\ b_{(m+1),1}& b_{(m+1),2}& ...& b_{(m+1),m}& c_{(m+1),(m+1)}\end{array}\right)---\left(10\right)$

所述转子气隙边界的对流刚度矩阵[N]_r为：

${\left[N\right]}_r=\left(\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& ...& a_{1,m}& b_{1,(m+1)}\\ a_{2,1}& a_{2,2}& ...& a_{2,m}& b_{2,(m+1)}\\ ...& ...& ...& ...& ...\\ a_{m,1}& a_{m,2}& ...& a_{m,m}& b_{m,(m+1)}\\ \frac{{\theta }_2}{{\theta }_1}{b}_{(m+1),1}& \frac{{\theta }_2}{{\theta }_1}{b}_{(m+1),2}& ...& \frac{{\theta }_2}{{\theta }_1}{b}_{(m+1),m}& \frac{{\theta }_2}{{\theta }_1}{c}_{(m+1),(m+1)}\end{array}\right)---\left(11\right)$

其中，m+1对应于新增气隙节点，矩阵元素a_i,j、b_i,j和c_i,j为：

$a_{i,j}=\underset{{\Gamma }_e}{\int }h·N_i·N_j·d\Gamma ,i,j=\mathrm{1.2.3}...m---\left(12\right)$

$b_{i,(m+1)}=b_{(m+1),i}=-\underset{{\Gamma }_e}{\int }h·N_i·d\Gamma ,i=1,2,3...m---\left(13\right)$

c_(m+1),(m+1)＝h·S_e(14)

其中，S_e为单元e的面积，N_i为节点i对应的单元插值基函数，N_j为节点j对应的 单元插值基函数，Γ_e为单元e所在区域。

本实施例中，有限元体单元的网格剖分以最常用的四面体单元来剖分，边界单元为三角 形，则得到定子气隙边界的对流刚度矩阵[N]_s为：

${\left[N\right]}_s=\left(\begin{array}{cccc}\frac{{\mathrm{hs}}_e}{6}& \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{12}& -\frac{{\mathrm{hs}}_e}{3}\\ \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{6}& \frac{{\mathrm{hs}}_e}{12}& -\frac{{\mathrm{hs}}_e}{3}\\ \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{6}& -\frac{{\mathrm{hs}}_e}{3}\\ -\frac{{\mathrm{hs}}_e}{3}& -\frac{{\mathrm{hs}}_e}{3}& -\frac{{\mathrm{hs}}_e}{3}& {\mathrm{hs}}_e\end{array}\right)---\left(15\right)$

转子气隙边界的对流刚度矩阵[N]_r为：

${\left[N\right]}_r=\left(\begin{array}{cccc}\frac{{\mathrm{hs}}_e}{6}& \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{12}& -\frac{{\mathrm{hs}}_e}{3}\\ \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{6}& \frac{{\mathrm{hs}}_e}{12}& -\frac{{\mathrm{hs}}_e}{3}\\ \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{12}& \frac{{\mathrm{hs}}_e}{6}& -\frac{{\mathrm{hs}}_e}{3}\\ -\frac{{\mathrm{hs}}_e}{3}·\frac{{\theta }_2}{{\theta }_1}& -\frac{{\mathrm{hs}}_e}{3}·\frac{{\theta }_2}{{\theta }_1}& -\frac{{\mathrm{hs}}_e}{3}·\frac{{\theta }_2}{{\theta }_1}& {\mathrm{hs}}_e·\frac{{\theta }_2}{{\theta }_1}\end{array}\right)---\left(16\right)$

4)，将非对称的定、转子气隙边界对流单元刚度矩阵叠加到整体刚度矩阵中，利用 非对称求解器求解整个几何模型区域上的温度分布。

根据上述方法，利用有限元程序自动生成系统生成求解线性稳态导热微分方程所需 要的基本程序。采用FEPG(FiniteElementProgramGenerator，有限元程序自动生成系统) 自动生成的有限元的程序仅能解决连续介质的稳态导热问题，无法将不连续的定子和转 子联系起来。但它提供了基本的程序框架、计算单元刚度矩阵的子程序和求解线性方程 组的程序。由于这些程序都是开源的，我们可以在其基础上进行修改和重新编译。首先 需要采用FEPG生成三维稳态热传导有限元问题基本程序框架，如图2-4所示为生成程 序所需的4个脚本文件。下面对它们进行详细说明：

图2为pde文件用于自动生成计算单元刚度矩阵和单元荷载向量等单元子程序。disp u表示求解问题为单变量问题，即每个节点的自由度只有一个为u；coorx,y,z表示求解 问题采用笛卡尔坐标系且为3维问题；shap％1％2表示单元类型符和节点数将在别的文 件中给出。gaus％3表示一个方向上高斯积分节点的数目，这里采用默认值。mate行定 义了材料参数，ekx，eky，ekz表示xyz三个方向上的导热率。ec为热熔，在此稳态问 题中用不到，eq为热源密度。dist行对应单元刚度矩阵，相当于式(6)中左端第一行；load 行表示体单元右端载荷项，相当于式(6)中右端第二项。

图3为fbc文件用于自动生成计算边界条件的单元刚度矩阵和载荷向量等子程序。 disp行的意义和pde文件中相同。coorx,y表示采用二维笛卡尔坐标系，这是因为对于 三维热传导问题其面单元为二维的三角形。shap、gaus和mate行的物理意义与pde文 件中相同，另外此处ea表示对流散热系数，eb表示环境温度。dist行对应边界面单元刚 度矩阵，相当于(6)中左端第二项；load行对应于边界面单元的右端载荷项，相当于(6) 中右端第一项。

图4为mdi文件用于描述坐标信息，给出场符、初值个数、自由度个数、自由度名 称等信息，并给出所选用的体单元和边界面单元的信息。其中3dxyzy表示该问题采用 3为笛卡尔坐标系，并且该问题有材料信息。#a01u对应的意义为：#场符，初值个数， 自由度个数，和自由度名字。pdeellw4表示所用pde文件的文件明为ell，采用w4单 元，w4表示剖分采用4面体单元。fbcellt3,q4表示所用fbc文件的文件明为ell，采用 t3和q4两种边界单元。其中t3表示对流边界的剖分采用三角形剖分，而q4则是为了 表示定、转子气隙边界与气隙节点相关联而形成的4节点单元，在主程序中留下调用该 子程序的接口。并不使用自动生成且q4单元子程序，而需要自己根据[N]_s和[N]_r的具体 内容进行修改。

图5为gcn文件用于描述求解场问题的所用的求解器和计算流程，此处为解决定转子对称 性不同的问题需要采用非对称的气隙对流矩阵，导致整体刚度矩阵不对称，因此需要非对称的 求解器。defi信息段的aell&表示求解问题为单场问题，用a表示场名.。STARTnin表示采用不 存单刚的方式进行初始化。SOLVnin表示采用不存单刚的方式求解，求解器采用非对称直接法。

将图2到图5中的四个文件放入一个同一个文件夹中，就可以通过有限元程序自动生成系 统生成基本的程序框架了。下面给出所得到计算所需几何信息文件的操作步骤：

1)按照上面描述的材料信息和边界条件设定相关参数的值；

2)通过前处理软件gid导入相关的几何模型，并且按图1所示为几何模型指定材料信息同 时添加边界条件；

3)采用四面体网格剖分的方式，得到剖分网格的几何信息，并导出得到coor0，id0，elem0 和disp0文件用于描述初始的到几何信息；

4)将剖分几何模型所产生的原始id0,coor0,elem0和disp0文件重命名为idint,coorint, elemint和dispint；

5)编写一个名为initial.exe的可执行文件，用与将idint,coorint,elemint和dispint转化成计 算需要的几何信息文件。initial.exe的内部计算流程如图6所示。

在执行完initial.exe后就完成从网格剖分得到几何信息文件到计算所需几何信息文件的转变。 首先在coor0中增加了一个用于表示气隙温度的节点，由于该节点的坐标被定义为负数所以在 将来后处理显示结果时将在图上显示不出该节点。在id0文件中，新增节点对应的规格书为1， 表示该节点时温度未知时待求解变量。在disp0文件，新增了一个节点初值为0，并将节点总数 加1。在elem0文件中，首先将气隙边界3和4剖分的单元由三角形单元变成了均与气隙节点关 联的四边形单元，并将其编号分别改为1和2。

下面需要修改生成的计算主程序elea.for和四边形单元子程序alq4.for。而自动生成的alt3.for 文件和aew4.for文件则分别对应三角形单元和四面体单元，这部分不需要修改。修改elea.for 文件时只需要将程序第304行的内容callagq4(r,coef,prmt,es,em,ec,ef,ne)修改成call alq4(r,coef,prmt,es,em,ec,ef,ne)即可。alq4.for需要按照[N]_s和[N]_r的内容编写相应的程序，该子 程序的功能就是对每个单元计算具体的单元刚度矩阵和载荷向量等内容。R为单元节点坐标数 组；coef为非线性系数数组；prmt为材料参数数组；es为单元刚度矩阵；em为单元质量矩阵； ec为单元阻尼矩阵；ef为单元载荷向量；ne为单元编号。其中前三个个数组为输入量，后面的 单元矩阵或向量都是通过这个子程序计算得到的输出量。所编制程序流程图如图7所示。

将elea.for连同ew4.foragt3.foralt3.for和alq4.for等单元子程序一起编译，就可以进行求解 得到了如图8所示的计算结果，这证明了所提方法是可以实现的。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员 来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也 应视为本发明的保护范围。