基于本征正交分解的气动优化设计外形数据挖掘∗

段焰辉 吴文华 范召林 罗佳奇

1)(中国空气动力研究与发展中心,计算空气动力学研究所,绵阳 621000)2)(北京大学工学院,北京 100871)

基于本征正交分解的气动优化设计外形数据挖掘∗

段焰辉1)†吴文华1)范召林1)罗佳奇2)‡

1)(中国空气动力研究与发展中心,计算空气动力学研究所,绵阳 621000)2)(北京大学工学院,北京 100871)

(2017年7月2日收到;2017年7月18日收到修改稿)

气动外形的全局优化设计会产生大量的过程数据,其中隐含的设计知识具有较高的挖掘价值.数据挖掘有助于获取直观、可定性描述的设计知识.本文采用基于本征正交分解的数据挖掘方法从气动优化设计的过程数据中获取设计知识,数据挖掘对象为跨音速压气机转子叶片NASA Rotor 37的优化过程数据,该数据由基于粒子群方法的绝热效率最大化优化设计产生.结果表明:基于本文数据挖掘方法获取的设计知识能够直接反映气动外形的变化规律,为叶片的气动外形设计提供参考;数据挖掘的设计知识成功地验证了优化设计结果的有效性.

数据挖掘,本征正交分解,气动优化设计,跨音速

1 引 言

气动外形优化设计是航空航天领域的重要学科之一,相较传统的设计方法,具有成本低、周期短的突出优点.随着气动外形优化设计技术的发展,随机搜索类方法因为能够克服传统优化方法局部性的缺陷,得到了日益广泛的应用.这类优化方法在优化设计中将产生大量的过程数据,但是只有最终的优化结果被保留应用.耗费大量计算资源的过程数据并非毫无用处,其中隐含了与优化相关的设计知识,即设计变量与目标函数的关系.若能获取这些数据中的设计知识并加以利用,有助于揭示设计空间内设计变量变化对目标函数影响的一般规律,更深层次地分析认识设计变量对目标函数的影响机制,为类似的设计优化积累设计信息和经验,对提高类似设计优化的效率与可靠性具有重要的指导意义.本质上,数据挖掘就是这样一种从大量数据中获取隐含知识的方法.本文将利用数据挖掘技术获取优化设计过程数据中的设计知识.

目前得到应用的数据挖掘方法较多,应用领域极广.在优化设计方面,数据挖掘技术多用于多学科、多目标优化设计的设计知识挖掘[1].Jeong等[2]、Chiba和Obayashi[3]采用总变差分析和自组织映射方法对两级入轨飞行器可飞回第一级的多学科、多目标优化设计结果进行了分析,实现了从非劣解集中选择最终结果;Oyama等[4,5]采用基于本征正交分解(proper orthogonal decomposition,POD)[6]的数据挖掘技术对跨音速翼型的多目标优化设计结果进行了分析,将优化结果分成了三类,并获取了每类结果的设计知识.郭振东等[7]对数据挖掘的应用进行了拓展,采用显著变量识别、总变差分析及自组织映射等数据挖掘技术对设计空间进行知识挖掘,得到的设计知识与优化设计结果进行了相互验证.汪伟等[8]将优化设计的过程数据作为数据挖掘对象,采用K-Means聚类分析、粗糙集属性重要性算法和决策树方法对优化过程数据进行处理和分析,得到了叶片优化设计变量与目标函数的隐含规则.综上可知,数据挖掘在气动外形优化设计中的应用尚处于探索阶段,分析方法和对象呈现多样化,具有重要的研究意义.

本文基于POD的数据挖掘技术应用于优化设计过程数据的处理,以获取设计知识.为了提高数据挖掘对象的可信度,也即数值精度,本文的优化过程数据都由计算流体力学(computational fluid dynamics,CFD)方法确定.为了缩短优化周期,采用同步并行粒子群方法(parallel synchronous particle swarm optimization,PSPSO)在大规模集群上进行优化.此外,当数据挖掘对象为气动优化设计的设计变量,结果的直观性普遍较差,因为多数参数化方法提取的设计变量都无法直接反映气动外形的特征.基于POD的数据挖掘可以直接对外形数据进行分析,数据挖掘的结果能直接反映外形的变化,提供直观的设计知识,因而本文采用POD进行数据挖掘,改善结果的直观性.

