高粘发酵体系不同搅拌桨的CFD模拟及发酵验证

董淑浩，朱萍，徐晓滢，李莎，姜永翔，徐虹



高粘发酵体系不同搅拌桨的CFD模拟及发酵验证

董淑浩，朱萍，徐晓滢，李莎，姜永翔，徐虹

南京工业大学食品与轻工学院 材料化学工程国家重点实验室，江苏 南京 211816

董淑浩, 朱萍, 徐晓滢, 等. 高粘发酵体系不同搅拌桨的CFD模拟及发酵验证. 生物工程学报, 2015, 31(7): 1099–1107.Dong SH, Zhu P, Xu XY, et al. Computational fluid dynamics simulation of different impeller combinations in high viscosity fermentation and its application. Chin J Biotech, 2015, 31(7): 1099–1107.

搅拌桨是高好氧高黏度微生物发酵实现高效反应必不可少的因素之一，不同搅拌桨组合对发酵过程的影响十分重要。威兰胶是由产碱杆菌在高耗氧高粘度发酵体系下合成的胞外微生物多糖，广泛应用于水泥、石油、油墨、食品等行业中。本研究借助于计算流体力学 (Computational fluid dynamics，CFD) 的方法，以威兰胶发酵液体系为研究体系，研究了6种不同搅拌桨组合在反应器内流体速率分布、剪切速率、和气含率等参数。将模拟效果较好的3种组合用于威兰胶发酵过程。研究表明MB-4-6搅拌桨组合对改善发酵罐内部的溶氧及流场分布效果最明显，威兰胶产量水平提高了13%。同时在该组合下威兰胶的产品粘度得到有效提高。

计算流体力学，非牛顿流体，威兰胶，搅拌桨组合

搅拌桨是生物反应器实现高效反应发酵工业中必不可少的因素之一，尤其高好氧高黏度微生物的发酵过程中，不同的搅拌桨组合对发酵过程的气液传质有着重要影响。随着发酵体系黏度的不断增大，不同桨叶对反应器内各处流体的速度分布矢量、溶氧浓度、剪切速率、营养物浓度等的影响大不相同[1]。威兰胶 (Welan gum) 是由产碱杆菌在高耗氧高粘度发酵体系下合成的胞外微生物多糖，该多糖易溶于水，其水溶液具有独特的假塑性流变学性质和良好的耐高温、耐酸碱性。可作为悬浮剂、增稠剂、稳定剂等广泛应用于水泥、石油、油墨、食品等行业中，因此受到各行业的广泛关注[2-3]。由于威兰胶发酵体系的特殊性，发酵过程总是面临溶氧传递难、物料分布不均匀、能量传质难等问题。因此，研究不同搅拌桨组合对生物反应器的影响对高粘度发酵体系尤为重要。

目前研究者发现不同搅拌桨组合对非牛顿流体发酵过程的传质影响各不相同：李艳等[4]发现抛物线圆盘桨＋翼型轴流桨＋翼型轴流桨耗能最少，罐体溶氧和气体分散性能明显好于其他组合。苏顺开等考察了高黏体系下双斜叶涡轮搅拌桨与圆盘涡轮搅拌桨对发酵体系转速和表观气速的影响，发现后者功率消耗降低35%左右，而氧传质效率提高超过24%[5]。在威兰胶发酵生产中研究者[6]发现高搅拌转速可以促进菌体生长，而中等转速有利于威兰胶的合成和分子量的增大，但关于新型搅拌桨型组合应用于威兰胶发酵的研究却少有报道。因此本文着重考察不同搅拌桨型组合对威兰胶发酵体系的影响。

计算流体力学 (CFD) 是以计算机为基础的一种数值模拟的方法，CFD技术能给出搅拌釜内的总体流动状况和局部流动信息，在搅拌桨的优化中，较传统的半经验半理论的方法有很大的优越性[7-8]。朱向哲等采用CFD对两种轴流式搅拌桨在聚苯乙烯反应器中的传热特性进行了对比研究并得出最佳搅拌桨[9]。Aubin等[10]通过CFD模拟对各种改进的径向流和轴向流搅拌桨在非牛顿流体中的操作进行了数值模拟。蒋啸靖等[11]通过对50 L生物反应器中不同的搅拌桨组合对搅拌流场等参数的CFD模拟来改善发酵罐内部的流场和气体分布。通过CFD技术对搅拌设备的开发研究已取得了显著的成果。

