CN110096760A - 一种工件热形变的数值模拟方法 - Google Patents

Abstract

本公开提供了一种工件热形变的数值模拟方法，属于热机械模拟技术领域。该工件热形变的数值模拟方法包括：建立工件的几何模型；网格划分；建立时间相关偏微分方程和力平衡方程；依次进行各计算步，任一计算步包括求解时间相关偏微分方程并根据结果确定力平衡方程在下一计算步中的初始条件，以及求解力平衡方程并根据结果确定时间相关偏微分方程在下一计算步中的初始条件。该工件热形变的数值模拟方法能够对工件热变形进行多物理问题耦合研究，便于提高工件的研发和评估效率。

Description

一种工件热形变的数值模拟方法

技术领域

本公开涉及热机械模拟技术领域，尤其涉及一种工件热形变的数值模拟方法。

背景技术

许多工件在高温条件下工作时，往往产生热形变，对工件的热形变进行热机械模拟是工件性能评估的重要方法。举例而言，航空发动机运转时，叶片在高温气流以及自身离心力的作用下产生热机械形变；叶片的热机械模拟是研发和评估航空发动机性能的重要环节。宏观热流问题由包含运输、扩散和热源项的时间相关偏微分方程描述，有限体积法尤其适合解此类偏微分方程；因此，现有的热流模拟数值方法均采用有限体积法。宏观叶片应力分布通过解力平衡方程获得，力平衡方程本质上为位移的二阶偏微分方程，这类问题的数值方法均采用有限元法。

然而，由于不同尺度和类型的模拟采用不同数值方法，难以在一个数值方法框架内解决全部问题，严重限制了热机械模拟在航空发动机研发和评估中的功能。各数值方法均有明显的优缺点，有限体积法难解非均质本征应变情况下的力平衡方程，而叶片的温度、组织和缺陷导致的本征应变分布往往是非均质的，因此，有限体积法无法作为应力分析的数值方法。有限元法的网格单元间连续性差，然而流体为连续体，因此，有限元法不适用于宏观热流问题。因此，现有的方法不能满足多物理问题的耦合研究。

所述背景技术部分公开的上述信息仅用于加强对本公开的背景的理解，因此它可以包括不构成对本领域普通技术人员已知的现有技术的信息。

发明内容

本公开的目的在于提供一种工件热形变的数值模拟方法，能够对工件热变形进行多物理问题耦合研究，便于提高工件的研发和评估效率。

为实现上述发明目的，本公开采用如下技术方案：

根据本公开的第一个方面，提供一种工件热形变的数值模拟方法，包括：

建模步骤：建立所述工件的几何模型；

网格划分步骤：对所述工件的几何模型进行网格划分，获得各个网格单元；

建立时间相关偏微分方程步骤：建立时间相关偏微分方程，并确定所述时间相关偏微分方程的约束条件条件和边界条件；

建立力平衡方程步骤：建立力平衡方程，并确定所述力平衡方程的约束条件和边界条件；

数值模拟步骤：按照预设顺序逐步进行多个计算步，其中，任一所述计算步包括：

基于当前计算步的时间相关偏微分方程的初始条件，通过求解所述时间相关偏微分方程获得各所述网格单元的单元中心值；

根据各所述网格单元的单元中心值确定下一所述计算步中所述力平衡方程的初始条件；

基于当前计算步的力平衡方程的初始条件，通过求解所述力平衡方程获得各所述网格单元的顶点值；

根据各所述网格单元的顶点值确定下一所述计算步中所述时间相关偏微分方程的初始条件。

在本公开的一种示例性实施例中，所述时间相关偏微分方程的边界条件包括第一个计算步的时间相关偏微分方程的初始条件；所述力平衡方程的边界条件包括第一个计算步的力平衡方程的初始条件。

在本公开的一种示例性实施例中，建立时间相关偏微分方程包括：

建立一般的时间相关偏微分方程，所述一般的时间相关偏微分方程为：

其中，ξ为要求解的场量，v为运输项的速度，Γ为扩散项的系数，S为源项；表示对时间的偏导；表示对空间的偏导。