本文首先介绍POD的基本原理及基于POD的数据挖掘方法.然后介绍本文采用的PSPSO气动外形优化设计方法,以及三维叶片的参数化方法和网格变形方法;并采用PSPSO方法对跨音速压气机转子叶片NASA Rotor 37进行气动优化设计,提取过程数据.最后采用POD对过程数据进行挖掘,分析所隐含的设计知识,并验证本文所提出的数据挖掘技术的有效性.

2 POD数据挖掘方法

数据挖掘方法一般可以分为三步:数据制造、数据预处理和数据挖掘.本文对优化设计的过程数据进行数据挖掘,数据制造对应优化设计,数据预处理对应制造数据的筛选和处理,数据挖掘对应优化气动外形的特征分析.本文的数据预处理主要是按目标函数的优劣对优化过程数据进行筛选,具体做法是:与原始外形的气动性能相比,将气动性能提高的过程数据筛选出来,剔除气动性能降低的过程数据.值得注意的是:优化设计中得到的过程数据包含气动外形、流场等信息,在数据挖掘时并非所有数据都要使用,因而在预处理阶段需要对这些数据予以区分.本文对气动外形进行数据挖掘,因此在预处理阶段只保留气动外形信息.

2.1 POD基本原理

原始POD方法[6]需要整个定义域中的解集,维度较高,对特征值求解的精度和稳定性提出了较大挑战,难以在工程中得到应用.目前多采用由Sirvoich提出的Snapshot POD方法[9],该方法用一组定义域中的快照代替整个定义域中的解集,提高了POD方法的稳定性,降低了计算量,目前在气动研究中得到了大量应用[10−15].下文中,若无特殊说明,所涉及的“POD方法”都为“Snapshot POD方法”.

假设在某一定义域内有一组快照:s1,s2,···,sm,其POD基为:φ1,φ2,···,φm,m为快照数.POD基的求解方法如下.

首先求解特征值问题,

其中R为快照集合的自相关矩阵,λ为特征值,a为特征向量.R定义为

由奇异值分解(singular value decomposition,SVD)确定R的特征值和特征向量后,POD基可由下式计算:

其中:φk表示第k个基,表示第k个特征值对应的特征向量的第i个元素.

2.2 基于POD的数据挖掘

本文基于POD的数据挖掘与应用主要采用文献[12,13]介绍的基于POD的参数化方法,在此基础上加入一些能够反映物理特征的量来挖掘设计知识.本文POD分析对象为优化过程中的气动外形.

假设被分析的气动外形由表面网格点组成,定义一组向量g1,g2,···,gm,对于优化设计的过程数据而言,这组向量应该具有某种共性,比如优化的目标函数都高于初始值.每个向量中的数据为气动外形表面网格点三维坐标的顺序排列,如gi={X,Y,Z}T,X,Y,Z依次为网格点的三个坐标值,并且按照网格点编号排成序列.若气动外形表面网格点数为n,则每个向量的维数为3n.由气动外形表面网格点坐标组成的向量可以确定POD所需的快照:

其中,g0可以选为网格点坐标向量的平均值,也可以是优化的初始外形.每个外形快照减去g0,便于分析所有气动外形的扰动特征.由SVD确定POD基后,所有快照都可由这些基表示:

数据挖掘的设计知识就隐含在POD基和系数(也称为基的权重)中.基能够反映气动外形变化的区域及相对大小,结合其系数就能得到气动外形变化的绝对大小,从而给出最终的气动外形变化规律.但是基较多(与快照数目相同),系数也较多,一般只需要对其中能够反映主要规律的基和相应系数进行分析即可.如何筛选基和系数,是基于POD数据挖掘方法的一个重点.

2.2.1 POD基及基系数的筛选

引入基的能量的概念[4]对基进行筛选,第i个POD基的能量Energy定义为

其中φi,j表示第i个POD基的第j个元素.由(6)式可知:能量其实就是每个外形的网格点位移大小的总体度量.再定义能量比重:

由能量比重可以确定每个POD基所占的比重,从而确定哪些基是获取设计知识的关键.

基系数的筛选要解决两个问题:分析哪些基的系数,分析哪些快照的系数.本文研究中,只对筛选得到的POD基的系数进行分析.快照则根据优化目标函数的大小,选择靠前的部分快照,这样可以保证得到的基的系数反映的都是最优外形的设计规律.

