基于桨叶调节的轴流转桨式水轮机模型飞逸数值模拟

周大庆，郭 优，姜德政

（1.河海大学能源与电气学院，江苏南京 211100；2.中国长江电力股份有限公司，湖北宜昌 443002）

基于桨叶调节的轴流转桨式水轮机模型飞逸数值模拟

周大庆1，郭 优1，姜德政2

（1.河海大学能源与电气学院，江苏南京 211100；2.中国长江电力股份有限公司，湖北宜昌 443002）

为了研究不同桨叶启闭规律对轴流转桨式水轮机飞逸过程的影响，采用三维非定常数值方法模拟了5种桨叶控制方式下的轴流转桨式水轮机模型飞逸过程，对比分析了转速、流量、力矩和压力脉动等参数随时间变化特性及桨叶表面压力分布和尾水管内流场演变规律。结果表明：以桨叶静止工况下的最大逸速为基准，在±10°内启闭桨叶对最大逸速影响范围为-6.6%～5.0%；在飞逸过程中打开桨叶会加剧外特性参数波动，尾水管中心部最大负压值可达初值的2.86倍，产生的偏心螺旋涡带诱发强烈低频脉动，不利于机组稳定；关闭桨叶可降低水流流速，减小压力脉动及改善尾水管流态，但需探究合理关闭方式以避免过大的转速最大上升值。

轴流转桨式水轮机；飞逸过程；桨叶调节；启闭规律；三维非定常数值模拟

水轮发电机组飞逸暂态过程是一种大波动不稳定过程，指水轮发电机组突然甩掉负荷，恰逢调速器失灵或导水机构拒动，机组转速迅速升高，直至转轮产生功率与机械损失功率相平衡时，转速达到该水头和导叶开度下的最大逸速。在飞逸状态下，机组动态附加载荷大增，易发生强烈的水力脉动和机组振动，会对水轮发电机组造成极大威胁。对于轴流转桨式水轮机组，在飞逸过程中往往导水机构和叶片操作机构同时失灵，两者间的协联关系遭到破坏，桨叶角度和导叶开度可以发生任意组合［1］。

在求解过渡过程问题时，传统一维方法主要采用两种求解方法：①外特性数值解法［2-4］，以完整综合特性曲线作为求解一维无黏流动的边界条件，然后通过特征线法求解；②内特性数值解法［5-7］，根据水力机械的基本动态参数推导出过渡过程的解析方程组，与有压非恒定管流基本方程组联立求解。上述两种方法在解决包括引水系统在内的整个水力系统过渡过程问题中发挥了重要作用［8-10］，但也存在依赖特性曲线和难以精细反映水轮机内流场细节等不足。

机组飞逸过程中流动的三维特性明显，用传统一维方法描述其过渡过程不够准确，有部分学者［11-14］采用三维非定常CFD方法对不同机型飞逸过程进行计算，所得机组宏观参数趋势与实验结果一致，验证了三维非定常CFD模拟的有效性，但由于机型或复杂动边界模拟方法限制，关于桨叶角度动态变化对飞逸过程影响的研究较少。

针对以上问题，本文采用一种非定常三维湍流数值模拟方法，研究在飞逸过程中轴流转桨式水轮机桨叶角度动态变化对过渡过程的影响，通过分析机组外特性和内部流场随时间变化规律，为寻求更有效的防飞逸措施和改善过渡过程品质提供参考。

1 数值计算方法

1.1 控制方程

针对飞逸过程中水轮机运行特点，流动区域采用的控制方程包括几何守恒方程、连续性方程和动量守恒方程。

几何守恒方程：

式中：n为控制体边界Ω处微元面积d S的外法向单位向量；D为边界移动速度；V为控制体体积。

连续性方程：

式中：uj为水流绝对运动速度；j为转轮区网格运动速度，水轮机其他过流部件由于网格不动，则其控制方程中j=0。

动量方程：

式中：fi为单位体积力；ρ为密度；p为水体压力；ν为运动黏性系数。

转轮旋转角速度由转动平衡方程确定：

式中：M为水流作用在转轮上的力矩；Mr为由机械和容积损失导致的阻力矩；J为机组转动惯量；ω为转轮旋转角速度。飞逸过程计算中，由式（4）计算出桨叶在每个时间步上的转速增量dω，与现时转速相加得出下一时间步上的转速，依此类推即可得到每个时间步上的转速。