在本公开的一种示例性实施例中，建立力平衡方程包括：

建立一般的力平衡方程，所述一般的力平衡方程为：

其中，C：表示弹性矩阵，ε^iel为非弹性应变，ε为总应变。

在本公开的一种示例性实施例中，通过求解所述时间相关偏微分方程获得各所述网格单元的单元中心值包括：

通过有限体积的方法求解时间相关偏微分方程，获得各个所述网格单元的单元中心值。

在本公开的一种示例性实施例中，任一所述计算步中，所述力平衡方程的初始条件包括各所述网格单元的积分点值；

根据各所述网格单元的单元中心值确定下一所述计算步中所述力平衡方程的初始条件包括：

将任一所述网格单元在当前计算步的单元中心值，作为该网格单元下一所述计算步中的积分点值。

在本公开的一种示例性实施例中，通过求解所述力平衡方程获得各所述网格单元的顶点值包括：

通过有限元的方法求解力平衡方程，获得各个网格单元的顶点值。

在本公开的一种示例性实施例中，所述网格单元的顶点值为所述网格单元的顶点的应力和应变值。

在本公开的一种示例性实施例中，获得任一所述网格单元的顶点值包括：

获取任一所述网格单元的各个顶点的应力和应变值；

计算该网格单元的各个顶点的应力和应变值的平均值，作为该网格单元的顶点值。

在本公开的一种示例性实施例中，任一所述计算步中，所述时间相关偏微分方程的初始条件包括各所述网格单元的单元中心值的初始值；

根据各所述网格单元的顶点值确定下一所述计算步中所述时间相关偏微分方程的初始条件包括：

根据任一所述网格单元的顶点的应力和应变值，计算该网格单元的单元中心值，并将该单元中心值作为该网格单元在下一所述计算步中的单元中心值的初始值。

本公开提供的工件热形变的数值模拟方法，能够根据任一计算步中的时间相关偏微分方程的计算结果确定下一计算步中力平衡方程的初始条件，根据任一计算步中的力平衡方程的计算结果确定下一计算步中时间相关偏微分方程的初始条件。因此，该工件热形变的数值模拟方法能够实现流体域和固体域计算结果的相互耦合计算，可用于流体、固体力学耦合问题，实现工件多物理问题的耦合研究。

附图说明

通过参照附图详细描述其示例实施方式，本公开的上述和其它特征及优点将变得更加明显。

图1是本公开实施方式的工件热形变的数值模拟方法流程示意图。

图2是本公开实施方式的网格划分示意图。

具体实施方式

现在将参考附图更全面地描述示例实施例。然而，示例实施例能够以多种形式实施，且不应被理解为限于在此阐述的范例；相反，提供这些实施例使得本公开将更加全面和完整，并将示例实施例的构思全面地传达给本领域的技术人员。所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中，提供许多具体细节从而给出对本公开的实施例的充分理解。

本公开提供一种工件热形变的数值模拟方法，用于对工件进行热机械模拟。该工件可以为航空发动机叶片、发动机轮盘或者内燃机叶片等，本公开对此不做特殊的限定。

如图1所示，该工件热形变的数值模拟方法包括：

建模步骤：建立工件的几何模型；

网格划分步骤：对工件的几何模型进行网格划分，获得各个网格单元；

建立时间相关偏微分方程步骤：建立时间相关偏微分方程，并确定时间相关偏微分方程的约束条件条件和边界条件；

建立力平衡方程步骤：建立力平衡方程，并确定力平衡方程的约束条件和边界条件；

数值模拟步骤：按照预设顺序逐步进行多个计算步，其中，任一计算步包括：

基于当前计算步的时间相关偏微分方程的初始条件，通过求解时间相关偏微分方程获得各网格单元的单元中心值；

根据各网格单元的单元中心值确定下一计算步中力平衡方程的初始条件；

基于当前计算步的力平衡方程的初始条件，通过求解力平衡方程获得各网格单元的顶点值；

根据各网格单元的顶点值确定下一计算步中时间相关偏微分方程的初始条件。