2.2.2 设计知识的获取方法

POD基及由POD基和基系数组合表示的气动外形,都是坐标的扰动信息,不易直接获得设计知识,需要将其变换成能够反映物理特征的量.气动外形的变化规律实际上就是每个网格点位置的变化规律,这种变化不仅有大小信息,还有方向信息.确定这种变化量的最佳办法是采用优化时的参数化方法,根据参数化方法中每个网格点位置扰动的施加方式(一般是在某个方向施加),反推网格点的位置变化.本文采用了前期研究中的参数化方法[16],扰动量沿着气动外形表面当地网格点的外法线方向施加,因此第i个POD基中第j个网格点的位置变化量(下文统称为r-变量)为

其中:n为该网格点外法线方向.下文的设计知识主要是通过对r-变量进行分析获取.

3 气动外形优化设计

3.1 同步并行粒子群方法

PSPSO方法是缩短基本PSO方法[9]优化周期的有效手段.图1为算法流程图,其中n为粒子数目,m为优化步数,i为粒子编号,k为当前迭代步,f和c分别为目标气动函数和约束函数,pi为第i个粒子进化到第k代时的最优粒子,为第k代中所有粒子中的最优粒子.粒子位置x和速度v的初始化方法与基本PSO相同,计算惯性权重因子采用基于差异系数(coefficient of variation,COV)的衰减方法[17].

图1 PSPSO优化流程图Fig.1.Flowchart of PSPSO-based optimization.

3.2 参数化和网格变形方法

本文参数化方法是将叶片分解为多个特征截面,采用Hicks-Henne型函数方法[18]对特征截面进行参数化,最后再通过插值得到整个三维叶片外形.本文参数化是在初始外形上叠加几何扰动量以产生新的外形,因此无需拟合初始外形,避免了在拟合初始外形时引入误差.网格变形采用代数插值方法,将叶片表面的变化映射到空间网格,得到新外形的网格.为了提高计算效率,网格变形计算只对叶片相邻的网格块进行.具体的参数化和网格变形方法见文献[16].

4 Rotor 37数据挖掘算例

4.1 CFD数值验证及优化问题描述

流场计算采用团队开发的大规模并行计算程序ASPAC,采用多块结构化网格.该程序求解雷诺平均Navier-Stokes方程,采用固连在转子上的旋转直角坐标系来描述流动.湍流模型选用一方程Spalart-Allmaras模型[19].无黏通量采用Roe格式离散,黏性项采用二阶中心差分格式离散.本文选择跨音速压气机转子叶片NASA Rotor 37进行研究,该模型试验数据比较详细,是叶轮机CFD领域常用的验证算例之一,试验测量站位、空间分布及结果见文献[20],设计参数见表1.

表1 Rotor 37设计参数Table 1.Design parameters of Rotor 37.

在出现旋转失速之前,转子叶片流场满足周期性边界条件,因此计算域仅选取一个叶片通道.图2为O4H拓扑的计算网格.叶高方向分布56个网格单元,叶顶间隙内径向网格单元为16.流向和周向网格单元分别为152和60,总网格单元约为90万.

图2 Rotor 37计算网格 (a)S1流面;(b)S2流面Fig.2.Grid of Rotor 37:(a)In S1 streamsurface;(b)in S2 streamsurface.

图3 (网刊彩色)Rotor 37特征工作线 (a)总压比;(b)绝热效率Fig.3. (color online)Performance characteristics:(a)Total pressure ratio;(b)adiabatic efficiency.

研究中还采用商业软件NUMECA来进行流场数值模拟.ASPAC堵塞流量为20.95 kg/s,NUMECA堵塞流量为20.93 kg/s,均在试验结果的误差范围内.图3对比了CFD和试验测量的总压、效率特征工作线.由图3可知:基于CFD计算的总压比均与试验较为接近,而CFD效率均与试验存在较明显的偏差.截至目前,很少有数值模拟的效率与试验结果符合较好.Denton在[21]对Rotor 37的“盲测”结果进行分析时指出:试验中端壁附近流动损失无法测量可能是造成这种偏差的主要原因.通过与NUMECA计算结果的对比分析可知:ASPAC特征工作线稍优于NUMECA,因而基于ASPAC的计算结果是合理的.另一方面,对于气动外形优化设计,更为重要的是能够正确反映出气动性能随设计变量的变化趋势.