本研究采用CFD技术对不同类型搅拌桨组合进行模拟研究，对宏观流动场速度、剪切速率等参数进行比较分析。在计算模拟得出最优结果的基础上，进行实际发酵验证，筛选得出有利于提高威兰胶产量的最佳搅拌桨组合，从而为威兰胶工业化生产提供理论指导和借鉴。

1 材料与方法

1.1 材料

1.1.1 菌种

产碱杆菌sp. NX-3，本实验室保藏[6]。

1.1.2 培养基

种子培养基 (g/L)：葡萄糖20，酵母膏1，牛肉膏3，K2HPO42，MgSO40.1，CaCO33, pH 7.0，121 ℃灭菌15 min。

发酵培养基 (g/L)：葡萄糖42.19，蛋白胨4.24，K2HPO43，MgSO40.35，CaCO33, pH 7.3，121℃灭菌30 min。

1.1.3 搅拌桨几何尺寸

搅拌桨的参数详见表1。

1.1.4 搅拌桨类型

六种搅拌桨组合及缩写分别为A：RT-6-2，直六叶+直六叶搅拌桨。B：ART-6，斜六叶+直六叶搅拌桨。C：SAT-6，斜六叶+弯曲非对称六叶搅拌。D：MA-4-6，四宽折叶+斜六叶搅拌桨。E：MR-4-6，四宽折叶+直六叶搅拌桨。F：MB-4-6，四宽折叶+弯曲非对称六叶搅拌桨。简图见图1。

1.2 方法

1.2.1培养方法

种子液：接种活化后的斜面菌苔2环于种子培养基中，装液量100 mL/500 mL，30 ℃、200 r/min摇瓶培养16 h。

摇瓶发酵培养：将上述培养好的种子液按10%的接种量接入装有100 mL发酵培养基的500 mL三角瓶中培养，30℃、200 r/ min摇瓶培养72 h。

7.5 L发酵罐培养 (BioFlo110，美国NBS公司，直径T=180 mm，罐高H=320 mm)：温度30 ℃，通气量4.5 L/min，装液量4.5 L，接种量8 %，转速 600 r/min，pH 6.5−7.0。

1.2.2 网格划分

对发酵罐建立3D模型，采用ICEM软件对高溶氧发酵罐划分网格，整体选非结构化的四面体网格，对罐内的罐体、搅拌桨、气体分布器和挡板进行局部加密。对网格非依赖性进行了考察，最终确定计算域最优网格数量为：桨叶区域30万，发酵罐180万。

1.2.3 参数设置

采用Euleri~Eulerian多相流模型模拟，模拟介质为粘度约2.4 Pa∙s的1%威兰胶水溶液，采用标准k-ε二方程湍流模型，气相为空气，通气量为1.1 vvm，压力为常压，温度30 ℃, 装液量为5 L。搅拌转速为600 r/min。罐体高为 320 mm，罐体内径为180 mm，下层搅拌桨距罐体底部10 mm，桨心距为95 mm。

表1 搅拌桨的几何参数

图1 六种不同搅拌组合桨型简图

1.2.4 葡萄糖含量的测定

葡萄糖含量由生物传感分析仪检测 (SBA-40C，山东省科学院生物研究所)。

1.2.5 威兰胶含量测定

称取适量体积的发酵液，蒸馏水稀释， 12 000 r/mim离心30 min，去除菌体细胞沉淀。上清液加入1/5体积的混合液三氯甲烷∶异戊醇=24∶1，/，12 000 r/min离心30 min，去除蛋白质，反复去除2次。上清液用乙醇沉淀威兰胶，60 ℃烘干，称重。

1.2.6 菌体干重测定

取5 mL发酵液并适当稀释，于已预先称重的干燥离心管中，10 000 r/min 离心20 min，弃上清，用去离子水洗涤2次，105 ℃烘箱内烘至恒重，称量即可。

1.2.7 威兰胶粘度的测定

威兰胶1%水溶液，NDJ－1型旋转黏度计 (上海恒平科学仪器有限公司 (4号转子))， 60 r/min，25 ℃测得。

2 结果与讨论

2.1 不同搅拌桨组合CFD模拟分析比较

2.1.1 不同搅拌桨组合剪切速率模拟分析