1.2 湍流模型

Spalart-Allmaras模型［15-16］相对于两方程模型鲁棒性好、计算量小、占用计算机资源较少、计算准确度较高，而且该模型方程在各节点处的求解不依赖其他节点的解，因此在复杂流场计算中不要求计算网格很精细，对工程应用领域特别是水轮机计算具有较高的适应性，本文也采用此湍流模型封闭控制方程组进行模拟计算。

1.3 网格技术

为适应模型的复杂几何外形，采用拓扑结构较自由的非结构化四面体网格划分计算区域。在非定常计算中，转轮体中的网格是旋转的，且发生变形和重构，本文采用变速滑移网格技术完成转轮区与毗邻区动静交界面的信息传递，应用非结构化动网格技术进行桨叶启闭过程的模拟，通过编写用户自定义程序（udf）实现飞逸过程中桨叶公转和自转结合的复合运动。

桨叶运动角速度公式：

式中：ωx、ωy、ωz分别为桨叶在x、y、z轴向分角速度；ωb为桨叶启闭角速度；θ为桨叶绕轴旋转角度。采用速度分解方法实现该运动。

图1为桨叶区网格示意图，在桨叶和轮毂之间的间隙处，需要将桨叶不接触轮毂部分的运动轨迹扫掠面设置为内部面（interior）型边界条件，承担网格变形重组和信息传递的任务。

图1 桨叶区网格示意图

1.4 控制方程离散和定解条件

采用有限体积法离散控制方程，压力项采用二阶中心差分格式，对流项采用二阶迎风差分格式，时间离散采用一阶向后差分隐式格式，用SIMPLEC算法进行速度压力耦合求解。

定解条件：进出口均为预设总压、湍动能、湍动耗损率的压力边界；固壁处采用无滑移边界，近壁区采用标准壁面函数；以选取工况的定常计算结果作为非稳态飞逸过程计算的初始条件。

1.5 算法流程

以FLUENT软件为基础，根据飞逸过程特点对程序进行二次开发和调试，实现飞逸过程中转速与桨叶角度的动态变化，计算步骤如下：

a.先对模型机组进行稳态工况定常计算，然后进行一段时间的稳态工况非定常计算。

b.将稳态工况流场及相关动态参数等赋为非定常计算的初始值，开始非稳态飞逸过程的计算。

c.通过用户自定义程序（udf）的自带宏，读取该时间步上桨叶所受合力矩，由转轮旋转力矩平衡方程，计算转轮区下一时间步的旋转速度。

d.根据计算结果调用相关函数，将桨叶转动角速度分解为x、y、z轴向分角速度，桨叶依据分角速度进行动作。

e.桨叶调节完成后，检查动网格重构后网格质量是否达到计算要求，如满足要求则在本时间步内迭代计算至收敛。

f.完成本时间步计算后，输出并保存预先监控的瞬态参数和一定时间间隔的流场信息。

g.判断时间t是否大于预设时间tmax，若否，重复步骤c至f；若是，结束计算。

2 算例分析

图2为本文计算模型的轴流转桨式水轮机模型结构示意图，转轮直径为0.46m，其原型水轮机转轮直径为10.2m，其他几何结构相似。相关参数如下：固定导叶数为16，活动导叶数为32，桨叶数为5，模型机组转动惯量J=16.8 kg·m2。由模型试验给定的参数如下：水头为0.492m，转速为270 r/min，初始桨叶角度为5°，导叶开度为30.8mm。

图2 轴流转桨式水轮机模型结构示意图

采用5种方案控制飞逸过程中的桨叶角度，各方案数值模拟过程时长均为21.1 s，在稳态工况非定常计算1.1 s期间保持桨叶角度为5°不变，平稳过渡到飞逸过程后控制桨叶匀速启闭或保持静止，在11.1 s时停止桨叶动作，并在该条件下持续计算至21.1 s。各方案桨叶角度变化规律如表1所示。

表1 各方案桨叶角度变化规律____________

为保证网格数量的合理性，设计3种网格划分形式，网格数分别为140万、220万和300万，在相同条件下分别对方案3进行预计算，得到的机组转速变化情况如图3所示。可见网格数量对于结果局部细节有影响，但不改变整体趋势。综合考虑计算成本和准确度后，选择计算网格数为220万，主要过流部件蜗壳、导水机构、转轮、尾水管的网格单元数分别为392363、719424、652974、438036。