4.2 基于粒子群的优化设计

本文在Rotor 37的设计转速近效率峰处进行外形优化设计,同时约束近效率峰的总压比和流量,优化问题可以由下式描述:

其中:η为等熵效率;m和m0分别为设计流量和目标流量,目标流量为原始外形的流量;π和π0分别为设计总压比和目标总压比.

选择位于100%,75%,50%,25%,0%叶高的五个特征截面进行参数化,选择中间三个截面进行设计.一般情况下,吸力面对绝热效率的影响更大,只在吸力面布置10个型函数,共30个设计变量.设计变量的变化范围定义为(−0.1,0.1).粒子群优化方法采用80个粒子,每个粒子使用8个CPU计算,一共使用640个CPU,共优化25代.

表2给出了优化前后的气动性能参数,绝热效率由0.8588提高到0.8692,提高了1.21%,总压比降低了0.46%,流量提高了0.43%,均在约束条件之内.

表2 优化前后叶片的气动性能Table 2.Aerodynamic performance of reference and optimized blades.

图4给出了叶片吸力面上的压力分布对比,50%叶高以上激波明显向叶片尾缘移动,从而影响激波诱导的分离.图5给出了吸力面附近的流线分布,图中蓝色流线为优化后的流线,由图可知:优化后的激波位置明显后移,相应的分离区域明显减小.

图5 (网刊彩色)叶片吸力面流线(相对速度):原始(红色);优化(蓝色)Fig.5.(color online)Streamlines(relative velocity)on the suction surface of reference(red)and optimized(blue)blades.

图6给出了25%,50%,75%叶高的气动外形和压力分布.25%叶高处,从弦长中部开始叶片厚度开始减小,在尾缘附近先增大后减小,整体上最大厚度后移,能延迟激波的产生.由压力分布也可以看出:优化后激波明显后移.50%叶高处,在尾缘附近叶片厚度明显增加,最大厚度后移,激波也随之后移,且激波后移更多.75%叶高处,从弦长中部开始叶片厚度开始增大,在尾缘附近先减小后增大,但变化非常小.相对于25%叶高,75%叶高的优化叶型虽然没有明显的最大厚度后移,但仍能延迟激波的产生.

由上述分析可知,Rotor 37绝热效率的提高主要源自激波位置的后移,以及吸力面流动分离区域的减小,导致这种变化的主要原因是叶型吸力面外形变化.

图6 (网刊彩色)不同叶高的气动外形和压力分布 (a)25%叶高;(b)50%叶高;(c)75%叶高Fig.6.(color online)Aerodynamic shape and pressure distributions at(a)25%span,(b)50%span,(c)75%span.

4.3 设计知识挖掘

优化共进行了25步,加上初始化粒子,共有2080个粒子,其中气动外形奇异和气动特性不满足优化要求的粒子在进行筛选时需要剔除.本文旨在获取与优化设计目标相关的设计知识,筛选数据的规则为:满足优化约束且绝热效率提高.筛选出有效粒子1734个,并按照绝热效率由高至低进行排列.图7给出了所有粒子的分类情况,其中:＜0%表示不满足规则.

对满足规则的数据按照(4)式进行预处理获得快照,并由SVD确定1734个POD基.POD基及特征值的特点是:特征值大的POD基包含最明显的特征,本文采用基的能量比重对有效基进行筛选,由(6)和(7)式计算的前10个POD基的能量比重如图8所示.由图可知,前三个基分别占了94.62%,3.65%和1.09%的总能量,99.36%的能量都集中在前三个基上面,因而下文主要对前三个基进行分析.

图7 (网刊彩色)所有粒子的分类情况Fig.7.(color online)Classi fication of proceeding data.

图8 前10个基的能量比重Fig.8.Energy ratio of the first ten POD modes.