本公开提供的工件热形变的数值模拟方法，能够根据任一计算步中的时间相关偏微分方程的计算结果确定下一计算步中力平衡方程的初始条件，根据任一计算步中的力平衡方程的计算结果确定下一计算步中时间相关偏微分方程的初始条件。因此，该工件热形变的数值模拟方法能够实现流体域和固体域计算结果的相互耦合计算，可用于流体、固体力学耦合问题，实现工件多物理问题的耦合研究。

下面，结合公式对本公开提供的工件热形变的数值模拟方法的各个步骤进行解释和说明。

其中，在建模步骤中，工件的几何模型可以为三维几何模型，该三维几何模型可以通过Solidworks、UG、AutoCAD、Maya等CAD(计算机辅助设计)软件建立，也可以通过三维扫描仪等建模设备对工件进行扫描而建立，或者通过其他可行的方式建立，本公开对此不做特殊的限定。

在网格划分步骤中，可以对整个工件的几何模型划分统一的网格，进而获得各个网格单元。该统一的网格划分，即适用于流体域的时间相关偏微分方程的计算和求解，也适用于固体域的力平衡方程的计算和求解。

在一实施方式中，可以按照图2所示方法进行网格划分。如图2所示，该图中包括3*3共9个网格单元，任一网格单元为黑色边框限定的区域。每个网格单元可以具有多个单元边界(黑色边框部分)，在网格单元的中心位置具有单元中心A，在单元边界的中心位置具有边界中心B，在网格单元的顶点位置具有单元顶点C，在单元顶点C与单元中心A之间具有积分点D。

可以采用MathLAB等计算工具对工件的几何模型进行网格划分，也可以采用ANSYS等网格划分软件对工件的几何模型进行网格划分，或者采用其他的方法对工件的几何模型进行网格划分，本公开对此不做特殊的限制。

在建立时间相关偏微分方程步骤中，时间相关偏微分方程的边界条件可以包括第一个计算步的时间相关偏微分方程的初始条件。

所建立的时间相关偏微分方程可以为一般的时间相关偏微分方程(1)：

其中，ξ为要求解的场量，v为运输项的速度，Γ为扩散项的系数，S为源项；表示对时间的偏导；表示对空间的偏导。该时间相关偏微分方程(1)表明，三种可能的原因导致了ξ的演变，其中，运输项表示在外部因素作用下的整体流动，扩散项表示ξ自身梯度作用下的相互流动，源项S表示ξ的直接增加或减小。

在数值模拟步骤中，在任一计算步中，基于当前计算步的时间相关偏微分方程的初始条件，可以通过有限体积的方法求解时间相关偏微分方程，获得各个网格单元的单元中心值，单元中心值为单元中心的场量值。

在求解时间相关偏微分方程时，可以通过有限体积法对一般的时间相关偏微分方程(1)进行离散化，获得离散化的时间相关偏微分方程。举例而言，在一实施方式中，可以通过如下方法对一般的时间相关偏微分方程(1)进行离散化。

该一般的时间相关偏微分方程(1)可以写成积分的形式，获得公式(2)，即

公式(2)可以在局部坐标系中逐渐离散，然后将离散过程推行至全局坐标系。

其中，时间项可离散为：

其中，V_e为网格划分后每个网格单元的体积；ξ_P为第P个网格单元的单元中心值，为待求解的场量值；为第P个网格单元的单元中心值的初始值；ΔV为单个网格单元的体积；Δt为时间步长。

运输项可离散为：

其中，v_P为第P个网格单元的单元中心的速度；ΔA为网格单元的单元边界的面积，n为网格单元的单元边界的外法向量，下标b指示单元边界；A_e为网格单元的表面积。

扩散项可离散为：

其中，Γ_p为第P个网格单元的扩散项系数。

源项可离散为：

其中，S_P为第P个网格单元的源项。

在离散后，一般的时间相关偏微分方程(1)变为线性方程：