图3 不同网格数计算结果对比

2.1 外特性变化规律

为监测压力变化情况，在机组流道内设置若干压力测点，如图4所示：测点1位于蜗壳进口附近半高处，测点2位于固定导叶进口半高处，测点3位于活动导叶出口半高处，测点4位于顶盖顶部，测点5～9设置在尾水管进口断面，测点5位于尾水管进口圆面中心，其余测点位于尾水管边壁圆周。

图4 水轮机流道压力测点布置

经过试算比较，将时间步长定为0.005 s，为减少稳态工况计算的不必要数据分析，将进入飞逸的时刻1.1 s定为0时刻，只考察飞逸的20 s过程，以该稳态工况的转速、流量、力矩作为基准值，将各个过程测得的参数瞬态值与基准值的比值绘成图5，以反映各个参数随时间变化的规律，其中转速基准值为270 r/min，流量基准值为0.252m3/s，力矩基准值为25.1N·m。

方案3中，机组模型经过9.4 s达到最大逸速399.03 r/min，是初始转速的1.478倍，与模型动态试验［17］实测最大逸速398 r/min基本相符，误差主要由计算中未考虑机械摩擦力矩和风阻力矩所致。

在各个数值模拟过程中，机组在0时刻突然甩掉负荷，假设某种原因使得活动导叶无法正常关闭，发生飞逸后，机组转速先是在水力矩作用下迅速上升，水力矩则呈快速下降趋势，随着水力矩的逐渐减小，转速上升率下降，直至达到最大逸速。比较不同桨叶控制方案下的转速变化规律（图5（a））：方案1～5最大逸速分别为基准值的1.547、1.541、1.478、1.417和1.376倍，分别在8.905 s、10.030 s、9.400 s、6.005 s和4.655 s达到，可见桨叶角度越小到达最大逸速时间越长，最大逸速值也越大，与桨叶静止工况相比，在±10°内的桨叶角度变化对最大逸速影响范围为-6.6%～5.0%。方案1表明，桨叶角度逐渐关小过程中，桨叶升力大于阻力，转速总体在增大，但是10 s以后由于流量过小，阻力大于升力，转速则稍有下降。

图5（b）表明，流量变化受桨叶启闭规律影响明显，以桨叶静止的方案3为参照，方案1和2中流量上升速度随桨叶关闭迅速减小，关闭速度越快流量峰值越小，到达峰值后流量下降趋势趋于线性，降速由关闭速度主导，在10 s后因桨叶停止动作，流量变化迅速平稳；在方案4和5中，桨叶打开速度越快则流量上升也越快，并且在桨叶停止动作后流量呈现不稳定的波动。

综合比较各个方案可知：在飞逸过程中打开桨叶会降低转速最大上升值，但加剧流量和力矩的波动，导致流场不稳定；关闭桨叶可以优化飞逸过程中的流态和减小外特性参数波动性，但会引起最大逸速上升，因此需探究合理关闭规律以改善飞逸过程品质。

图5 飞逸过程中各参数相对值随时间变化曲线

图6 流道压力测点静压随时间变化曲线

图6为各测点静压值随时间变化情况，取对比明显的方案1、3、5分析，在飞逸过程中因导叶开度固定，各过程中转轮前测点1～4的压强均较为稳定，变化规律与流量变化规律相符，压力随流速变化稍有起伏。对位于尾水管进口断面的测点5～9，桨叶角度变化对静压值影响凸显：机组刚进入飞逸工况时，各测点静压受到因水流冲角变化引起的尾水管内反向水击影响，曲线产生“突点”，测点5因处于流动性较弱的尾水管低压区，水击导致的压力突变最为明显，其阶跃极值为0.9 kPa，同理在10 s时，因桨叶动作突变，方案1和5的压力变化曲线也发生了类似现象；断面中心测点5的压力变化主要受尾水涡带影响，各方案中随着涡带逐渐形成，测点压力在阶跃后快速降低，在1.2～1.8 s间达到最大负压约-2 kPa，随后因轴向回流作用缓慢上升，并在2～3 s间达到动态平衡，在方案3的14 s和方案5的9 s时，因涡带强度增大导致的涡带断裂使得压力曲线呈快速波动式下降，最终转为波幅较稳定的规律性波动；在尾水涡带强度较低时，边壁测点6～9的压力变化主要由快速旋转的转轮出流主导，压力值主要呈现高频脉动，随着涡带半径增大其对壁面的撞击增强，压力脉动转为有一定相位差的低频周期性波动，在涡带强度最大的方案5中，峰值可达0.7 kPa。