POD基的分析对象为基的r-变量,r-变量的正负并不直接反映叶片外形厚度的增大或减小,还要考虑基系数的影响.本文将r-变量的正负描述为正向变化和负向变化.由上文所述,从基能够获得的设计知识为气动外形变化量的相对大小和所处区域.为了更好地描述外形变化量的位置,本文将叶片吸力面划分为9个区域,划分方法为:沿叶高选择25%,50%,75%三个截面,以1,2,3命名;沿弦向将这三个截面划分为三个大致相等的区域,以A,B,C命名.图9给出了9个区域的r-变量分布.由图9可知:外形变形量大的区域位于B1,C1和B2,C2区域;变形复杂的区域位于B3,C3区域;A区域整体变形都很小.需要注意的是:r-变量的比较分析都是基于相同的原始外形几何尺寸.根据上述设计知识可知:对于此类优化设计问题,设计变量应重点分布在B,C区域,设计变量的设计空间在截面1,2处应较大,在截面3处应较小.

接下来分析基系数的影响.由于前三个POD基所占能量已经超过99%,本文用前三个基来描述目标函数排序靠前的100个快照.图10给出了快照的前三个基系数.由图可知:三个基系数的绝对值依次减小,与基所占能量的大小变化趋势一致;此外,基系数都在某一固定值附近波动,该固定值应该能够反映基系数中隐含的主要设计规律.

将基系数分解为固定值和扰动的迭加形式:

其中αm表示前100个快照的第m个基的系数组成的向量,为该向量的固定值,为该向量扰动值.采用线性回归的方法确定上述公式中的固定值,可得:

其中n表示第n个快照.

图10中不同颜色的虚线对应各个基系数的线性回归结果.虽然每个快照的固定值略有不同,但是差别非常小,可以忽略,因而

此外,三组基系数的扰动量最大值分别为:

相对于固定值,前两个基系数扰动量的最大值较小,可忽略其对基权重的影响.第三个基系数扰动量与固定值量级相同,但是考虑第三个基的权重仅占总能量的1.09%,其对外形的影响很小,基本可以忽略第三个基系数的扰动量的影响.因而,可以仅由POD基及基系数的固定值确定“拟合优化外形”,如下式所示:

其中φcomb表示POD基组合,也即拟合优化外形.

图9 (网刊彩色)r-变量分布云图 (a)第一个基;(b)第二个基;(c)第三个基Fig.9.(color online)Contours of r-variation on:(a)The 1st mode;(b)the 2nd mode;(c)the 3rd mode.

图10 (网刊彩色)前三个基的系数分布Fig.10.(color online)Coefficient distributions of the if rst three basis modes.

图11给出了拟合优化外形及优化外形的r-变量分布.由图可知:截面1,2,3上两者的变化趋势基本相同;截面1上,C1区域中部气动外形变化的最大值略有差别;截面2上,A2区域优化外形减薄区域较拟合优化外形的小,气动外形变化的最大值虽有差别,但都接近零,C2区域气动外形增厚的最大值略有差别;截面3上,B3区域气动外形增厚的最大值略有差别,C3区域气动外形增厚、减薄的最大值也略有差别.

由前三个POD基和基系数的固定值线性拟合得到的r-变量的变化规律与真实优化外形的r-变量变化规律基本一致,对优化结果和数据挖掘结果起到了相互验证的作用.

5 结 论

1)数据挖掘得到的设计知识与优化设计结果基本一致,证实了本文采用的数据挖掘方法和优化设计结果均是可靠的,验证了本文所采用的优化设计方法和CFD程序;

2)基于POD的数据挖掘方法所获取的设计知识能够直观反映气动外形的变化规律,具体到本文的叶片算例,明确地反映了设计变量的布置位置和变化幅度,同类型叶片均可按此规律进行气动外形的优化设计,具有较大的指导意义;

3)数据挖掘的目标还可以进一步扩展,不局限于优化设计目标,比如:可以从效率优化过程数据中提取压比增大的过程数据,对压比增大的气动外形进行数据挖掘.整体上,提高数据挖掘目标的多样性,能够进一步发掘优化过程数据的隐含价值.

[1]Jeong S,Shimoyama K 2011Proc.Inst.Mech.Eng.Part G:J.Aerosp.Eng.225 469

[2]Jeong S,Chiba K,Obayashi S 2005J.Aeros.Comp.Inf.Com.2 452

[3]Chiba K,Obayashi S 2008J.Spacecraft Rockets45 975

[4]Oyama A,Nonomura T,Fujii K 2010J.Aircraft47 1756

[5]Oyama A,Verburg P,Nonomura T,Harry W M,Fujii K 2010AIAA Paper2010 1140

