一种适用于汽轮机叶片流热固耦合分析方法与流程

1.本发明属于发电领域，具体涉及一种适用于汽轮机叶片流热固耦合分析方法。


背景技术：

2.汽轮机进口参数的不断提高将导致前端叶片承受的热载荷不断增加，叶片长期运行在高温、高压和高转速环境下时，会诱使材料发生热疲劳和高温蠕变。而在诸如启停机、变负荷等瞬态运行过程中，高温部件内部热应力和热变形会进一步加剧并引起低周疲劳，进而影响高温区域叶片的运行可靠性和寿命。由于该区域的流动特性具有极其的三维特征，压力、温度以及速度场分布呈现极不均匀状态，因此准确分析高压高温环境下汽轮机叶片的流动特性、受力特征及应变情况，对高性能叶片开发设计及安全性分析具有非常重要的意义。
3.目前对高压高温环境下汽轮机叶片的模拟技术通常采用流动、传热及应力场无耦合独立仿真的计算方法，即先通过cfd软件进行流动计算获取气动边界条件后，再将结果单向传递至fem计算软件进行后续的固体域温度场及应力场计算。这样的简化计算方式使得计算模型的边界条件在一定程度上与真实模型不符，模拟结果将产生不确定的插值误差。


技术实现要素：

4.本发明的目的在于克服上述不足，提供一种适用于汽轮机叶片流热固耦合分析方法，可极大程度的提高模型计算方法的真实性及计算结果可靠性。
5.为了达到上述目的，本发明包括以下步骤：
6.s1，根据汽轮机的流通结构，建立叶片几何模型，并依据实际工况建立流体及固体计算域；
7.s2，分别对流体计算域及固体计算域进行网格刨分，生成多块结构化网格；
8.s3，确定流体及固体计算模型的边界条件；
9.s4，在叶片几何模型上进行流体计算域与固体计算域间的热
‑
流耦合计算，对流体域和固体域之间的传热进行分析，流体域和固体域内同步求解并实时数据传递；
10.s5，在流体域中通过流动求解器数值求解雷诺时均navier
‑
stokes方程组，并引入boussinesq湍流模型假设使湍流计算雷诺时均navier
‑
stokes方程组封闭，通过计算获取流体计算域压力和温度及流速；
11.s6，在固体域中只对热传导方程进行求解，在热连接面建立流体域和固体域之间的温度分布曲线，并通过热平衡方程进行反复的迭代求解以保证两侧热通量守恒传递，通过计算获取节点温度；
12.s7，在叶片几何模型上进行流体计算域与固体计算域间的热
‑
固耦合计算，热
‑
固耦合分析以传递并施加在固体表面的流体域压力结果作为边界条件，并读入s4中热
‑
流耦合计算分析得到的节点温度作为热载荷，施加相应的约束条件以及转速，最终对模型进行综合应力分析。
13.s4中，固体域中对热传导方程进行求解，控制体离散热传导方程为：
[0014][0015]
其中，ω为单元控制体体积，为控制体面上法向量，κ为固体导热系数，ρ为固体材料密度，cp为比热。
[0016]
s4中，流体域和固体域之间的热连接面采用完全非匹配热连接技术，由fnmb分解算法得到的单元进行热耦合方程求解得到每个单元独立的连接面温度和热通量，而将两侧单元上的面温度通过加权平均计算获得单元体温度。
[0017]
s4中，热耦合方程如下：
[0018][0019][0020]
其中，单元面完全未连接时w
i
＝1，当单元面部分连接时0＜w
i
＜1，当完全连接时w
i
＝0。
[0021]
s4中，s5中，通过流动求解器数值求解雷诺时均navier
‑
stokes方程组，具体方法如下：
[0022][0023][0024][0025]
其中，i和j指标取值范围为1、2、3，ρ为密度，s
m
为动量方程源项，s
e
为能量方程源项，τ为剪切应力的张量，为雷诺应力项，为湍流通量，为粘性做功项。
[0026]
s4中，引入boussinesq湍流模型假设使湍流计算雷诺时均navier
‑
stokes方程组封闭，具体方法如下：
[0027][0028][0029]
其中，μ
eff
为有效粘性。
[0030]
s7中，应力具有六个分量，其中三个正应力和三个剪切应力，其矩阵形式为：
[0031][0032]
s7中，应力的变化矩阵形式为：
[0033][0034]
s7中，弹性体应变由应力和温度变化两部分引起：
[0035]
{ε}＝{ε
th
}+{ε
el
}
[0036]
其中{ε
el
}为应力引起的应变向量，{ε
t
}为温度导致的应变向量，弹性体应变对于各项同性材料的表达式为：
[0037]
{ε
th
}＝δt[a a a 0 0 0]
t
[0038]
式中：a为材料线膨胀系数，δt＝t
‑
t
ref
为当地温度与参考温度的温度差。
[0039]
与现有技术相比，本发明首先在所建立的模型上进行流体计算域与固体计算域间的热
‑
流耦合计算，为了实现流体域和固体域热通量和温度的守恒传递，二者之间采用热连接面进行数据传递，流体域和固体域内同步求解并实时数据传递。接着在所建立的模型上进行流体计算域与固体计算域间的热
‑
固耦合计算，其中热
‑
固分析以传递并施加在固体表面的cfd压力结果作为边界条件，并读入热分析得到的节点温度作为热载荷，同时施加相应的约束条件以及转速，从而实现了模型的综合应力分析，可极大程度的提高模型计算方法的真实性及计算结果可靠性。
附图说明
[0040]
图1为本发明实施例的计算模型；
[0041]
图2本发明实施例的叶片流体域网格；
[0042]
图3本发明实施例的叶片固体域网格；
[0043]
图4为本发明实施例的叶片温度场计算结果；
[0044]
图5为本发明实施例的叶片应力场计算结果；其中，(a)为叶片离心应力，(b)为叶片弯扭应力，(c)为叶片热应力，(d)为叶片综合等效应力；
[0045]
图6为本发明的流程图。
具体实施方式
[0046]
下面结合附图对本发明做进一步说明。
[0047]
参见图1至图6，本发明包括以下步骤：
[0048]
步骤一，参照汽轮机真实通流结构，根据几何图纸通过三维建模软件建立叶片几何模型，并依据实际工况建立流体及固体计算域；
[0049]
步骤二，采用商用软件分别对流体计算域及固体计算域进行网格刨分，生成多块结构化网格；
[0050]
步骤三，根据物理实际确定流体及固体计算模型的边界条件；
[0051]
步骤四，在所建立的模型上进行流体计算域与固体计算域间的热
‑
流耦合计算，对流体域和体域之间的传热进行分析，为了实现流体域和固体域热通量和温度的守恒传递，二者之间采用热连接面进行数据传递，流体域和固体域内同步求解并实时数据传递。
[0052]
固体传热分析计算中采用的格子中心有限体积法。由于固体内部没有流动现象，其传热方程可以简化为下式，用于求解基于流固耦合面热通量守恒而引起的给定时间步长内固体温度变化：
[0053][0054]
式中：ω为单元控制体体积/m3；为控制体面上法向量；κ为固体导热系数/w
·
(m
·
k)
‑1；ρ为固体材料密度/kg
·
m
‑3；c
p
为比热/j(kg
·
℃)
‑1。
[0055]
热连接面采用完全非匹配热连接技术，由fnmb分解算法得到的单元进行热耦合方程求解可以得到每个单元独立的连接面温度和热通量，而将两侧单元上的面温度通过加权平均计算可以进一步获得单元体温度。
[0056][0057]
式中当单元面完全未连接时w
i
＝1，当单元面部分连接时0＜w
i
＜1，当完全连接时w
i
＝0。
[0058]
步骤五，在流体域中通过流动求解器数值求解雷诺时均navier
‑
stokes方程组，并引入boussinesq湍流模型假设使湍流计算雷诺时均navier
‑
stokes方程组封闭，通过计算获取流体计算域压力、温度及流速等重要气动参数；
[0059]
通过流动求解器数值求解雷诺时均navier
‑
stokes方程组，其具体形式如下：
[0060][0061][0062][0063]
式中：i和j指标取值范围为(1，2，3)ρ为密度/kg
·
m3；s
m
为动量方程源项；s
e
为能量方程源项；τ为剪切应力的张量；为雷诺应力项；为湍流通量；为粘性做功项。
[0064]
并引入boussinesq湍流模型假设，使湍流计算雷诺时均navier
‑
stokes方程组封闭，其具体形式如下：
[0065][0066][0067]
式中：μ
eff
为有效粘性
[0068]
步骤六，在固体域中只对热传导方程进行求解，在热连接面建立流体域和固体域之间的温度分布曲线，并通过热平衡方程进行反复的迭代求解以保证两侧热通量守恒传递，通过计算获取节点温度等边界参数。热传导方程形式如下：
[0069]
h
f
(t
w
‑
t
f
)＝
‑
h
s
(t
w
‑
t
s
)
[0070]
式中：h
f
为流体侧导热系数/w
·
(m
·
k)
‑1；h
s
为固体侧导热系数/w
·
(m
·
k)
‑1；t
w
为连接面温度/k；t
f
为流体侧第一层网格中心温度/k；t
s
为固体侧第一层网格中心温度/k。
[0071]
固体侧传热系数表达式为：
[0072][0073]
式中：κ为当地热导率/w
·
(m
·
k)
‑1；δy为第一层网格中心到壁面距离/m。
[0074]
而当第一层网格没有在粘性底层内部时需要采用对流热通量，此时流体域和固体域之间的热量传递适合采用结合壁面函数的高雷诺数湍流模型，其传热系数采用下式求解：
[0075][0076]
在粘性底层中(y
+
＜13.2)：
[0077]
t
+
＝pr y
+
[0078]
在湍流层中：
[0079][0080]
步骤七，在所建立的模型上进行流体计算域与固体计算域间的热
‑
固耦合计算，热
‑
固耦合分析以传递并施加在固体表面的流体域压力结果作为边界条件，并读入步骤四热
‑
流耦合计算分析得到的节点温度作为热载荷，进一步施加相应的约束条件以及转速，最终对模型进行综合应力分析。
[0081]
三维问题中弹性体内任意点的应力具有6个分量：3个正应力和3个剪切应力，其矩阵形式为：
[0082]
[0083]
应变的矩阵形式为：
[0084][0085]
根据热弹性理论，弹性体应变由应力和温度变化两部分引起：
[0086]
{ε}＝{ε
th
}+{ε
el
}
[0087]
其中{ε
el
}为应力引起的应变向量，{ε
th
}为温度导致的应变向量，其对于各项同性材料的表达式为：
[0088]
{ε
th
}＝δt[a a a 0 0 0]
t
[0089]
式中：a为材料线膨胀系数/m
·
k
‑1；δt＝t
‑
t
ref
为当地温度与参考温度的温度差。