图7 飞逸过程中桨叶压力面（左）和吸力面（右）压力分布变化过程（单位：Pa）

大桨叶角度下尾水管水力稳定性急剧恶化，尾水涡带诱发幅值随桨叶角度增大而增大的低频脉动，甚至引起上游测点1～4压力波动，说明其传播距离远，对机组影响大，是诱发机组振动的主因；转轮的变速旋转引起的高频脉动在流态较复杂的尾水管边壁附近作用明显，这种压力脉动频率高，传播距离较短，影响区域较小。

2.2 桨叶压力分布变化过程

因各桨叶压力分布相似，故只取单片桨叶压力图进行分析，仍选具有代表性的方案1、3、5对比，分别列出2 s时、达到最大逸速时、15 s时的桨叶压力面和吸力面压力分布情况，如图7所示。

各方案中，在机组达到最大逸速前，随着转速升高，水流进口角逐渐减小，从正冲角变成负冲角，桨叶吸力面进水侧近叶缘处发生水流撞击，压力面进水侧叶缘处出现脱流，因此产生对应的高压和低压区，其影响范围与转速正相关，并随桨叶角度减小而增大。压力极值产生于机组达到最大逸速前后，大小主要由流量、转速和瞬时流态决定，且小桨叶角度对应较大压力极值，极大值约22 kPa，极小值为-30 kPa，与所记录飞逸初始时刻压力极大值9.2 kPa、极小值-10.2 kPa比较，分别可达初值的2.39和2.94倍。

在飞逸过程趋于稳定后，桨叶压力分布取决于对应桨叶角度下的流态，在桨叶静止的方案3中压力分布基本与达到最大逸速时一致；而在方案1和5中因达到最大逸速后桨叶继续动作，压力分布仍有变化。比较各方案可知：桨叶角度大时，由于流量大，桨叶表面的静压值相对变小。

2.3 尾水管内流场变化规律

图8为方案1、3、5在2 s时、达到最大逸速时、15 s时的尾水管内流态，辅以-1 kPa压力等值面表示尾水管内涡带形状和低压区位置。比较各过程中尾水管进口断面压力变化可见：大桨叶角度下流量增大，进而引发进口水流圆周分速度增大，水流主流向边壁方向移动，周向压力分布不均匀性加大，尾水管中心负压随之上升，形成与水流旋向一致的低压涡带。

在方案1中，关闭桨叶使进口水流流速降低，尾水管内低压区面积随之不断减小，达到最大逸速时，尾水管中心管状涡带直径相比于2 s时已明显缩小，在飞逸趋于稳定的15 s时，涡带形状变为偏向顺时针的柱状涡带，流态较好；而在打开桨叶的方案5中，流量增大使得进口水流周向旋转加剧，引发的轴向回流使低压区由边壁向尾水管中心发展，中心部最大负压可达-3 kPa，为所记录飞逸初始时刻最大负压值的2.86倍，形成的管状涡带直径和影响范围不断增大，顺时针转动趋势越发明显，最终演变为极不稳定的不规则偏心螺旋涡带。

以桨叶静止的方案3为基准比较各方案可知：大桨叶角度下的飞逸过程中尾水管内易出现低压区扩散、水流凌乱等现象，继而引发剧烈压力波动；而关闭桨叶对稳定尾水压力脉动和优化流态有显著作用。

3 结 论

a.采用基于动网格的CFD数值模拟方法实现5种桨叶控制方案下的轴流转桨式水轮机模型的飞逸过渡过程模拟，经数值计算得到的最大逸速与试验值接近，并可以准确获取外特性参数和内流场的变化规律。

b.最大逸速反比于桨叶角度；飞逸过程中增大桨叶角度将加剧外特性参数波动和压力脉动，不利于流场稳定；减小桨叶角度可优化流态并限制机组振动。