[6]Holmes P,Lumley J L,Berkooz G 1997Q.J.Roy.Meteor.Soc.123 2500

[7]Guo Z D,Song L M,Li J,Li G J,Feng Z P 2015J.Propul.Technol.36 207(in Chinese)[郭振东,宋立明,李军,李国君,丰镇平2015推进技术36 207]

[8]Wang W,Mo R,Zhang Y 2013Comput.Eng.Appl.49 11(in Chinese)[汪伟,莫蓉,张岩2013计算机工程与应用49 11]

[9]Sirvoich L,Kirby M 1987Quart.Appl.Math.45 561

[10]Duan Y H,Cai J S,Li Y Z 2012AIAA J.50 968

[11]LeGresley P,Alonso J 2000AIAA Paper2000 2545

[12]Toal D J J,Bresslo ffN W,Keane A J,Holden C M E 2010AIAA J.48 916

[13]Ghoman S,Wang Z,Chen P,Kapania K 2012AIAA Paper2012 1808

[14]Luo J,Duan Y,Tang X,Liu F 2015ASME Paper2015 42876

[15]Luo J Q,Duan Y H,Xia Z H 2016Acta Phys.Sin.65 124702(in Chinese)[罗佳奇,段焰辉,夏振华2016物理学报65 124702]

[16]Duan Y H,Wu W H,Fan Z L,Chen T 2016ASME Paper2016 56861

[17]Venter G,Sobieszczanski-Sobieski J 2003AIAA J.41 1583

[18]Hicks R M,Henne P A 1987J.Aircraft15 407

[19]Spalart P R A,Allmaras S 1992AIAA Paper1992 0439

[20]Reid L,Moore R D 1978NASA TP1978 1337

[21]Denton J D 1998J.Therm.Sci.6 1

PACS:02.60.Pn,47.85.Gj,47.11.—j,29.85.CaDOI:10.7498/aps.66.220203

*Projects supported by The National Nature Science Foundation of China(Grant Nos.51676003,51206003).

†Corresponding author.E-mail:duanyanhui@foxmail.com

‡Corresponding author.E-mail:jiaqil@pku.edu.cn

Proper orthogonal decomposition-based data mining of aerodynamic shape for design optimization∗

Duan Yan-Hui1)†Wu Wen-Hua1)Fan Zhao-Lin1)Luo Jia-Qi2)‡

1)(China Aerodynamic Research and Development Center,Computational Aerodynamics Research Institute,

Mianyang 621000,China)2)(College of Engineering,Peking University,Beijing 100871,China)

2 July 2017;revised manuscript

18 July 2017)

Global optimization methods are becoming more and more important in aerodynamic shape optimization.A large number of proceeding data will be generated during design optimization,from which the implicit but valuable design knowledge can be extracted.The design knowledge can then be used to help the designers to acquire the e ff ects of geometric variations on the aerodynamic performance changes.In this paper,we strive to extract the implicit design knowledge from proceeding data by a data mining method based on proper orthogonal decomposition(POD),by which the design knowledge more enriched and more visualized than those obtained from other data mining methods can be obtained.Proceeding data for data mining are ingathered from aerodynamic shape optimization of a transonic compressor rotor blade,NASA Rotor 37.The design optimization attempts to maximize the adiabatic efficiency of Rotor 37 under the operation condition near peak efficiency with the constrains of mass flow rate and total pressure ratio.The parallel synchronous particle swarm optimization method is employed to search for the optimization in the design space.The particles with improved adiabatic efficiency,while within the optimization constrain tolerances are picked up from the design optimization,which are then used for data mining.The geometric coordinates of the aerodynamic shape with respect to the ingathered particles are regarded as the snapshots.Then the POD modes of the aerodynamic shape can be obtained by singular value decomposition on the snapshots.The results show that the universal rules of geometry variations for the optimization maximizing the adiabatic efficiency of Rotor 37 can be directly visualized by the design knowledge extracted from the proceeding data by POD-based data mining technique.Furthermore,the optimization results are also veri fied by the design knowledge extracted by data mining.

data mining,proper orthogonal decomposition,aerodynamic optimization design,transonic

10.7498/aps.66.220203

∗国家自然科学基金(批准号:51676003,51206003)资助的课题.

†通信作者.E-mail:duanyanhui@foxmail.com

‡通信作者.E-mail:jiaqil@pku.edu.cn