其中，局部坐标系为建立几何模型时的坐标系，全局坐标系为网格划分时的坐标系，局部坐标系和全局坐标系之间可以通过坐标系转换矩阵进行转换。因此，可以利用坐标系转换矩阵实现将局部坐标系中的离散过程推行至全局坐标系，具体方法为本领域常规技术手段，本公开对此不做详述。

在线性方程(7)中，既有网格单元的单元中心值，又有网格单元的单元边界上的值(如等)。也就是说，未知数比线性方程的数量更多，线性方程不可解。为了使得线性方程可解，可以用单元中心值近似取代单元边界上的值。(v_P)_b和(Γ_P)_b可近似为：

(v_P)_b＝av_P+(1-a)v_B (9)

(Γ_P)_b＝aΓ_P+(1-a)Γ_B (10)

其中，a是权重系数，为第P个网格单元的相邻网格单元的单元中心值的初始值；v_B为第P个网格单元的相邻网格单元的单元中心的速度；Γ_B为第P个网格单元的相邻网格单元的扩散项系数。

可以近似为：

其中，r_P为第P个网格单元的单元中心的坐标；r_B为第P个网格单元的相邻网格单元的单元中心的坐标。

在对和作近似后，线性方程(7)可以进一步写为离散的时间相关偏微分方程(12)的形式：

其中，n_b为第P个网格单元的单元边界的外法向量；ΔA_b第P个网格单元的单元边界的面积。

结合当前计算步中时间相关偏微分方程的初始条件，通过对离散的时间相关偏微分方程(12)的求解，可以获得当前计算步中各个网格单元的单元中心值。

可以理解的是，上述离散的时间相关偏微分方程(12)为通过一种示例性的离散化方法对一般的时间相关偏微分方程(1)离散化的结果，技术人员可以通过其他方式对一般的时间相关偏微分方程(1)进行离散化并获得其他的离散的时间相关偏微分方程。

在一实施方式中，任一计算步中，力平衡方程的初始条件包括各网格单元的积分点值。

如此，可以将各个网格单元在当前计算步的单元中心值，作为同一网格单元在下一计算步中的积分点值，实现对力平衡方程在下一计算步中的初始条件的更新，将流体域的计算结果耦合至固体域的计算中。

在建立力平衡方程步骤中，力平衡方程的边界条件可以包括第一个计算步的力平衡方程的初始条件。

可以建立一般的力平衡方程(13)：

其中，

其中，C：表示弹性矩阵，ε^iel为非弹性应变，即塑性应变，是已知量，u为位移，是未知量，ε为总应变，ε^el为弹性应变。

在数值模拟步骤中，可以通过有限元的方法求解力平衡方程，获得各个网格单元的顶点值，网格单元的顶点值指的是网格单元的顶点的值。

在一实施方式中，网格单元的顶点值可以包括单元顶点的应力和应变值。

在一实施方式中，任一网格单元具有多个顶点，可以通过求解力平衡方程获得任一网格单元的各个顶点的顶点值，然后通过求平均值的方法，获得同一网格单元的各个顶点的顶点值的平均值，并将该平均值作为该网格单元的顶点值。

在一实施方式中，可以求解力平衡方程获得网格单元的位移，并根据网格单元的位移确定网格单元的应力和应变值，作为网格单元的顶点值。

在另一实施方式中，可以求解力平衡方程获得任一网格单元的各个顶点的位移，并计算出同一网格单元的各个顶点的应力和应变值，计算出各个顶点的应力和应变值作为该网格单元的顶点值。

举例而言，在一实施方式中，可以通过如下方法求解各个网格单元的位移：

(1)选取矩形网格单元将待求区域进行离散，矩形网格单元有四个节点，选取适当的坐标轴建立笛卡尔坐标系，每个节点i有x，y两个方向的节点位移因此该网格单元有8个自由度。网格单元内部任一点的位移都可假设为双线性的，位移模式为：

u_e，1＝a₁+a₂x+a₃y+O₄xy

u_e，2＝a₅+a₆x+a₇y+a₈xy (15)

其中u_e，1代表某点处位移的x分量，u_e，2表示某点处位移的y分量；方程组(15)可以通过式(16)的形式表示。

定义：

则[u_e]＝[M][a]。

