CN115146402A - 高温合金的氧化过程模拟方法及装置 - Google Patents

高温合金的氧化过程模拟方法及装置 Download PDF

Info

Publication number
CN115146402A
CN115146402A CN202210666206.2A CN202210666206A CN115146402A CN 115146402 A CN115146402 A CN 115146402A CN 202210666206 A CN202210666206 A CN 202210666206A CN 115146402 A CN115146402 A CN 115146402A
Authority
CN
China
Prior art keywords
oxidation
equation
concentration
oxidation process
finite element
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210666206.2A
Other languages
English (en)
Inventor
王心美
李磊
岳珠峰
金鹏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202210666206.2A priority Critical patent/CN115146402A/zh
Publication of CN115146402A publication Critical patent/CN115146402A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/10—Analysis or design of chemical reactions, syntheses or processes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Chemical & Material Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Analytical Chemistry (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本公开提供了一种高温合金的氧化过程模拟方法，属于材料数值模拟技术领域。该高温合金的氧化过程模拟方法可以包括如下步骤：基于热力学、经典氧化动力学以及静力学，建立高温合金的氧化力学‑化学的耦合模型。根据耦合模型，确定高温合金的氧化力学‑化学的有限元模型。根据有限元模型模拟高温合金的氧化过程。在高温合金的氧化模拟过程中，无需进行试验，即可获取高温合金氧化过程中的各种数据。且由于该模拟方法是基于高温合金的氧化力学‑化学的有限元模型，从而能够降低科研人员的工作量，提高计算效率，进而能够完整和准确地模拟出高温合金的氧化过程。

Description

高温合金的氧化过程模拟方法及装置
技术领域
本公开涉及材料数值模拟技术领域，尤其涉及一种高温合金的氧化过程模拟方法及装置。
背景技术
高温合金是航空发动机的关键材料，常用于制造涡轮叶片。模拟高温合金的氧化过程，对于航空发动机的涡轮叶片材料的优化设计具有重要意义。
目前，主要通过氧化试验获得高温合金的氧化行为数据，进而分析高温合金的氧化过程。但是对于氧化试验而言，由于在高温环境下高温合金的氧化数据获取极为困难，而高温合金的氧化过程中材料变换包括微观组织演化、氧化层厚度等数据则更加难以获取，进而无法准确和完整地模拟高温合金的氧化过程。
所述背景技术部分公开的上述信息仅用于加强对本公开的背景的理解，因此它可以包括不构成对本领域普通技术人员已知的现有技术的信息。
发明内容
本公开的目的在于提供一种高温合金的氧化过程模拟方法及装置，无需模拟试验，即可完整和准确的模拟高温合金的氧化过程。
为实现上述发明目的，本公开采用如下技术方案：
根据本公开的第一个方面，提供一种高温合金的氧化过程模拟方法，所述氧化过程模拟方法包括：
基于热力学、经典氧化动力学以及静力学，建立高温合金的氧化力学-化学的耦合模型；
根据所述耦合模型，确定高温合金的氧化力学-化学的有限元模型；
根据所述有限元模型模拟所述高温合金的氧化过程。
在本公开的一种示例性实施例中，基于热力学、经典氧化动力学以及静力学，建立高温合金的氧化力学-化学的耦合模型包括：
基于热力学、经典氧化动力学以及静力学，建立第一方程组，所述第一方程组包括位移场的本构方程和浓度场的本构方程；
所述位移场的本构方程为：
Figure BDA0003691697730000021
所述浓度场的本构方程为：
Figure BDA0003691697730000022
其中，σij为应力，Dijk1为刚度系数，εk1为应变张量，下标i、j、k、1分别表示自由指标，下标s为金属离子和氧离子，下标p为氧化物，Δ为梯度算子，ηs为金属离子和氧离子的化学膨胀系数，cs为金属离子和氧离子的浓度，ηp为氧化物的化学膨胀系数，cp为氧化物的浓度，δk1为克罗内克符号，Js为金属离子和氧离子的扩散通道，Ds为金属离子和氧离子的扩散系数，Fs为常数，tr(ε)为应变的迹，
Figure BDA0003691697730000026
为偏导数，ε为应变，X为位移梯度因子，J为离子扩散通道。
在本公开的一种示例性实施例中，所述ηs通过第一预设公式确定，所述第一预设公式为：
Figure BDA0003691697730000023
所述Fs通过第二预设公式确定，所述第二预设公式为：
Figure BDA0003691697730000024
其中，vm为摩尔体积，
Figure BDA0003691697730000025
为金属离子和氧离子的摩尔体积，E为弹性模量，ν为泊松比，R为玻尔兹曼常数，T为温度。
在本公开的一种示例性实施例中，根据所述耦合模型，确定高温合金的氧化力学-化学的有限元模型还包括：
根据所述位移场的本构方程、位移场的控制方程、力的边界条件建立第一基本方程的弱形式；
根据所述浓度场的本构方程、浓度场的控制方程，浓度场的边界条件建立第二基本方程的弱形式；
根据所述第一基本方程的弱形式和所述第二基本方程的弱形式确定所述有限元模型。
在本公开的一种示例性实施例中，所述位移场的控制方程为：
σij,j+fi＝0；
所述力的边界条件为：
σijnj－ti＝0；
所述第一基本方程的弱形式为：
Figure BDA0003691697730000031
其中，v为求解的单元体区域，fi为体力，s为单元体表面，ti为面力，nj为单元体表面的外法线，δ为变分符号，d为微分，u为单元节点位移，t为时间。
在本公开的一种示例性实施例中，所述浓度场的控制方程为：
Figure BDA0003691697730000032
Figure BDA0003691697730000033
所述浓度场的边界条件为：
q＝-n·J；
所述第二基本方程的弱形式为：
Figure BDA0003691697730000034
其中，q为通过表面上的物质通量，RP为氧化反应速率，n为表面的外法线。
在本公开的一种示例性实施例中，根据所述第一基本方程的弱形式和所述第二基本方程的弱形式确定所述有限元模型包括：
根据所述第一基本方程的弱形式和所述第二基本方程的弱形式，建立第二方程组，所述第二方程组为：
Figure BDA0003691697730000035
Figure BDA0003691697730000041
其中，Iu为位移变化引起的不平衡力，
Figure BDA0003691697730000042
为金属离子和氧离子的浓度变化引起的不平衡力，Nu为位移场形函数矩阵，Nc为浓度场形函数矩阵，Bu为位移场应变矩阵，Bc为浓度场应变矩阵，C为浓度，上标T表示对矩阵的转置。
在本公开的一种示例性实施例中，根据所述第一基本方程的弱形式和所述第二基本方程的弱形式确定所述有限元模型还包括：
根据所述第二方程组建立有限元模型，所述有限元模型为：
Figure BDA0003691697730000043
Figure BDA0003691697730000044
其中，
Figure BDA0003691697730000045
为金属离子浓度变化所引起的不平衡力，
Figure BDA0003691697730000046
为氧离子浓度变化所引起的不平衡力，
Figure BDA0003691697730000047
为氧化物浓度的变化所引起的不平衡力，Δu为位移增量，ΔcA为金属离子浓度增量，ΔcO为氧离子浓度增量，ΔcP为氧化物浓度增量，Δt为时间增量，CA为金属离子浓度，CO为氧离子浓度，CP为氧化物浓度，Kuu、
Figure BDA0003691697730000048
Figure BDA0003691697730000049
分别为单元刚度矩阵的子矩阵。
在本公开的一种示例性实施例中，根据所述有限元模型模拟所述高温合金的氧化过程还包括：
根据所述有限元模型建立用户自定义单元子程序；
根据所述用户自定义单元子程序获取分析数据；
根据所述分析数据以模拟所述高温合金的氧化过程。
根据本公开的第二个方面，提供一种高温合金的氧化过程模拟装置，应用于上述的氧化过程模拟方法，其特征在于，所述氧化过程模拟装置包括：
建立模块，用于基于热力学、经典氧化动力学以及静力学，建立高温合金氧化力学-化学的耦合模型；
确定模块，根据所述耦合模型，确定高温合金的氧化力学-化学的有限元模型；
模拟模块，根据所述有限元模型模拟所述高温合金的氧化过程。
本公开实施方式的高温合金的氧化过程模拟方法，在模拟过程中，无需进行试验，即可获取高温合金氧化过程中的各种数据。且由于该模拟方法是基于高温合金的氧化力学-化学的有限元模型，从而能够降低科研人员的工作量，提高计算效率，进而能够完整和准确地模拟出高温合金的氧化过程。
附图说明
为了更清楚地说明本申请实施例或传统技术中的技术方案，下面将对实施例或传统技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本公开的一种实施方式的高温合金的氧化过程模拟方法的流程图。
图2是本公开的一种实施方式的镍基单晶高温合金氧化初始截面形貌的结构图。
图3是本公开的一种实施方式的镍基单晶高温合金在1000个单位时间氧化后的截面形貌图。
图4是本公开的一种实施方式的镍基单晶高温合金在1000个单位时间氧化后的氧离子分布图。
图5是本公开的一种实施方式的镍基单晶高温合金在1000个单位时间氧化后的氧化物分布图。
图6是本公开的一种实施方式的镍基单晶高温合金在1000个单位时间氧化后的应力分布图。
图7是本公开的一种实施方式的镍基单晶高温合金在1000个单位时间氧化后的应变分布图。
图8是本公开的一种实施方式的镍基单晶高温合金在10000个单位时间氧化后的截面形貌图。
图9是本公开的一种实施方式的镍基单晶高温合金在10000个单位时间氧化后的氧离子分布图
图10是本公开的一种实施方式的镍基单晶高温合金在10000个单位时间氧化后的氧化物分布图。
图11是本公开的一种实施方式的镍基单晶高温合金在1000个单位时间氧化后的应力分布图。
图12是本公开的一种实施方式的镍基单晶高温合金在10000个单位时间氧化后的应变分布图。
图13是本公开的一种实施方式的高温合金氧化过程模拟装置的结构示意图。
图中主要元件附图标记说明如下：
1、建立模块；2、确定模块；3、模拟模块。
具体实施方式
现在将参考附图更全面地描述示例实施例。然而，示例实施例能够以多种形式实施，且不应被理解为限于在此阐述的范例；相反，提供这些实施例使得本公开将更加全面和完整，并将示例实施例的构思全面地传达给本领域的技术人员。所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中，提供许多具体细节从而给出对本公开的实施例的充分理解。
所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中，提供许多具体细节从而给出对本公开的实施例的充分理解。然而，本领域技术人员将意识到，可以实践本公开的技术方案而没有所述特定细节中的一个或更多，或者可以采用其它的方法、组元、材料等。在其它情况下，不详细示出或描述公知结构、材料或者操作以避免模糊本公开的主要技术创意。
当某结构在其它结构“上”时，有可能是指某结构一体形成于其它结构上，或指某结构“直接”设置在其它结构上，或指某结构通过另一结构“间接”设置在其它结构上。
用语“一个”、“一”、“所述”用以表示存在一个或多个要素/组成部分/等；用语“包括”和“具有”用以表示开放式的包括在内的意思并且是指除了列出的要素/组成部分/等之外还可存在另外的要素/组成部分/等。用语“第一”和“第二”等仅作为标记使用，不是对其对象的数量限制。
高温合金是航空发动机的关键材料，常用于制造涡轮叶片。模拟高温合金的氧化过程，对于航空发动机的涡轮叶片材料的优化设计具有重要意义。
目前，主要通过氧化试验获得高温合金的氧化行为数据，进而分析高温合金的氧化过程。但是对于氧化试验而言，由于在高温环境下高温合金的氧化数据获取极为困难，而高温合金的氧化过程中材料变换包括微观组织演化、氧化层厚度等数据则更加难以获取，进而无法准确和完整地模拟高温合金的氧化过程。
本公开实施方式中提供一种高温合金的氧化过程模拟方法，如图1所示，该高温合金的氧化过程模拟方法可以包括步骤S100、步骤S110和步骤S120，其中：
步骤S100，基于热力学、经典氧化动力学以及静力学，建立高温合金的氧化力学-化学的耦合模型。
步骤S110，根据耦合模型，确定高温合金的氧化力学-化学的有限元模型。
步骤S120，根据所述有限元模型模拟所述高温合金的氧化过程。
本公开实施方式的高温合金的氧化过程模拟方法，在模拟过程中，无需进行试验，即可获取高温合金氧化过程中的各种数据。且由于该模拟方法是基于高温合金的氧化力学-化学的有限元模型，从而能够降低科研人员的工作量，提高计算效率，进而能够完整和准确地模拟出高温合金的氧化过程。
下面对本公开实施方式提供的高温合金的氧化过程模拟方法的各步骤进行详细说明：
在步骤S100中，基于热力学、经典氧化动力学以及静力学，建立高温合金的氧化力学-化学的耦合模型。
该高温合金可以指以铁、镍、钴为基，能在600℃以上的高温下长期工作的一类金属材料，其可以用于航空发动机。举例而言，该高温合金可以为Ni3Al基合金。步骤S100可以包括：基于热力学、经典氧化动力学以及静力学，建立第一方程组，该第一方程组包括位移场的本构方程和浓度场的本构方程。
其中：
位移场的本构方程为：
Figure BDA0003691697730000081
浓度场的本构方程为：
Figure BDA0003691697730000082
在上述方程组中，σij为应力，Dijk1为刚度系数，εk1为应变张量，下标i、j、k、1分别表示自由指标(张量采用下标记号法，i，j，k，l取值取1，2，3，表示x，y，z3个空间方向)，下标s为金属离子和氧离子，下标p为氧化物，Δ为梯度算子，ηs为金属离子和氧离子的化学膨胀系数，cs为金属离子和氧离子的浓度，ηp为氧化物的化学膨胀系数，cp为氧化物的浓度，δk1为克罗内克符号，Js为金属离子和氧离子的扩散通道，Ds为金属离子和氧离子的扩散系数，Fs为常数，tr(ε)为应变的迹，
Figure BDA0003691697730000083
为偏导数，ε为应变，X为位移梯度因子，J为离子扩散通道。
上述的ηs通过第一预设公式确定，该第一预设公式为：
Figure BDA0003691697730000084
上述的Fs通过第二预设公式确定，所述第二预设公式为：
Figure BDA0003691697730000085
在上述第一预设公式和第二预设公式中，vm为摩尔体积，
Figure BDA0003691697730000086
为金属离子和氧离子的摩尔体积，E为弹性模量，ν为泊松比，R为玻尔兹曼常数，T为温度。
在步骤S110中，根据所述耦合模型，确定高温合金的氧化力学-化学的有限元模型还包括如下步骤：
步骤S200，根据位移场的本构方程、位移场的控制方程、力的边界条件建立第一基本方程的弱形式。
步骤S210，根据浓度场的本构方程、浓度场的控制方程，浓度场的边界条件建立第二基本方程的弱形式。
步骤S220，根据第一基本方程的弱形式和第二基本方程的弱形式确定有限元模型。
可选的，在步骤S200中，位移场的控制方程为：
σij,j+fi＝0；
力的边界条件为：
σijnj－ti＝0；
其中，根据位移场的控制方程和力的边界条件可以确定位移场微分方程的等效积分弱形式，根据位移场微分方程的等效积分弱形式和位移场的本构方程可以确定第一基本方程的弱形式，其详细步骤如下：
对位移场的控制方程和力的边界条件的权函数取基本变量的变分δui以及其边界值(取负值)，则可得到第三预设公式，该第三预设公式为：
∫vδui(σij,j+fi)dv-∫sδui(σij-ti)ds＝0；
根据弹性力学变形协调条件可得第四预设公式，该第四预设公式为：
Figure BDA0003691697730000091
对第四预设公式进行分部积分可以得到第五预设公式，该第五预设公式为：
Figure BDA0003691697730000092
将第五预设公式代入第三预设公式可以得到第六预设公式，该第六预设公式为位移场微分方程的等效积分弱形式，该第六预设公式为：
∫v(-εεijσij+δuifi)dv+∫sδuitids＝0；
该第六预设的公式的矩阵形式为第七预设公式，该第七预设公式为：
∫vσ:δεdv＝∫st·δuds+∫vf·δudv；
将位移场的本构方程代入第七预设公式中可以得到第一基本方程的弱形式，该第一基本方程的弱形式为：
Figure BDA0003691697730000101
在第一基本方程的弱形式中，v为求解的单元体区域，fi为体力，s为单元体表面，ti为面力，nj为单元体表面的外法线，δ为变分符号，d为微分，u为单元节点位移，t为时间。
可选的，在第一基本方式的弱形式中，可以忽略体力或者加速度作用。即f＝0，则可以得到第八预设公式，该第八预设公式为：
Figure BDA0003691697730000102
可选的，在步骤S210中，浓度场的控制方程为：
Figure BDA0003691697730000103
Figure BDA0003691697730000104
浓度场的边界条件为：
q＝-n·J；
其中，根据浓度场的控制方程和浓度场的边界条件可以确定第二基本方程的弱形式，其详细步骤如下：
对于浓度场的控制方程和浓度场的边界条件，其权函数取基本变量的微分δcs以及其边界值(取负值)可以得到第九预设公式，该第九预设公式为：
Figure BDA0003691697730000105
对第九预设公式进行分部积分可以得到第十预设公式，该第十预设公式为：
Figure BDA0003691697730000106
对上述第十预设公式应用散度定理可以得到第十一预设公式，该第十一预设公式为：
Figure BDA0003691697730000111
将浓度场的本构方程代入上述第十一预设公式即可以得到第二基本方程的弱形式，该第二基本方程的弱形式为：
Figure BDA0003691697730000112
其中，q为通过表面上的物质通量，RP为氧化反应速率，n为表面的外法线。
可选的，对于氧化物的浓度，如果仅考虑其氧化反应，根据第二基本方程的弱形式，则可以得到其第十二预设公式，该第十二预设公式为：
Figure BDA0003691697730000113
在步骤S220中，根据第一基本方程的弱形式和第二基本方程的弱形式确定有限元模型还包括：
根据第一基本方程的弱形式和第二基本方程的弱形式，建立第二方程组，该第二方程组为：
Figure BDA0003691697730000114
Figure BDA0003691697730000115
在上述第二方程组中，Iu为位移变化引起的不平衡力，
Figure BDA0003691697730000116
为金属离子和氧离子的浓度变化引起的不平衡力，Nu为位移场形函数矩阵，Nc为浓度场形函数矩阵，Bu为位移场应变矩阵，Bc为浓度场应变矩阵，C为浓度，上标T表示对矩阵的转置。
下面以一种实施例的方式对上述第二方程组的建立进行详细的说明，具体的：
首先对于位移场来说，对虚位移进行插值可以得到第十三预设公式，该十三预设公式为：
Figure BDA0003691697730000117
上述第十三预设公式中的
Figure BDA0003691697730000121
表示单元节点处的位移。
根据第十三预设公式可以确定其相应的虚应变，该虚应变为：
Figure BDA0003691697730000122
此外，trε可以通过第十四预设公式确定，该第十四预设公式为：
Figure BDA0003691697730000123
则根据第十四预设公式可以确定第十五预设公式，该第十五预设公式为：
Figure BDA0003691697730000124
对上述的第一基本方程的弱形式进行插值，则可以得到该第十五预设公式的离散形式：
Figure BDA0003691697730000125
值得注意的是，由于
Figure BDA0003691697730000126
的任意性，根据第十五预设公式可以得到第十六预设公式，该第十六预设公式为：
Figure BDA0003691697730000127
接着对于浓度场来说，对虚浓度插值可以得到第十七预设公式，该第十七预设公式为：
Figure BDA0003691697730000128
根据第十七预设公式可以确定其相应的虚浓度梯度：
Figure BDA0003691697730000129
对第二基本方程的弱形式进行插值，即可以得到其离散形式，该离散形式为：
Figure BDA00036916977300001210
由于
Figure BDA00036916977300001211
的任意性，根据上述的离散形式可以得到第十八预设公式，该第十八预设公式为：
Figure BDA0003691697730000131
由此，上述的第十六预设公式和第十八预设公式可以组成第二方程组，介于此，根据第一基本方程的弱形式和第二基本方程的弱形式建立第二方程组详细过程已推导完成。
可选的，基于第十八预设公式，对于时间微分问题，利用隐式欧拉时间积分法进行求解，则可以确定在t+Δt时刻的第十九预设公式，该第十九预设公式为：
Figure BDA0003691697730000132
在上述第十九预设公式中，cs为当前时刻的浓度，而
Figure BDA0003691697730000133
为前一时刻该处的浓度。
其中，dt可以通过第二十预设公式确定，该第二十预设公式为：
dt＝tn+1-tn；
可选的，由于对于氧化物来说，其微分方程形式简单，且值依赖于该处金属离子和氧离子浓度，此处省去了其推导过程，其积分方程由第二十一预设公式确定，该第二十一预设公式为：
Figure BDA0003691697730000134
可选的，在步骤S220中还可以包括根据所述第二方程组建立有限元模型，该有限元模型为：
Figure BDA0003691697730000135
Figure BDA0003691697730000136
其中，
Figure BDA0003691697730000137
为金属离子浓度变化所引起的不平衡力，
Figure BDA0003691697730000138
为氧离子浓度变化所引起的不平衡力，
Figure BDA0003691697730000139
为氧化物浓度的变化所引起的不平衡力，Δu为位移增量，ΔcA为金属离子浓度增量，ΔcO为氧离子浓度增量，ΔcP为氧化物浓度增量，Δt为时间增量，CA为金属离子浓度，CO为氧离子浓度，CP为氧化物浓度，Kuu、
Figure BDA0003691697730000141
Figure BDA0003691697730000142
分别为单元刚度矩阵的子矩阵。
在上述有限元模型中的，Kuu、
Figure BDA0003691697730000143
Figure BDA0003691697730000144
可以通过第三方程组确定。通过对第十八预设公式、第十九预设公式、第二十一预设公式取关于基本变量的微分，则可以建立第三方程组，该第三方程组为：
Figure BDA0003691697730000145
Figure BDA0003691697730000146
Figure BDA0003691697730000147
Figure BDA0003691697730000148
Figure BDA0003691697730000149
Figure BDA00036916977300001410
Figure BDA00036916977300001411
Figure BDA00036916977300001412
Figure BDA00036916977300001413
Figure BDA00036916977300001414
Figure BDA00036916977300001415
Figure BDA00036916977300001416
Figure BDA00036916977300001417
在步骤S120中，根据有限元模型模拟高温合金的氧化过程还包括如下步骤：
步骤S300，根据有限元模型建立用户自定义单元子程序。
步骤S310，根据用户自定义单元子程序获取分析数据。
步骤S320，根据分析数据以模拟高温合金的氧化过程。
可选的，在步骤S300之前，还可以包括如下步骤：
步骤S400，根据有限元模式确定模型输入文件。
步骤S410，根据模型输入文件确定模型修改文件。
具体的，在步骤S400中，以镍基单晶高温合金为模拟对象，在有限元分析软件中建立模拟模型。该有限元分析软件可以为ABAQUS、CAE(Computer-Aided-Engineering，计算机辅助设计)、ANSYS(analysis-system分析系统)等。该模拟模型的建立依据镍基单晶高温合金的材料特点，充分考虑其两相结构，将二维壳分割成具有γ/γ’两相结构的CAE模型，如图1所示。在有限元软件中，分析步选择温度-位移耦合分析步，其中分析步总时间为氧化模拟的总时间；定义边界条件，其中需要在初始分析步内设置预定义场，分别对前一步所划分的两相结构赋予不同的温度，温度的设置根据材料的属性设置，分别代表γ相和γ’相中的Al含量。划分网格，生成INP(Input，输入模型文件)，最后以浓度场代替ABAQUS中温度-位移耦合分析步中的温度场。
在步骤S410中，以镍基单晶高温合金氧化模拟的需要对步骤S400中的模型输入文件进行修改。并根据氧化模拟模型的需要，按照用户自定义单元子程序的使用规则和耦合模型对输入文件进行修改，包括自定义单元类型、节点数、模拟所需参数、状态变量以及节点自由度编号，其中自由度包括位移和温度。
接着根据相场模型模拟需求对输出模型文件进行修改，获取模型修改文件还包括：建立一层虚拟单元，将计算单元计算的结果映射到所述虚拟单元上，以使后处理可直接在有限元分析软件上进行。其中，虚拟单元和实际的计算单元一一对应，仅单元编号不同。
其中，ABAQUS的用户手册中有关于UEL子程序的使用方法和规则。例如，UMAT子程序在实际中应用比较多，而且用户可以在ABAQUS的property模块告诉软件在计算时要用用户定义的材料属性。但是当用户在使用UEL子程序时并不能在ABAQUS的Mesh模块告诉软件在计算时使用用户定义的单元类型，必须通过修改INP文件，告诉软件计算中使用用户定义的单元类型。这就是在使用UEL子程序时必须修改INP文件的第一个原因，这也是UEL子程序使用的一个重要规则。此外由于使用UEL子程序无法直接进行后处理，本公开实现后处理的方法是将UEL子程序计算的结果映射到一层“虚拟”单元上，所以需要在INP文件中定义一层“虚拟”单元，这就时本公开中需要修改INP文件的第二个原因。
其中，INP文件中用户定义单元的接口如下所示：
*User element，nodes＝4，type＝U1001，properties＝0，coordinates＝2，VARIABLES＝52；
其中，*User element为关键词，表示后续单元为自定义单元，如nodes＝4表示单元为四节点单元，type＝U1001表示单元类型为U1001，properties＝0表示单元参数有0个(材料参数可以在INP文件中定义，也可以直接在程序里定义，本示例直接在程序中定义)，coordinates＝2表示计算模型是2维模型，VARIABLES＝52表示每个单元有52个状态变量。
由于使用UEL子程序时无法直接后处理，所以此方法需要建立一层“虚拟单元”，将计算单元计算的结果映射到“虚拟单元”上，以便后期的后处理可直接在ABAQUS上进行。具体方法是在INP文件生成的单元信息后面，定义“虚拟单元”，虚拟单元为ABAQUS中单元库自带单元，虚拟单元和实际的计算单元一一对应，只有单元编号不同。如下所示：
*Element,type＝U1001
1,1,66,873,92
2,66,67,874,873
3,67,68,875,874
4,68,2,876,875
……
*Element,type＝CPE4
3137,1,66,873,92
3138,66,67,874,873
3139,67,68,875,874
3140,68,2,876,875
……
由于选择了温度-位移耦合分析步，所以模型中要至少有一个单元的单元类型的自由度是同时包含位移和温度的，如CPE4T。具体定义如下：
*Node
999996，0.0，0.0
999997，0.01，0.0
999998，0.01，0.01
999999，0.0，0.01
*Nset，nset＝extraElement
999996，999997，999998，999999
*Element，Type＝CPE4T
999999，999996，999997，999998，999999
为了避免报错，其中节点编号和单元编号必须要和前面的单元区分开，不能重复。
在步骤S300中，根据有限元模型建立用户自定义单元子程序可以通过如下步骤实现：
根据所述有限元模型确定RHS，AMATRX，SVARS三个数组，其中RHS是有限元模型的右端项矩阵，AMATRX是有限元模型的单元刚度矩阵，SVARS是有限元模型的状态变量。
下面以一种实施例的方式给出了UEL子程序的Fortran子程序接口，其中RHS和AMATRX两个数组是必须定义的两个数组，其他数组则是根据实际的需求去定义，比如ENERGY数组，如果模型需要计算能量，那就定义，如果不需要则不需要定义。本公开中用到了SVARS数组来存储和传递各个状态变量，比如位移，应力，浓度等。
下面给了UEL子程序的Fortran子程序接口，其中RHS和AMATRX两个数组是必须定义的两个数组，其他数组则是根据实际的需求去定义，比如ENERGY数组，如果模型需要计算能量，那就定义，如果不需要则不需要定义。本公开中用到了SVARS数组来存储和传递各个状态变量，比如位移，应力，浓度等。
编写UEL子程序，ABAQUS为用户提供了固定格式的子程序接口，用户必须按照接口的格式编写程序，子程序接口如下：
Figure BDA0003691697730000181
其中RHS为有限元模型的右端项矩阵，该右端项矩阵可以记为Re，AMATRX为有限元模型的单元刚度矩阵，该单元刚度矩阵可以记为ke，SVARS为状态变量。需要注意的是，浓度场的边界条件也需要在UEL程序中定义，在本实施例中为：
Figure BDA0003691697730000182
在UEL中，也需要添加UVARM子程序用于计算结果的后处理，ABAQUS提供的FORTRAN接口为：
Figure BDA0003691697730000183
在步骤S310中，根据用户自定义单元子程序获取分析数据通过如下步骤实现：
通过工作模块提交在步骤S410中获得的模型修改文件且对步骤S300中建立的用户自定义单元子程序进行分析，并获取结果文件。根据模型修改文件采用用户自定义单元子程序计算的结果传递给自定义输出变量子程序进行计算，以获取分析数据。
示例的，在JOB模块用第一步修改好的INP文件和第二步编写好的的用户自定义单元子程序提交分析，计算完成后即可后处理。由于修改了INP文件，无法将修改后的文件在ABAQUS的前处理中导入，所以直接在ABAQUS的JOB模块的模型文件选择修改好的INP文件，调用调试好的子程序，提交计算。
在步骤S310中，根据分析数据以模拟高温合金的氧化过程包括通过有限元分析软件的后处理模块对分析数据进行后处理，以得到高温合金的氧化过程。
下面结合一种实施例对本公开的步骤S310做更进一步的详细说明。
以模拟镍基单晶高温合金在10000个单位时间内的氧化过程为例，具体的参数设置依据不得材料而定。如图2所示为镍基单晶高温合金氧化初始截面形貌；如图3、图4、图5所示分别为镍基单晶高温合金在1000个单位时间氧化后的截面形貌、氧离子分布、氧化物分布图；如图6、图7所示分别为镍基单晶高温合金在1000个单位时间氧化后的应力和应变分布图；如图8、图9所示分别为镍基单晶高温合金在10000个单位时间氧化后的截面形貌、氧离子分布、氧化物分布图；如图10、图11所示分别为镍基单晶高温合金在10000个单位时间氧化后的应力和应变分布图。
以上仅展示了3个不同时刻内的一些镍基单晶高温合金的计算数据，可以看出，本公开的高温合金氧化过程的模拟方法可以很好地模拟出镍基单晶高温合金的氧化过程，以及包括镍基单晶高温合金微观组织的演化过程，氧化层生长情况和氧化过程中的应力应变分布情况等。这样可以很好的弥补氧化过程试验中的不足，从而能够为材料的设计提供理论参考，尤其可以为涡轮叶片的材料设计提供理论参考。
本实施方式还提供一种高温合金的氧化过程的模拟装置，应用于上述的高温合金的氧化过程模拟方法，如图13所示，该高温合金氧化过程的模拟装置可以包括建立模块1，确定模块2，模拟模块3。其中：
该建立模块1用于基于热力学、经典氧化动力学以及静力学，建立高温合金氧化力学-化学的耦合模型；该确定模块2用于根据耦合模型，确定高温合金的氧化力学-化学的有限元模型；该模拟模块3，用于根据有限元模型模拟高温合金的氧化过程。
本公开实施方式的高温合金氧化过程的模拟装置，在模拟过程中，无需进行试验，即可获取高温合金氧化过程中的各种数据。且由于该模拟模块3基于模拟方法，而模拟方法是基于高温合金的氧化力学-化学的有限元模型，从而能够降低科研人员的工作量，提高计算效率，进而能够完整和准确地模拟出高温合金的氧化过程。
需要说明的是，尽管在附图中以特定顺序描述了本公开中方法的各个步骤，但是，这并非要求或者暗示必须按照该特定顺序来执行这些步骤，或是必须执行全部所示的步骤才能实现期望的结果。附加的或备选的，可以省略某些步骤，将多个步骤合并为一个步骤执行，以及/或者将一个步骤分解为多个步骤执行等，均应视为本公开的一部分。
应可理解的是，本公开不将其应用限制到本说明书提出的部件的详细结构和布置方式。本公开能够具有其他实施方式，并且能够以多种方式实现并且执行。前述变形形式和修改形式落在本公开的范围内。应可理解的是，本说明书公开和限定的本公开延伸到文中和/或附图中提到或明显的两个或两个以上单独特征的所有可替代组合。所有这些不同的组合构成本公开的多个可替代方面。本说明书的实施方式说明了已知用于实现本公开的最佳方式，并且将使本领域技术人员能够利用本公开。

Claims (10)

1.一种高温合金的氧化过程模拟方法，其特征在于，所述氧化过程模拟方法包括：
基于热力学、经典氧化动力学以及静力学，建立高温合金的氧化力学-化学的耦合模型；
根据所述耦合模型，确定高温合金的氧化力学-化学的有限元模型；
根据所述有限元模型模拟所述高温合金的氧化过程。
2.根据权利要求1所述的氧化过程模拟方法，其特征在于，基于热力学、经典氧化动力学以及静力学，建立高温合金的氧化力学-化学的耦合模型包括：
基于热力学、经典氧化动力学以及静力学，建立第一方程组，所述第一方程组包括位移场的本构方程和浓度场的本构方程；
所述位移场的本构方程为：
Figure FDA0003691697720000011
所述浓度场的本构方程为：
Figure FDA0003691697720000012
其中，σij为应力，Dijk1为刚度系数，εk1为应变张量，下标i、j、k、1分别表示自由指标，下标s为金属离子和氧离子，下标p为氧化物，Δ为梯度算子，ηs为金属离子和氧离子的化学膨胀系数，cs为金属离子和氧离子的浓度，ηp为氧化物的化学膨胀系数，cp为氧化物的浓度，δk1为克罗内克符号，Js为金属离子和氧离子的扩散通道，Ds为金属离子和氧离子的扩散系数，Fs为常数，tr(ε)为应变的迹，
Figure FDA0003691697720000013
为偏导数，ε为应变，X为位移梯度因子，J为离子扩散通道。
3.根据权利要求2所述的氧化过程模拟方法，其特征在于，所述ηs通过第一预设公式确定，所述第一预设公式为：
Figure FDA0003691697720000014
所述Fs通过第二预设公式确定，所述第二预设公式为：
Figure FDA0003691697720000021
其中，vm为摩尔体积，
Figure FDA0003691697720000022
为金属离子和氧离子的摩尔体积，E为弹性模量，ν为泊松比，R为玻尔兹曼常数，T为温度。
4.根据权利要求2所述的氧化过程模拟方法，其特征在于，根据所述耦合模型，确定高温合金的氧化力学-化学的有限元模型还包括：
根据所述位移场的本构方程、位移场的控制方程、力的边界条件建立第一基本方程的弱形式；
根据所述浓度场的本构方程、浓度场的控制方程，浓度场的边界条件建立第二基本方程的弱形式；
根据所述第一基本方程的弱形式和所述第二基本方程的弱形式确定所述有限元模型。
5.根据权利要求4所述的氧化过程模拟方法，其特征在于，所述位移场的控制方程为：
σij,j+fi＝0；
所述力的边界条件为：
σijnj－ti＝0；
所述第一基本方程的弱形式为：
Figure FDA0003691697720000023
其中，v为求解的单元体区域，fi为体力，s为单元体表面，ti为面力，nj为单元体表面的外法线，δ为变分符号，d为微分，u为单元节点位移，t为时间。
6.根据权利要求5所述的氧化过程模拟方法，其特征在于，所述浓度场的控制方程为：
Figure FDA0003691697720000024
Figure FDA0003691697720000025
所述浓度场的边界条件为：
q＝-n·J；
所述第二基本方程的弱形式为：
Figure FDA0003691697720000031
其中，q为通过表面上的物质通量，RP为氧化反应速率，n为表面的外法线。
7.根据权利要求6所述的氧化过程模拟方法，其特征在于，根据所述第一基本方程的弱形式和所述第二基本方程的弱形式确定所述有限元模型包括：
根据所述第一基本方程的弱形式和所述第二基本方程的弱形式，建立第二方程组，所述第二方程组为：
Figure FDA0003691697720000032
Figure FDA0003691697720000033
其中，Iu为位移变化引起的不平衡力，
Figure FDA0003691697720000034
为金属离子和氧离子的浓度变化引起的不平衡力，Nu为位移场形函数矩阵，Nc为浓度场形函数矩阵，Bu为位移场应变矩阵，Bc为浓度场应变矩阵，C为浓度，上标T表示对矩阵的转置。
8.根据权利要求7所述的氧化过程模拟方法，其特征在于，根据所述第一基本方程的弱形式和所述第二基本方程的弱形式确定所述有限元模型还包括：
根据所述第二方程组建立有限元模型，所述有限元模型为：
Figure FDA0003691697720000035
Figure FDA0003691697720000036
其中，
Figure FDA0003691697720000037
为金属离子浓度变化所引起的不平衡力，
Figure FDA0003691697720000038
为氧离子浓度变化所引起的不平衡力，
Figure FDA0003691697720000039
为氧化物浓度的变化所引起的不平衡力，Δu为位移增量，ΔcA为金属离子浓度增量，ΔcO为氧离子浓度增量，ΔcP为氧化物浓度增量，Δt为时间增量，CA为金属离子浓度，CO为氧离子浓度，CP为氧化物浓度，Kuu、
Figure FDA0003691697720000041
Figure FDA0003691697720000042
分别为单元刚度矩阵的子矩阵。
9.根据权利要求1所述的氧化过程模拟方法，其特征在于，根据所述有限元模型模拟所述高温合金的氧化过程还包括：
根据所述有限元模型建立用户自定义单元子程序；
根据所述用户自定义单元子程序获取分析数据；
根据所述分析数据以模拟所述高温合金的氧化过程。
10.一种高温合金的氧化过程模拟装置，应用于如权利要求1～9任意一项所述的氧化过程模拟方法，其特征在于，所述氧化过程模拟装置包括：
建立模块，用于基于热力学、经典氧化动力学以及静力学，建立高温合金氧化力学-化学的耦合模型；
确定模块，根据所述耦合模型，确定高温合金的氧化力学-化学的有限元模型；
模拟模块，根据所述有限元模型模拟所述高温合金的氧化过程。
CN202210666206.2A 2022-06-13 2022-06-13 高温合金的氧化过程模拟方法及装置 Pending CN115146402A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210666206.2A CN115146402A (zh) 2022-06-13 2022-06-13 高温合金的氧化过程模拟方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210666206.2A CN115146402A (zh) 2022-06-13 2022-06-13 高温合金的氧化过程模拟方法及装置

Publications (1)

Publication Number Publication Date
CN115146402A true CN115146402A (zh) 2022-10-04

Family

ID=83407824

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210666206.2A Pending CN115146402A (zh) 2022-06-13 2022-06-13 高温合金的氧化过程模拟方法及装置

Country Status (1)

Country Link
CN (1) CN115146402A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116306046A (zh) * 2023-05-23 2023-06-23 北京云道智造科技有限公司 一种确定燃烧仿真中组分浓度的方法及装置

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116306046A (zh) * 2023-05-23 2023-06-23 北京云道智造科技有限公司 一种确定燃烧仿真中组分浓度的方法及装置
CN116306046B (zh) * 2023-05-23 2023-10-03 北京云道智造科技有限公司 一种确定燃烧仿真中组分浓度的方法及装置

Similar Documents

Publication Publication Date Title
Michopoulos et al. Modeling and simulation of multiphysics systems
Kamakoti et al. Fluid–structure interaction for aeroelastic applications
Jagota et al. Finite element method: an overview
CN104133933B (zh) 一种高超声速飞行器热环境下气动弹性力学特性分析方法
Farhat et al. Robust and provably second‐order explicit–explicit and implicit–explicit staggered time‐integrators for highly non‐linear compressible fluid–structure interaction problems
Palacios et al. Stanford university unstructured (su 2): an open-source integrated computational environment for multi-physics simulation and design
Vidanović et al. Aerodynamic–structural missile fin optimization
Timme et al. Transonic aeroelastic stability analysis using a Kriging-based Schur complement formulation
Slone et al. A finite volume unstructured mesh approach to dynamic fluid–structure interaction: an assessment of the challenge of predicting the onset of flutter
Noor et al. Advances and trends in computational structural mechanics
Shershnev et al. HyCFS, a high-resolution shock capturing code for numerical simulation on hybrid computational clusters
McDaniel et al. HPCMP CREATETM-AV Kestrel new capabilities and future directions
Zhao et al. Transonic wing flutter predictions by a loosely-coupled method
Thomas et al. Review of model order reduction methods and their applications in aeroelasticity loads analysis for design optimization of complex airframes
CN115146402A (zh) 高温合金的氧化过程模拟方法及装置
CN115525980A (zh) 一种再入飞行器气动外形的优化方法和优化装置
Da Ronch et al. Modeling of unsteady aerodynamic loads
Jonsson et al. Development of flutter constraints for high-fidelity aerostructural optimization
Poole et al. Aerofoil design variable extraction for aerodynamic optimization
CN104008234B (zh) 密集模态含阻尼结构模型修正方法
Thomas et al. Static/dynamic correction approach for reduced-order modeling of unsteady aerodynamics
Neumann et al. Coupling strategies for large industrial models
Farhat et al. Recent advances in reduced-order modeling and application to nonlinear computational aeroelasticity
Gupta et al. Systems identification approach for a computational-fluid-dynamics-based aeroelastic analysis
Dasgupta Locking-free compressible quadrilateral finite elements: Poisson’s ratio-dependent vector interpolants

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination