井-岩结合的EGS热流耦合数值模拟

Numerical Simulation of Thermal Flow Coupling in Well-Rock Combined EGS

摘要

现有研究忽视了井-岩耦合过程及深井筒传热效应，难以准确评估增强型地热系统（EGS）在干热岩中的采热性能.为此，本文基于COMSOL（多物理场仿真软件）构建了井-岩结合的三维热流耦合模型，分析了各因素对EGS采热性能的影响.模拟结果表明，EGS的采热过程主要依赖裂隙和井筒换热，其中裂隙贡献73%，井筒的换热占27%；在一注一采模式下，注采井间距为500 m时采热效果最佳，出口温度可达145.5 ℃；增加注水流量可提升出口温度，但会缩短热储寿命.通过调整井布局、优化流动传热路径发现：一注一采布局寿命长且出口温度高，而一注四采模式能获得最高采热量，是一注一采的3.1倍.

Abstract

Existing studies overlook the well-rock coupling process and deep wellbore heat transfer effects， making it difficult to accurately evaluate the heat extraction performance of enhanced geothermal system （EGS） in hot dry rock. To address this issue， a three-dimensional thermal flow coupling model for well-rock combination based on COMSOL was developed， and the impact of various factors on the heat extraction performance of the EGS was analyzed. The simulation results show that the heat extraction process of EGS primarily relies on heat transfer through fractures and wellbores， with heat exchange through fractures contributing 73% and that through wellbore accounting for 27%. In a single injection and production mode， the best heat extraction occurs at a 500 m spacing between wells， with an outlet temperature reaching 145.5 ℃. Increasing the injection flow rate raises the outlet temperature， but shortens the heat reservoir’s lifespan. Optimizing well layout and flow heat transfer paths reveals that the single injection and production layout has a longer lifespan and higher outlet temperature， while the single injection and four production mode achieves the highest heat extraction， 3.1 times that of the single injection and production.

Graphical abstract

关键词

干热岩 / 增强型地热系统 / 热流耦合 / 采热性能 / 出口温度

Key words

hot dry rock / enhanced geothermal system / thermal flow coupling / heat extraction performance / outlet temperature

鉴于能源需求的持续增长、气候变化与环境问题的严峻性，以及可持续发展的紧迫性，亟须探索替代能源方案，以有效应对当前的能源危机^［1］.干热岩（HDR）作为一种深层地热资源，拥有巨大的热量储备，具有成为可持续发展关键替代能源的潜力.据统计，即便仅开采干热岩资源中的2%，也足以满足我国2024年全年能源消费量的2 993倍.增强型地热系统（EGS）是目前广泛应用于干热岩热提取领域的重要手段，其中，人工制造的裂缝性储层是热岩与冷工质之间天然的“热交换器”，而注采井是“热交换器”的入口和出口^［2］.

学者们通过调用地质勘探的实测数据模拟得到了不同地质背景下的EGS采热性能，且在干热岩开采过程中的多物理场耦合方面作了大量研究^［3-5］.Chen等^［6］研究发现EGS地下区域的流动传热（TH）在抽热过程中发挥着关键作用，并采用局部非平衡热理论研究地下3 800 m深处500 m范围热储在4种井布置下的采热性能；Huang等^［7］研究了地下4 km深处的500 m范围干热岩热储内的热流耦合情况，并基于一注一采布置下的储层孔隙度分布对EGS采热情况进行了模拟，结果显示该系统在10年后出口温度仍可达173 ℃；Lu等^［8］将注采井等效为出入口边界条件，研究了井眼布置和注采方案对系统热流耦合换热性能的影响，结果表明一注两采三角布置最好，20年后取热量仍能达到1.3 MW/a；Yao等^［9］考虑了岩石基质和裂缝中的变形、渗流和传热过程及其相互作用，建立了热-水-力耦合的三维随机裂缝模型，表明出口温度受进口温度及注采井间距影响较大；单丹丹等^［10］开展了二维井筒-随机裂隙网络储层热流耦合数值模拟研究，分析了井开孔段长度及井筒壁热传导系数对EGS产能与寿命的影响.然而，上述研究都忽略了多物理场耦合过程中的井-岩相互作用；此外，研究对象仅针对深层干热岩小范围热储内的高热区域，这忽略了长深井筒的热损失，对于准确评估实际EGS工程的采热性能存在不足^［11-12］，目前关于井-岩热流耦合模型的井布置方案对HDR产热影响鲜有报道.

本文探讨了井-岩热流耦合作用下多井布置EGS的潜在可行性，采用COMSOL软件建立了考虑长深井筒与围岩及裂隙储层耦合流动传热的三维深层瞬态模型，开展了不同井间距、注水流量下EGS采热效果的研究，并对比分析了单井及多井布置时EGS的采热能力.本文旨在为EGS实际工程提供参考及指导.

1 井-岩流动传热耦合模型建立

1.1 几何模型

选取地下4 500 m深处的辽河盆地沉积型干热岩地热资源作为本文背景，地温梯度35.1 ℃/km.考虑辽宁省实际地质情况及计算效率，建立700 m×700 m×4 500 m模型，热储区域位于地下3 800 m处，热储尺寸为500 m×500 m×500 m.以流动传热理论为基础建立井-岩热流耦合模型，其概念模型如图1所示.由于实际地质情况过于复杂，因此对模型进行适当简化，以节省计算成本：

1） 岩石储层被认为是经过压裂等手段改造过的储层区域，岩层为各向同性连续多孔介质，热物理性质恒定；

2） 水在热储内流动为单向水流，采用连续介质下的双孔隙多孔介质模型，其中基岩及热储区域为超低渗透率、超低孔隙度的多孔介质区域，裂隙层区域的渗透率远高于围岩区域；

3） 不考虑流体损失及流体相态变化；

4） 忽略传热过程热储内水与岩石之间的化学反应和热储内应力作用；

5） 仅考虑热对流及热传导，不考虑辐射作用.

1.2 流动传热控制方程

本文公式中字母的含义如表1所示.

岩石基质内渗流场控制方程：

\frac{\partial }{\partial t}\left({\epsilon }_{\mathrm{s}}{\rho }_{\mathrm{w}}\right)+\nabla \left({\rho }_{\mathrm{w}}{\mathit{u}}_{\mathbf{s}}\right)=\mathrm{0}.(1)

其内流动满足达西定律（考虑重力作用）：

{\mathit{u}}_{\mathrm{s}}=-\frac{K_{\mathrm{s}}}{{\mu }_{\mathrm{w}}}\left(\nabla p-{\rho }_{\mathrm{w}}g\right).(2)

式中：\nabla p为压力的梯度.

岩石基质内多孔介质温度场控制方程：

{\left(\rho c_p\right)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\frac{\partial \theta }{\partial t}+\left({\rho }_{\mathrm{w}}{c}_{\mathrm{w}}\right){\mathit{u}}_{\mathrm{s}}\nabla {\theta }_{\mathrm{s}}+\nabla \left(-{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla {\theta }_{\mathrm{s}}\right)=Q_{\mathrm{T}},(3)

(\rho c_p{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}={\epsilon }_{\mathrm{s}}{\rho }_{\mathrm{w}}{c}_{\mathrm{w}}+\left(\mathrm{1}-{\epsilon }_{\mathrm{s}}\right){\rho }_{\mathrm{s}}{c}_{\mathrm{s}},(4)

{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}={\epsilon }_{\mathrm{s}}{\lambda }_{\mathrm{w}}+\left(\mathrm{1}-{\epsilon }_{\mathrm{s}}\right){\lambda }_{\mathrm{s}}.(5)

式中，Q_T为热源项.

热储裂隙内渗流场及温度场控制方程：

d_{\mathrm{f}}\frac{\partial }{\partial t}\left({\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{w}}\right)+\nabla \left(d_{\mathrm{f}}{\rho }_{\mathrm{w}}{\mathit{u}}_{\mathrm{f}}\right)=\mathrm{0},(6)

{\mathit{u}}_{\mathrm{f}}=-\frac{K_{\mathrm{f}}}{{\mu }_{\mathrm{w}}}\left(\nabla p_{\mathrm{f}}-{\rho }_{\mathrm{w}}g\right),(7)

\begin{array}{c}{d}_{\mathrm{f}}{\left(\rho c_p\right)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\frac{\partial \theta }{\partial t}+d_{\mathrm{f}}\left({\rho }_{\mathrm{w}}{c}_{\mathrm{w}}\right){\mathit{u}}_{\mathrm{f}}\nabla {\theta }_{\mathrm{f}}+d_{\mathrm{f}}\nabla \left(-{\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla {\theta }_{\mathrm{f}}\right)=\\ -d_{\mathrm{f}}{Q}_{\mathrm{T}}.\end{array}(8)

式中，\nabla p_{\mathrm{f}}为裂隙压力梯度.

将井筒简化为一维线性单元，质量方程和动量方程如下：

\frac{\partial \left(A_{\mathrm{c}}{\rho }_{\mathrm{w}}\right)}{\partial t}+\nabla \left(A_{\mathrm{c}}{\rho }_{\mathrm{w}}{u}_{\mathrm{h}}\right)=\mathrm{0},(9)

{\rho }_{\mathrm{w}}\frac{\partial u_{\mathrm{h}}}{\partial t}=-\nabla p-\frac{\mathrm{1}}{\mathrm{2}}{f}_{\mathrm{D}}\frac{{\rho }_{\mathrm{w}}}{d_{\mathrm{h}}}\left|u_{\mathrm{h}}\right|u_{\mathrm{h}}+F.(10)

式中：达西摩擦系数f_{\mathrm{D}}采用Churchill^［13］关系式，该公式可用于层流、湍流及层流和湍流之间的过渡区域，由管道直径和表面粗糙度决定，表示为雷诺数的函数；F为体积力.

不可压缩流体在管道中流动的能量方程为

\begin{array}{c}{\rho }_{\mathrm{w}}{A}_{\mathrm{c}}{c}_{\mathrm{w}}\frac{\partial \theta }{\partial t}+{\rho }_{\mathrm{w}}{A}_{\mathrm{c}}{c}_{\mathrm{w}}{u}_{\mathrm{h}}\nabla \theta =\nabla ·\left(A_{\mathrm{c}}{\lambda }_{\mathrm{w}}\nabla \theta \right)+\\ \frac{\mathrm{1}}{\mathrm{2}}{f}_{\mathrm{D}}\frac{{\rho }_{\mathrm{w}}{A}_{\mathrm{c}}}{d_{\mathrm{h}}}\left|u_{\mathrm{h}}\right|{u_{\mathrm{h}}}^{\mathrm{2}}+Q_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}.\end{array}(11)

其中Q_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}为管道与周围热交换产生的热源项，在岩层与井筒耦合中，井筒与周围地层存在径向传热，井筒与周围基岩单位长度的径向热传导为

\frac{\partial Q_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}}{\partial z}=\frac{\mathrm{2}\mathrm{\pi }{\lambda }_{\mathrm{s}}\left(\theta -{\theta }_{\mathrm{h}}\right)}{f\left(t\right)}.(12)

式中：z为深度；f(t)为瞬态导热系数，选用文献［14］定义的瞬态导热函数关系式：

f\left(t\right)=\mathrm{l}\mathrm{n}\left(\frac{\mathrm{4}\sqrt[]{\alpha t}}{d_{\mathrm{h}}}\right)-\mathrm{0.288}\mathrm{6},(13)

\alpha =\frac{{\lambda }_{\mathrm{s}}}{{\rho }_{\mathrm{s}}{c}_{\mathrm{s}}}.(14)

抽热比^［15］\gamma \left(t\right)定义为累计抽热量除以储存或最大可抽热量.

\gamma \left(t\right)=\frac{{\displaystyle \underset{V}{\overset{}{\iiint }}}\frac{T_{\mathrm{s}}-T_{\mathrm{s}}\left(t\right)}{T_{\mathrm{s}}-T_{\mathrm{i}\mathrm{n}\mathrm{j}}}\mathrm{d}V}{{\displaystyle \underset{V}{\overset{}{\iiint }}}\mathrm{d}V}.(15)

式中：T_{\mathrm{s}}为初始岩层温度；T_{\mathrm{i}\mathrm{n}\mathrm{j}}为初始注入温度.

井-岩热流耦合关系如图2所示，在生产过程中，井筒内流体温度不断变化，导致流体密度和黏度发生改变，影响井筒内流动状态.同时，井筒内流动会为热对流提供速度场，进而影响热对流强度.通过耦合变量p，T数据的实时传递，可以实现井筒内流动与传热的双向耦合.此外，由于井筒初始注入流体温度较低，在温差作用下，井筒周围地层会持续向井筒进行热量补给.而井筒内流体与周围地层的传热会改变地层内温度场的分布，引起地层内流体物性的变化，产生浮升力，形成自然对流，通过改变速度场影响热对流强度.基于双孔隙度的地层渗流通过热对流提供速度场影响传热强度.取热过程中，储层内流体温度的变化导致流体的密度和黏度发生改变，进而影响其渗流状态.此外，地层内流动传热不会直接影响井筒中流体，而是通过影响井内传热过程来间接影响井内流体流动.

1.3 边界条件与初始条件

温度场边界条件与初始条件：模型上表面温度为20 ℃，模型四周绝热，温度梯度为35.1 ℃/km，注水井入口指定注入温度20 ℃，出口为热流出；采出井底温度与热储区域温度数据耦合，出口为热流出.渗流场边界条件与初始条件：模型上表面压力边界条件设置为大气压，模型四周为不可渗透边界条件.采用定流量边界条件，注水井入口流量20 kg/s，注水井底压力与热储区域达西压力数据耦合，采出井口为大气压.本研究以地热开采40 a时间作为研究期限，计算时间步长设置为0.5 a，相对容差为10^-7，作为数值计算的收敛条件.采用全耦合算法求解数学方程.

1.4 参数设置

数值模拟中用到的物理参数见表2，其中岩层数据参考辽宁省地质条件设定^［16］，工质水采用COMSOL内置，其物性随温度发生变化.

2 模型及网格无关性验证

采用扫掠法对热储部分建立三棱柱体网格，对基岩生成四面体网格，其中针对注采井及裂缝区域进行网格加密，提高计算精度.为获得计算结果准确的网格，以40 a内出口温度变化作为评估指标，结果如图3所示，确定使用的网格单元数为612 898，网格质量为0.7，符合计算要求.

采用Ghassemi等^［18］提出的裂隙面温度分布解析式来验证数值解的有效性.500 m×500 m裂隙面平均温度随裂隙长度分布的数值解与解析解对比如图4所示，结果表明本文模型具有较高的准确性.

3 结果分析

3.1 一注一采井岩耦合温度场及渗流场分析

以表2给定的参数进行数值模拟得出系统运行1，10，20，40 a时井所在面及裂隙面温度云图如图5所示.由图5可知，随着年限增加，注水井内排出温度相对较低的冷流体会迅速冷却临近注水井井眼的干热岩体，在热传导作用下，冷区域不断由注水井向生产井扩散，形成低温区.对比之下，生产井周围的温度变化不明显，说明系统可以在一定时间跨度内保持稳定产热，具有较长的热采寿命.随时间推移，更多热量从储存库中被提取出来，冷区域向生产井眼靠近，达西渗流场呈现如图5e~5h所示的规律.但井身周围温度变化范围很小（图5a~5d），其热区域半径仅为10 m左右，因此该系统主要取热的通道是裂隙，其次是井筒.

如图6所示，EGS取热通道依次为注水井—深层裂隙—采出井.其中，注水井内初始温度20 ℃的工质水与高温地层实现间接换热，到达井底的水温最高达到70 ℃，可见以往研究仅将注水井作为入口条件，而不考虑其热收入的情况过于理想.换热工质经注水井排出进入深层裂隙换热，换热后水温能达到160 ℃；最后从采出井采出时，由于向上输出过程中地温梯度递减，工质热量部分散失，出口温度最高达到143 ℃.随时间进行，注水井及采出井内温差减小，这是由于井筒与周围地层间热传导削弱，冷区域扩散.因此，整个换热过程依靠裂隙及井筒的换热作用进行，其中裂隙贡献率近73%，井筒换热占比为27%.

3.2 井间距、注水流量的影响

3.2.1 注采井间距影响

采用100~500 m井间距分析其对采热的影响如图7a所示，随注采井间距增加，出口温度不断提高，采热率也随之缓慢提升.200~400 m井间距在40 a后热量衰减比例依次为25%，9%和2%，500 m井间距下生产温度始终保持缓慢上升趋势，几乎不衰减.井间距为100 m时，开采寿命仅为2 a，随后生产井内温度迅速下降，最后趋于平缓下降；随注采井间距由100 m增加到500 m时，开采寿命延长了38 a.同时由图7b可知，采热率随井间距逐渐升高.500 m间距时采热效果最好，出口温度最高达到145.5 ℃，证明井间距越长越有利于高温热开采及提高热储寿命.

3.2.2 注入流量影响

本节选定500 m作为井间距，研究10~50 kg/s注水流量变化对采热的影响.如图8a所示，当注水流量为10 kg/s时，出口温度在40 a时为131 ℃并始终保持缓慢上升趋势，而注水流量50 kg/s情况下出口温度在第10年达到最高153.63 ℃，此后开始下降.出口温度在前20 a间随注水流量增加而线性增加；20 a后，20~50 kg/s情况下出现不同程度的下降，其中50 kg/s下降最快，但其采热率在同等水平下最高，如图8b所示.为了充分利用地热资源，建议采用20~30 kg/s的注水流量为宜.

3.3 井布置方案影响

多井布置能更充分地利用热储资源.由图5可知，一注一采布置下，井筒内工质流入热储，沿裂隙多方向流动散失，冷区域以注水井眼为中心扩散，只有部分回收至生产井，随时间增加，井身周围热扩散半径增大，但最大不超过10 m.因此，研究多井布置时，井与井间可视为互不影响.研究井布置对于优化流场分布及提高采热率起着至关重要的作用.为验证井布置对采热的影响，在现有一注一采基础上，设计了一注两采、一注三采、一注四采3种井布置方案，如图9所示.采用注水流量20 kg/s，40 a时的井截面及裂隙面温度场云图如图10所示，冷区域总是优先从注水井向采出井扩散，其渗流场流线随井眼位置变化呈现出规律的对称分布.

对比40 a的生产温度及采热量来评价多井布局的性能如图11所示.其中一注一采的出口温度始终保持最高且热开采寿命最长，而多采情况下，寿命衰减速度加快，稳定开采时间不足10 a，这是由于一注一采模式可以保持最佳的井间距，有利于冷工质与热储间的充分换热.但多采情况受热储区域的限制，井间间距也受限.采热量随采出井数增加而升高，随时间呈先升高后降低的趋势，其中一注四采情况下，40 a时年采热量仍能达到5 MW以上，是一注一采热量的3.1倍.结合云图变化可知，一注两采冷区域扩散范围最大，但更多的热量散失于岩体中无法回收，而一注四采虽然扩散范围最小，却均匀地扩散至4座生产井所在位置，因此该情况下采热量最高.

4 结 论

1） EGS的主要取热通道是裂隙，其次是井筒，整个采热过程依靠裂隙及井筒的换热作用进行.其中，裂隙发挥了73%的作用，井筒换热占27%.在后续研究中，建议同时对注水井筒的换热及采出井筒的保温作用进行强化.

2） 在一注一采情况下，井间距为500 ｍ时采热效果最好，出口温度达到145.5 ℃，在40 a间几乎不衰减.可见井间距长有利于高温热开采及提高热储寿命；注水流量对EGS影响较大，注入流量越高，越有利于在短时间内充分开发HDR储层热能，但只追求快速完全开发不利于HDR资源的可持续利用.

3） 井布置方式影响冷工质与热岩间的热交换能力.一注一采布置具有最长的热开采寿命和较高的出口温度.而一注四采具有最高的采热量，是一注一采时的3.1倍.通过优化EGS生产井布局，可以减少优先流动路径的注水，使大部分热储层的流体参与循环，提取足够多的理想热量.

参考文献

[1] Chu H Q， Huang Z， Zhang Z K， et al. Integration of carbon emission reduction policies and technologies： research progress on carbon capture， utilization and storage technologies［J］. Separation and Purification Technology， 2024， 343： 127153.

[2] Shan D D， Yan T， Li W， et al. A review of multi-field coupling simulation of wellbore and heat reservoir in enhanced geothermal system［J］. IOP Conference Series： Earth and Environmental Science， 2019， 300（2）： 022076.

[3] Ma W W， Yang C， Ahmed S F， et al. Effects of thermophysical parameters of fracturing fluid on hot dry rock damage in hydraulic fracturing［J］. Geomechanics for Energy and the Environment， 2022， 32： 100405.

[4] Sun Y X， Zhang X， Liu Y S， et al. Prediction and optimization of productivity and lifespan in multi-well enhanced geothermal system［J］. Applied Thermal Engineering， 2023， 234： 121155.

[5] Wang Z　Ｗ， Zhang B， Zhu X　Ｙ， et al. How does the space of production wells influence the heat extraction efficiency？ Defined by a 3D modeling work in enhanced geothermal system （EGS） ［J］. Frontiers in Energy Research， 2024， 12：1361290.

[6] Chen J L， Jiang F M. Designing multi-well layout for enhanced geothermal system to better exploit hot dry rock geothermal energy［J］. Renewable Energy， 2015， 74： 37-48.

[7] Huang W B， Cao W J， Jiang F M. Heat extraction performance of EGS with heterogeneous reservoir： a numerical evaluation［J］. International Journal of Heat and Mass Transfer， 2017， 108： 645-657.

[8] Lu X， Tong X L， Du X P， et al. Effect of wellbore layout and varying flow rate on fluid flow and heat transfer of deep geothermal mining system［J］. Thermal Science and Engineering Progress， 2023， 42： 101870.

[9] Yao J， Zhang X， Sun Z X， et al. Numerical simulation of the heat extraction in 3D-EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model［J］. Geothermics， 2018， 74： 19-34.

[10] 单丹丹， 闫铁， 李玮， 等. 单裂隙热储热流耦合数值模拟分析［J］. 当代化工， 2020， 49（4）： 716-719， 723.

[11] Shan Dan-dan， Yan Tie， Li Wei， et al. Numerical simulation and analysis of thermal-hydraulic coupling in a single-fracture thermal reservoir［J］. Contemporary Chemical Industry， 2020， 49（4）： 716-719， 723.

[12] Aliyu M D， Archer R A. A thermo-hydro-mechanical model of a hot dry rock geothermal reservoir ［J］. Renewable Energy， 2021， 176： 475-493.

[13] Zhao K， Song W J， Wang X Y， et al. Study on the influence law of well location and water injection displacement on heat extraction performance of EGS［J］. Energy Science & Engineering， 2023， 11（7）： 2571-2584.

[14] Churchill S W. Friction-factor equation spans all fluid-flow regimes ［J］. Chemical Engineering， 1977， 84（24）： 91-92.

[15] Keller H H， Couch E J， Berry P M. Temperature distribution in circulating mud columns［J］. Society of Petroleum Engineers Journal， 1973， 13（1）： 23-30.

[16] Jiang F M， Luo L， Chen J L. A novel three-dimensional transient model for subsurface heat exchange in enhanced geothermal systems［J］. International Communications in Heat and Mass Transfer， 2013， 41： 57-62.

[17] 于小磊， 张晓辉， 刘振伟， 等. 辽宁省干热岩资源赋存条件分析及靶区初选［J］. 地质与资源， 2023， 32（5）： 608-615.

[18] Yu Xiao-lei， Zhang Xiao-hui， Liu Zhen-wei， et al. Hot dry rocks in Liaoning Province： occurrence conditions and primary selection of target areas［J］. Geology and Resources， 2023， 32（5）： 608-615.

[19] Zeng Y C， Su Z， Wu N Y. Numerical simulation of heat production potential from hot dry rock by water circulating through two horizontal wells at Desert Peak geothermal field ［J］. Energy， 2013， 56： 92-107.

[20] Ghassemi A， Tarasovs S， Cheng A H D. Integral equation solution of heat extraction-induced thermal stress in enhanced geothermal reservoirs［J］. International Journal for Numerical and Analytical Methods in Geomechanics， 2005， 29（8）： 829-844.

基金资助

辽宁省科技重点研发项目(2022JH2/101300109)