其中，在节点处：

定义：

则

因此，进而可知令[M][A]^-1＝[N_e]，则其中，

[N_e]为形状函数矩阵，其建立起单元节点位移与网格单元内任意一点位移的关系。其中，各形状函数为：

i表示网格单元的第i个节点，Δ表示单元面积，a_i、b_i、c_i、d_i为系数，它们都与网格单元的节点坐标有关，则引入节点坐标矩阵[Λ]，

节点按照逆时针顺序排列，a_i、b_i、c_i、d_i分别为[Λ]矩阵第i行各元素对应的代数余子式，由此可得到形状函数矩阵[N_e]。

(2)在有限元中，通过虚位移原理可以证明单元节点力[P]^e与单元节点位移之间的关系，

其中，[K]^e为单元刚度矩阵，其计算公式为：

[B_e]为应变矩阵，是形状函数矩阵经微分算子矩阵作用所得的结果，[B_e]^T是其转置矩阵，[C]为弹性矩阵，表征应力和应变之间的关系；其中，E是材料的弹性常数，v是泊松比。

单元刚度矩阵[K]^e是8×8的矩阵，将[K]^e按照2×2的子矩阵进行划分，得到

根据公式(26)可以得到：

其中，

由于p_i和都有x和y方向两个分量(下标i表示单元局部节点编号，i＝1、2、3或者4)，每个网格单元的刚度矩阵[K]^e算出后，将单元局部节点编号统一改为整体节点编号，在统一的编号体系下按照节点对应的位置进行叠加，即把相同节点处的单元刚度矩阵元素相加，得到整体刚度矩阵[K]，此时， 为所有节点的位移，[P]为整体的节点载荷矩阵。通过已知的边界条件和约束条件可以得到矩阵[P]，将边界条件和约束条件代入上式，可解得所有网格单元的单元节点位移，即得到代入式中可以得到单元位移[u_e]。

在一实施方式中，任一计算步中，时间相关偏微分方程的初始条件包括各网格单元的单元中心值的初始值。由于同一网格单元的单元中心值与单元顶点值之间具有确定的关系，因此在数值模拟步骤中，可以根据各网格单元的顶点值确定下一计算步中时间相关偏微分方程的初始条件。

举例而言，可以根据网格单元的顶点值，例如网格单元的应力和应变值，计算出该网格单元的单元中心值，然后将计算出来的单元中心值作为该网格单元在下一计算步中的单元中心值的初始值。如此，可以实现对下一计算步的时间相关偏微分方程的初始条件的更新，将固体域计算结果耦合至流体域计算结果中。

由于在求解时间相关偏微分方程时，主要借助有限体积方法；在求解力平衡方程时，主要借助有限元的方法。因此，本公开提供的工件热形变的数值模拟方法可以结合有限体积法和有限元法的优点，实现固体域和流体域的耦合计算，能够有效地应用于多物理问题耦合的分析，提高对工件热变形的分析和评估效率，提供工件的研发和评估效率。

可以理解的是，本公开提供的工件热形变的数值模拟方法，还可以与傅里叶变换法等问题相结合，应用于微观尺寸领域，实现微观应力分析、组织和缺陷的演变分析等领域。

可以理解的是，本公开提供的工件热形变的数值模拟方法适用于各种边界条件，例如可以适用于循环边界条件、自由边界条件、Dirichlet边界条件、Neumann边界条件。本领域技术人员可以通过各类可行的边界条件对所建立的时间相关偏微分方程、力平衡方程等进行边界约束。

需要说明的是，尽管在附图中以特定顺序描述了本公开中方法的各个步骤，但是，这并非要求或者暗示必须按照该特定顺序来执行这些步骤，或是必须执行全部所示的步骤才能实现期望的结果。附加的或备选的，可以省略某些步骤，将多个步骤合并为一个步骤执行，以及/或者将一个步骤分解为多个步骤执行等，均应视为本公开的一部分。