不同类型搅拌桨对流体会产生不同的剪切速率，从而造成反应器内不同区域的流动状态和混合状态出现显著性差异。而剪切速率对高粘发酵的表观粘度有决定性影响。因此对6种不同搅拌桨组合的剪切速率分析对威兰胶发酵的优化具有重要参考意义。图2是不同搅拌桨组合发酵罐内剪切速率的直观分布，从图中可看出，高粘发酵体系发酵罐液面和底部有部分区域无明显剪切作用，且有效剪切区域主要分布在搅拌桨桨叶附近。四种单桨型的剪切速率区别非常明显。单直六叶搅拌桨 (RT-6) 剪切体积最大，斜六叶搅拌桨 (ST-6) 最小。这主要因为径向桨桨叶垂直轴，有效剪切区域更大，而轴向桨与轴呈一定角度，剪切力部分垂直轴并且轴向分布。图2A中最大剪切速率的RT-6组合剪切速率分布梯度变化非常明显，搅拌桨之间并没有形成很好的剪切传递，区域化明显。图2D为上下均为轴向桨，图2C、F上桨为径向弯曲桨，但由于各个桨的有效剪切面积小没能使上下形成均匀的剪切区间，混合效果不佳。整体均匀性效果最好的是图2E的MR-4-6组合。这说明底部使用四宽叶轴向流搅拌桨有利于提高发酵罐内径向及轴向的搅拌，主体混合更均匀，能减少区域化现象，增加传质效果。

图2 不同桨型液相剪切速率分布

2.1.2 不同搅拌桨组合气含率模拟分析

气含率是反应器中气相所占的体积分率，是评判反应器内的流型、气液接触面积等流体力学参数有重要影响[12]。图3为相同通气量、不同桨型组合下发酵罐内气含率的直观分布情况。从整体来看，六种搅拌桨组合气含率变化非常明显，图3中C、D、F的搅拌效果图远高于A、E组合。对比图3A、B可知底层为轴向搅拌桨时可显著提高整体气含率，这与郝志刚等[13]径向流型搅拌桨对气体的分散作用优于轴流式搅拌桨的结论相反。一般而言搅拌气液混合体系中气体的分散是通过桨叶后方的气穴来分散的，径向搅拌桨后方的气穴中耗散的能量要大于轴流桨，形成的气泡较小，且在槽的底部形成了良好气体再循环，从而使气含率略高于轴流式底桨的各组合桨。模拟出的结果与理论不同，这可能因为轴向搅拌桨拥有较大的叶片覆盖面积率，可以更多地遮挡上升气体，阻止部分气体从搅拌器直接穿过，延长气液接触时间。另外倾斜的六叶可以对分布器进入的气体产生向下的推力，增加了气体的停留时间，从而增大气含率。但这也是从模拟结果提出的解释，为了更确切地进行研究有必要对A、B进行发酵验证。对比图3中B、C搅拌桨组合，在底部轴向桨型相同的情况下，上方为弯六叶搅拌桨时气含率更大。气液混合更加均匀。图3F可以看出上方为弯六叶搅拌桨时相对于图3D两桨之间的混合更加均匀。这主要因为弯六叶搅拌桨的抛物型叶面更符合流体力学，弯曲的叶片上下不对称，阻止了更多的气体直接溢出，与后边上升的气体相混合使得两桨之间气含率明显提高。同时抛物型桨叶叶后不易形成较大的气穴，因此有更多的气体被打散进入流体循环中，从而使得图3C、F的气体分散性能较好。

图3 不同桨型气含率分布

2.1.3 不同搅拌桨速度矢量图