c.小桨叶角度对应较大的桨叶压力极值，模拟得到的桨叶压力极大值、极小值可达初值的2.39倍和2.94倍；桨叶角度增大导致负压区面积扩大，压力梯度分布不均。

d.尾水管真空度与桨叶角度正相关，随桨叶角度增大水流周向旋转加剧，引起的轴向回流使涡带由柱状发展为螺旋状，诱发低频高幅的尾水压力脉动。

［1］金建范.轴流转桨式水轮机协联调整改进［J］.水电站机电技术，1993（1）：1-6.（JIN Jianfan.Improvement on coupled relationship for Kaplan turbines［J］.Mechanical Electrical Technique of Hydropower Station，1993（1）：1-6.（in Chinese））

［2］杨琳，陈乃祥.水泵水轮机转轮全特性与蓄能电站过渡过程的相关性分析［J］.清华大学学报（自然科学版），2003，43（10）：1424-1427.（YANG Lin，CHEN Naixiang. Analysis of relation between characteristics of pumpturbine runner and transition of pumped-storage power station［J］.Journal of Tsinghua University（Science and Technology），2003，43（10）：1424-1427.（in Chinese））

［3］潘熙和，程远楚，王杰飞.水泵水轮机全特性曲线处理的研究进展［J］.长江科学院院报，2014，31（12）：117-123.（PAN Xihe，CHENG Yuanchu，WANG Jiefei. Research progress of the processing of pump turbine characteristic curve［J］.Journal of Yangtze River Scientific Research Institute，2014，31（12）：117-123.（in Chinese））

［4］杨建东，曾威，杨威嘉.水泵水轮机飞逸稳定性及其与反S特性曲线的关联［J］.农业机械学报，2015，46（4）：59-64.（YANG Jiandong，ZENGWei，YANGWeijia. Runaway stabilities of pump-turbines and its correlations with S-shaped characteristic curves［J］.Transactions of the Chinese Society for Agricultural Machinery，2015，46（4）：59-64.（in Chinese））

［5］刘延泽，常近时.灯泡贯流式水轮机装置甩负荷过渡过程基于内特性解析理论的数值计算方法［J］.中国农业大学学报，2008，13（1）：89-93.（LIU Yanze，CHANG Jinshi.Numerical method based on internal character for load rejection transient calculation of a bulb turbine installation［J］.Journal of China Agricultural University，2008，13（1）：89-93.（in Chinese））

［6］邵卫云.含导叶不同步装置的水泵水轮机全特性的内特性解析［J］.水力发电学报，2007，26（6）：116-119.（SHAO Weiyun.Complete characteristic of reversible pump turbine with misaligned guide vanes［J］.Journal of Hydroelectric Engineering，2007，26（6）：116-119.（in Chinese））

［7］李小芹，常近时，李长胜.改善水轮机过渡过程品质的导叶异步关闭方式［J］.水力发电学报，2014，33（1）：202-206.（LI Xiaoqin，CHANG Jinshi，LI Changsheng. Guide vane asynchronous closure mode for improving the transient quality of hydroturbine［J］.Journal of Hydroelectric Engineering，2014，33（1）：202-206.（in Chinese））

［8］常近时.轴流转桨式水轮机装置甩负荷过渡过程的合理控制方式［J］.水利学报，1986（7）：29-36.（CHANG Jinshi.A rational controlmode for load-rejection transient process of Kaplan turbine installation［J］.Journal of Hydraulic Engineering，1986（7）：29-36.（in Chinese））

［9］李小芹，常近时，李长胜.抽水蓄能电站中球阀协同导叶关闭的水力瞬变过程控制方式［J］.水利学报，2014，45（2）：220-226.（LI Xiaoqin，CHANG Jinshi，LI Changsheng.Hydraulic transient controlling rules coupled with ball-valve and guide-vane in pumped storage power station［J］.Journal of Hydraulic Engineering，2014，45（2）：220-226.（in Chinese））

［10］杨建东.导叶关闭规律的优化及对水力过渡过程的影响［J］.水力发电学报，1999（2）：78-86.（YANG Jiandong.Optimization ofwicket closure rule and influence on hydraulic transients for water power station［J］.Journal of Hydroelectric Engineering，1999（2）：78-86.（in Chinese））