应可理解的是，本公开不将其应用限制到本说明书提出的部件的详细结构和布置方式。本公开能够具有其他实施方式，并且能够以多种方式实现并且执行。前述变形形式和修改形式落在本公开的范围内。应可理解的是，本说明书公开和限定的本公开延伸到文中和/或附图中提到或明显的两个或两个以上单独特征的所有可替代组合。所有这些不同的组合构成本公开的多个可替代方面。本说明书所述的实施方式说明了已知用于实现本公开的最佳方式，并且将使本领域技术人员能够利用本公开。

Claims (9)

1.一种工件热形变的数值模拟方法，其特征在于，包括：

建模步骤：建立所述工件的几何模型；

网格划分步骤：对所述工件的几何模型进行网格划分，获得各个网格单元；

建立时间相关偏微分方程步骤：建立时间相关偏微分方程，并确定所述时间相关偏微分方程的约束条件条件和边界条件；

建立力平衡方程步骤：建立力平衡方程，并确定所述力平衡方程的约束条件和边界条件；

数值模拟步骤：按照预设顺序逐步进行多个计算步，其中，任一所述计算步包括：

基于当前计算步的时间相关偏微分方程的初始条件，通过求解所述时间相关偏微分方程获得各所述网格单元的单元中心值；

根据各所述网格单元的单元中心值确定下一所述计算步中所述力平衡方程的初始条件；

基于当前计算步的力平衡方程的初始条件，通过求解所述力平衡方程获得各所述网格单元的顶点值；

根据各所述网格单元的顶点值确定下一所述计算步中所述时间相关偏微分方程的初始条件。

2.根据权利要求1所述的工件热形变的数值模拟方法，其特征在于，所述时间相关偏微分方程的边界条件包括第一个计算步的时间相关偏微分方程的初始条件；所述力平衡方程的边界条件包括第一个计算步的力平衡方程的初始条件。

3.根据权利要求1所述的工件热形变的数值模拟方法，其特征在于，建立时间相关偏微分方程包括：

建立一般的时间相关偏微分方程，所述一般的时间相关偏微分方程为：

其中，ξ为要求解的场量，v为运输项的速度，Γ为扩散项的系数，S为源项；表示对时间的偏导；表示对空间的偏导。

4.根据权利要求1所述的工件热形变的数值模拟方法，其特征在于，建立力平衡方程包括：

建立一般的力平衡方程，所述一般的力平衡方程为：

σ＝C：(ε-ε^ie1)

其中，C：表示弹性矩阵，ε^ie1为非弹性应变，ε为总应变。

5.根据权利要求1所述的工件热形变的数值模拟方法，其特征在于，通过求解所述时间相关偏微分方程获得各所述网格单元的单元中心值包括：

通过有限体积的方法求解时间相关偏微分方程，获得各个所述网格单元的单元中心值。

6.根据权利要求5所述的工件热形变的数值模拟方法，其特征在于，任一所述计算步中，所述力平衡方程的初始条件包括各所述网格单元的积分点值；

根据各所述网格单元的单元中心值确定下一所述计算步中所述力平衡方程的初始条件包括：

将任一所述网格单元在当前计算步的单元中心值，作为该网格单元下一所述计算步中的积分点值。

7.根据权利要求1所述的工件热形变的数值模拟方法，其特征在于，通过求解所述力平衡方程获得各所述网格单元的顶点值包括：

通过有限元的方法求解力平衡方程，获得各个网格单元的顶点值。

8.根据权利要求7所述的工件热形变的数值模拟方法，其特征在于，获得任一所述网格单元的顶点值包括：

获取任一所述网格单元的各个顶点的应力和应变值；

计算该网格单元的各个顶点的应力和应变值的平均值，作为该网格单元的顶点值。

9.根据权利要求8所述的工件热形变的数值模拟方法，其特征在于，任一所述计算步中，所述时间相关偏微分方程的初始条件包括各所述网格单元的单元中心值的初始值；

根据各所述网格单元的顶点值确定下一所述计算步中所述时间相关偏微分方程的初始条件包括：

根据任一所述网格单元的顶点的应力和应变值，计算该网格单元的单元中心值，并将该单元中心值作为该网格单元在下一所述计算步中的单元中心值的初始值。