图4可以看出，所有搅拌桨的最大速度出现在叶片边缘，边缘的流体由于挡板和罐内壁的作用产生垂直方向的分量，从而将流体分成几个不同的区间。直六叶搅拌桨 (RT-6) 及弯曲非对称六叶搅拌桨 (BT-6) 在搅拌桨桨叶上下边缘形成了两个涡流，而斜六叶搅拌桨 (ST-6)和四宽折叶搅拌桨 (BP-4) 轴向流的搅拌桨则只形成单一的涡流。这主要因为径向搅拌桨的竖直方向叶片边缘与流体接触后水平撞击罐壁分成上下两个速度方向，并且分别沿着罐壁经过循环后再从桨叶的上边缘进入叶轮。而倾斜的桨叶减少了流体在水平方向上的径向流动，在搅拌桨旋转过程中，流场的变化更为平缓规则。图4A中罐内虽然上下形成4个漩涡，但每个涡流的速度矢量都很小，并且各个漩涡流动性不强[14]。这种组合在微生物多糖的发酵后期，会使得大量液体处于迟滞区，流动非常缓慢，不利于传质，造成死区。图4D两轴向搅拌桨组合的速度矢量值最大，上下混合更均匀。这主要因为这种组合有明显的轴向速度，二者可以使流体循环迅速，将更多底部液体裹挟到整体大循环中，从而缩短混合时间[15]。但这往往会降低气体在罐体内的滞留时间。降低了氧气的利用率。图4C、F的速度矢量整体效果最好，既不会降低气体滞留时间，同时较大的流体速度能提高发酵液传质。对于高粘度发酵搅拌桨所能驱使的液体越多对发酵越有利[16]。

图4 不同桨型液相速度矢量图

2.2 搅拌桨组合发酵参数分析

从模拟实验的3组数据分析得出图3中C、D、F组合的气含率效果较好，E效果最差。图4D组合的速度分布场效果最佳，其死区面积最小。3组数据中C、F两种组合模拟效果最为接近。为了进一步验证模拟结果，本研究对A、B、D、F组合进行发酵实验验证。

四种不同搅拌桨组合对威兰胶发酵过程的影响曲线如图5所示。图5A为不同搅拌桨组合 DCW的变化情况。在发酵进行到50 h时4组搅拌桨DCW趋向稳定，最大DCW分别为 (8.48±0.25)、(8.23±0.22)、(8.06±0.25)和(7.63±0.23) g/L。MA-4-6搅拌桨组合有最大菌体量，说明上下均为轴向搅拌桨有利于提高发酵罐内气含率，特别是在粘度较低的对数生长期时能让整个罐体溶氧更均匀，利于菌体生长。图3D可以很好说明这点。从图5B中不难发现在发酵前30 h时ART-6和RT-6-2组合粘度均高于MA-4-6与MB-4-6，而当发酵的粘度达到3 Pa∙s时后者出现了反超并一直高于前者，两种组合的最大值分别稳定在 (4.75±0.15) Pa∙s和 (4.56±0.15) Pa∙s。说明在发酵后期直六叶桨相对叶片较少的四叶桨、斜六叶桨及弯曲不对称六叶桨表现出更大的剪切稀化，对胶体也产生了严重的机械损伤，降低了发酵液的粘度[17]。图5C表明不同搅拌桨组合对糖消耗也有着显著的影响，ART-6组合对糖的消耗较低 (剩余残糖(19±1.2) g/L)，而同样是轴向与径向组合的MB-4-6组合却有着最高的糖消耗 (剩余残糖(8±1.1) g/L)。这是由于上桨为弯曲非对称六叶搅拌较直六叶桨能产生向下循环的流体。使更多的气体能进行强烈的上下混合，同时抛物线型桨也减少了气穴形成的可能[18]。在转子周围区域能具有较高的气含率，更利于菌体生长和胶的产生，增加了糖耗。图5D为4种搅拌桨组合威兰胶产量的变化情况。MB-4-6组合产量最高，达到了 (27.3±0.6) g/L。这主要因为发酵前期由于四折叶桨采用下压式，可以将上升的气泡再次向下翻动，有利于氧的传质，使得菌体DCW有一个较高的值。而在发酵中后期，此组合上桨为六抛物线桨能形成上下混合又增加了气泡在液体中的滞留时间，搅拌桨能更多地驱动发酵液上下混合使罐体传质更加均匀。而DCW值最高的MA-4-6组合产量并不是最高，这表明，菌体量的增多并不一定会引起威兰胶合成量的增多[19]。这主要因为发酵前期MA-4-6组合有较好的上下气液传质混合效果，促进菌体生长，菌体量累积达到了一个较大的值。而在发酵中后期，特别是粘度升高后，传质水平开始下降，罐内气体径向分布明显会降低，不利于气体在罐内停留。造成溶氧不足，不利于胞外多糖的分泌，从图2D中的气含率分布上也很好地说明了这一点。

图5 不同搅拌桨组合对威兰胶菌体量(A)、粘度(B)、葡萄糖消耗量(C)和产量(D) 的发酵曲线影响