［11］周大庆，吴玉林，刘树红.轴流式水轮机模型飞逸过程三维湍流数值模拟［J］.水利学报，2010，41（2）：233-238.（ZHOU Daqing，WU Yulin，LIU Shuhong.Threedimensional CFD simulation of the runaway transients of a propeller turbine model［J］.Journal of Hydraulic Engineering，2010，41（2）：233-238.（in Chinese））

［12］夏林生，程永光，张晓曦.灯泡式水轮机飞逸过渡过程3维CFD模拟［J］.四川大学学报（工程科学版），2014，46（5）：35-41.（XIA Linsheng，CHENG Yongguang，ZHANG Xiaoxi.3D CFD simulation of the runaway transients of bulb turbine［J］.Journal of Sichuan University（Engineering Science Edition），2014，46（5）：35-41.（in Chinese））

［13］李金伟，刘树红，周大庆.混流式水轮机飞逸暂态过程的三维非定常湍流数值模拟［J］.水力发电学报，2009，27（1）：178-182.（LI Jinwei，LIU Shuhong，ZHOU Daqing.Three-dimensional unsteady simulation of the runaway transient of the Francis turbine［J］.Journal of Hydroelectric Engineering，2009，27（1）：178-182.（in Chinese））

［14］黄剑峰，张立翔，何士华.混流式水轮机全流道三维定常及非定常流数值模拟［J］.中国电机工程学报，2009，29（2）：87-94.（HUANG Jianfeng，ZHANG Lixiang，HE Shihua.Numerical simulation of 3-D steady and unsteady flows in whole flow passage of a Francis hydro-turbine［J］.Proceedings of the CSEE，2009，29（2）：87-94.（in Chinese））

［15］ZHANG Qiang，YANG Yong.A new simpler rotation/ curvature correction method for Spalart-Allmaras turbulence model［J］.Chinese Journal of Aeronautics，2013，26（2）：326-333.

［16］侯祎华，齐学义，张静.Spalart-Allmaras湍流模型在混流式长短叶片转轮流场计算中的应用［J］.水力发电学报，2010，29（1）：152-157.（HOU Weihua，QI Xueyi，ZHANG Jing.Application of Spalart-Allmaras turbulence model to flow calculation of Francis runner with splitter blades［J］.Journal of Hydroelectric Engineering，2010，29（1）：152-157.（in Chinese））

［17］刘大恺.葛洲坝水电厂ZZ500水轮机过渡过程的动态模型试验报告［R］.南京：河海大学，1992.

Numerical simulation of runaway transients of Kap lan turbinemodel based on blade regulation

ZHOU Daqing1，GUO You1，JIANG Dezheng2
（1.College ofEnergy and Electrical Engineering，HohaiUniversity，Nanjing 211100，China；2.China Yangtze Power Co.，Ltd.，Yichang 443002，China）

In order to study the effect of different laws regarding the opening and closing of the blade on the runaway transient of a Kaplan turbine，the runaway transient of a Kaplan turbine model in five blade controlmodes was simulated using the three-dimensional unsteady numericalmethod.Changes of the rotation speed，flow rate，torques，and pressure fluctuation with time were analyzed，and the pressure distribution on the blade surface and inner flow patterns in the draft tube were examined.The results show that，compared with themaximum runaway speed under the fixed blade conditions，the variation of themaximum runaway speed ranges from-6.6%to 5.0%when the blade angle is changed from-10°to 10°.Opening blades during the runaway transient increases characteristic parameter fluctuations，causes the maximum negative pressure at the central area of the draft tube to reach 2.86 times the initial value，and produces an eccentric spiral vortex rope，which induces strongly low-frequency fluctuation and is adverse to the unit stability.Closing blades can reduce water velocity，alleviate pressure fluctuation，and improve the flow pattern in the draft tube，but reasonable closuremodes should be investigated in order to avoid an excessive increase in the runaway speed.

Kaplan turbine；runaway transient；blade regulation；opening and closure law；three-dimensional unsteady numerical simulation

TK733+.5；TV741

A

1006-7647（2016）04-0013-07

10.3880/j.issn.1006-7647.2016.04.003

2015- 06 19 编辑：熊水斌）

国家自然科学基金（51339005，51106042）

周大庆（1976—），男，副教授，博士，主要从事水利水电动力工程与流体机械研究。E-mail：zhoudaqing@hhu.edu.cn