3 结论

不同桨叶组合对高粘度反应器系统内的液体宏观流场分布有较大的差异，径向与轴向搅拌桨组合使得发酵罐内气液混合更均匀。对6种搅拌桨组合进行CFD模拟及发酵实验验证，得出MB-4-6组合 (即下为四宽折叶上为弯曲非对称六叶搅拌桨) 对威兰胶发酵效果最好。在发酵前期该组合搅拌均匀能促进菌体生长，后期该组合桨利用下压式将上升气泡再次下压到桨叶下形成多个循环，增加了气泡在发酵罐内时间。提高了溶氧及罐体内传质[20]。在对优化的搅拌桨进行发酵验证时，MB-4-6组合发酵参数整体效果最佳，相对于直六叶搅拌桨产量提高了13%。CFD模拟结果和实际发酵结果相一致。因此采用CFD模拟对搅拌桨进行优化是可行的。

[1] Kawase, Moo-Young M. Mixing time of bioreactors. J Chem Technol Biot, 1989, 44(1): 63−75.

[2] Xu HM, Zhao L, Wu H, et al. Experimental research on the new foam drilling fluid system. Adv Mater, 2013, 28(57): 781−784.

[3] Prajapati VD, Jani GK, Zala BS, et al. An insight into the emerging exopolysaccharide gellan gum as a novel polymer. Carbohyd Polym, 2013, 92(54): 670−678.

[4] Li Y , Zhan XB, Zheng ZY, et al. Study on oxygen transfer characteristics of combined impeller in non-Newtonian fluid. Food Mchnr, 2008, 24(4): 93−97 (in Chinese).李艳, 詹晓北, 郑志永, 等. 新型搅拌桨组合在非牛顿流体中的传氧性能研究. 食品与机械, 2008, 24(4): 93−97.

[5] Su SK, Ji XY, Zheng ZY, et al. Computational fluid dynamics simulation of gas-liquid mass transfer in xanthan gum solution by novel impeller. Chem Eng, 2010, 38(10): 172−177 (in Chinese).苏顺开, 季新跃, 郑志永, 等. 新型搅拌桨用于黄原胶溶液气液传质的计算流体力学模拟. 化学工程, 2010, 38(10): 172−177.

[6] Li H, Xu H, Li S, et al. Effects of dissolved oxygen and shear stress on the synthesis and molecular weight of welan gum produced fromsp. CGMCC2428. Process Biochem, 2011, 46(11): 1172−1178.

[7] Liu L, Wang M, Du G, et al. Enhanced hyaluronic acid production ofby an intermittent alkaline-stress strategy. Lett Appl Microbiol, 2008, 46(3): 383−388.

[8] Huang WC, Chen SJ, Chen TL. The role of dissolved oxygen and function of agitation in hyaluronic acid fermentation. Biochem Eng J, 2006, 32(3): 239−243.

[9] Zhu XZ, Wang WQ, Ma WT, et alStudy on thermal characteristic of two axial flow impellers by CFD. Petro-Chem Equipment, 2007, 36(4): 20−32 (in Chinese).朱向哲, 王卫强, 马文涛, 等. 两种轴流式搅拌桨传热特性的CFD研究. 石油化工设备, 2007, 36(4): 20−32.

[10] Aubin J, Fletcher DF, Xuereb C. Modeling turbulent flow in stirred tanks with CFD: the influence of the modeling approach, turbulence model and numerical scheme. Exp Therm Fluid Sci, 2004, 28(5): 431−445.

[11] Jiang XJ, Xia JY, Zhao J, et alCFD Simulation of mixing process in stirred bioreactor and its application in fermentation. Chem Bioeng, 2008, 25(7): 54−57 (in Chinese).蒋啸靖, 夏建业, 赵劼, 等. 生物搅拌反应器内混合情况的CFD模拟及在发酵中的应用. 化学与生物工程, 2008, 25(7): 54−57.

[12] Hu FT, Huo WD, Shi SD, et alDetermition methods of hydrodynamic parameters in the loop reactor. Sci Technol Chem Ind, 2007, 15(1): 42−45 (in Chinese).胡发亭, 霍卫东, 史士东, 等. 环流反应器流体力学参数测定技术研究. 化工科技, 2007, 15(1): 42−45.

[13] Hao ZG, Bao YY, Gao JK, et alGas-Liquid dispersion in a multi-impeller stirred tank. J Chem Eng Chin Univ, 2004, 5(18): 547−552 (in Chinese).郝志刚, 包雨云, 高正明, 等. 多层组合桨搅拌槽内气-液分散特性的研究. 高校化学工程学报, 2004, 5(18): 547−552.

[14] Wu GJ, Chen JP, Huang JK, et al. Effect of agitation on oxygen transfer inHY96-2 culture broth and CFD simulation. J East China Univ Sci Technol: Nat Sci Ed, 2009, 35(3): 339−345 (in Chinese).吴高杰, 陈剑佩, 黄建科, 等. 搅拌对多黏类芽孢杆菌发酵氧传递过程的影响及CFD模拟. 华东理工大学学报: 自然科学版, 2009, 35(3): 339−345.

[15] Scargiali F, Busciglio A, Grisafi F, et al. Simplified dynamic pressure method for kLa measurement in aerated bioreactors. Biochem Eng J, 2010, 49(6): 165−172.

[16] Kogan G, Soltes L, Stern R, et al. Hyaluronic acid: a natural biopolymer with a broad range of biomedical and industrial applications. Biotechnol Lett, 2007, 29(1): 17−25.

[17] Gao HJ, Chen J, Du GC, et al. Effect of agitation and mixing on hyaluronic acid produced by. J Chem Ind Eng, 2003, 54(3): 350−356 (in Chinese).高海军, 陈坚, 堵国成, 等. 搅拌与混合对兽疫链球菌发酵生产透明质酸的影响. 化工学报, 2003, 54(3): 350−356.

[18] Liu L, Du G, Chen J, et al. Enhanced hyaluronic acid production by a two-stage culture strategy based on the modeling of batch and fed-batch cultivation of. Bioresource Technol, 2008, 99(17): 8532−8536.

[19] Li S, Xu H, Li H, et al. Optimizing the production of welan gum byNX-3 using statistical experiment design. Afr J Biotechnol, 2010, 9(7): 1024−1030.

[20] Gu XH, Duan XJ, Tan WS, et al. Effect of agitation on hyaluronic acid produced byby using computational fluid dynamics. Chin J Biotech, 2009, 25(11): 1671−1678 (in Chinese).顾小华, 段须杰, 谭文松, 等. 采用计算流体力学技术研究搅拌对兽疫链球菌发酵生产透明质酸的影响. 生物工程学报, 2009, 25(11): 1671−1678.

(本文责编 郝丽芳)

Computational fluid dynamics simulation of different impeller combinations in high viscosity fermentation and its application

Shuhao Dong, Ping Zhu, Xiaoying Xu, Sha Li, Yongxiang Jiang, and Hong Xu

State Key Laboratory of Materials-Oriented Chemical Engineering, College of Food Science and Light Industry, Nanjing Tech University, Nanjing211816, Jiangsu, China

Agitator is one of the essential factors to realize high efficient fermentation for high aerobic and viscous microorganisms, and the influence of different impeller combination on the fermentation process is very important. Welan gum is a microbial exopolysaccharide produced bysp. under high aerobic and high viscos conditions. Computational fluid dynamics (CFD) numerical simulation was used for analyzing the distribution of velocity, shear rate and gas holdup in the welan fermentation reactor under six different impeller combinations. The best three combinations of impellers were applied to the fermentation of welan. By analyzing the fermentation performance, the MB-4-6 combination had better effect on dissolved oxygen and velocity. The content of welan was increased by 13%. Furthermore, the viscosity of production were also increased.

computational fluid dynamics, Non-Newtonian fluids, welan gum, impeller combinations

December 26, 2014;

February 10, 2015

National Basic Research Program of China (973 Program) (No. 2013CB733603), National High Technology Research and Development Program of China (863 Program) (No. 2013AA020301), National Key Technologies Research and Development Program of China (No. 2011BAD23B04), National Nature Science Foundation of China (No. 21106062).

Hong Xu. Tel/Fax: +86-25-58139433; E-mail: xuh@njtech.edu.cn

国家重点基础研究发展计划(973计划) (No. 2013CB733603)，国家高技术研究发展计划(863计划) (No. 2013AA020301)，国家科技攻关计划(No. 2011BAD23B04)，国家自然科学基金(No. 21106062) 资助。