塔中地区构造应力场数值模拟研究

邱登峰，郑孟林，张瑜，张仲培

(中国石化石油勘探开发研究院构造与沉积储层实验室，北京100083)

塔中地区构造应力场数值模拟研究

邱登峰，郑孟林，张瑜，张仲培

(中国石化石油勘探开发研究院构造与沉积储层实验室，北京100083)

构造应力场研究对于油气运聚分析及储层特征描述具有重要的理论和实际意义。本文根据塔中隆起断裂特征及其演化史，结合钻井资料，选择中下奥陶统碳酸盐岩为主要研究对象，应用大型有限元分析软件ANSYS，在处理断层问题时运用非连续的接触算法，模拟了中奥陶世末、中泥盆世末、三叠纪末、新近纪末的构造应力场，研究了该区不同时期的古应力分布特征。研究表明，最大主压应力受东部车尔臣－星星峡走滑断裂影响显著，在隆起区沿断裂带呈条带状低值分布，在不同时期塔中Ⅰ号断裂带无一例外地位于最大主压应力的最低区，反映了塔中断裂带尤其是塔中Ⅰ号断裂带是油气运移的有利指向区。最小主压应力沿断裂带呈条带状高值分布，在断裂的上盘、断裂走向发生变化和断裂的倾末端是张应力高值区，为张裂缝发育的重要构造部位。最小主压应力和最大剪应力受车尔臣－星星峡走滑断裂的影响明显减弱。

塔中;构造应力场;有限元;数值模拟;非连续

构造应力场就是在一个空间范围内构造应力的分布(万天丰，1988)，各种构造迹象，包括矿物颗粒在三度空间排列方位的规律性，是岩石在构造应力场中应力作用的反映(李四光，1973)。它不仅控制着裂缝的形成和演化(Treagus，1981;单家增等，2004)，也是油气运移的重要动力之一(Perrodon，1992;李明诚，1994;王毅等，2005)。因此，开展古构造应力场研究，尤其是古构造应力场的定量研究，对进行含油气盆地储集层特征描述及油气运聚分析具有重要的理论和实际意义(孙晓庆，2008)。构造应力场研究的基本方法为利用岩石中已存在的构造形变特征来反推形变作用发生前后的构造应力状态，除地质研究方法之外，还有模拟研究方法，即包括各种物理模拟与数值模拟(万天丰，1988)。有限单元法是构造应力场数值模拟的重要方法之一。目前在构造应力场数值模拟及其与油气运聚关系及裂缝预测方面的研究取得了一批非常有价值的成果(蒋有录等，2005;杨奎锋等，2006;刘翠荣，2002;刘聪等，2008;佟彦明，2007;张胜利等，2007;桑广森等，2010)。不过，由于计算条件和理论方法的限制，以往的多数工作往往采用连续变形算法，在处理断层问题时，常采用软弱带来模拟断层的物理性质。这对断裂构造的模拟有很大的不足(郑勇等，2007;雷显权和陈运平，2011)，而断裂构造是造成地壳岩体中应力发生复杂变化的主要因素之一(苏生瑞，2002)。不连续体概念和不连续变形理论的逐渐成熟使得利用数值模拟方法研究断层真实动力学状态成为可能(石根华，1993;郑勇等，2007)。本文以塔中隆起为研究对象，根据塔中隆起断裂特征及其演化史，结合钻井资料，应用非连续的接触模型模拟断层，对中下奥陶统不同时期的构造应力场进行了有限元数值模拟研究，并对研究区内构造应力场与油气运聚关系进行了初步探讨。

1 塔中地区的构造地质概况

塔中地区位于塔里木盆地中央隆起的中段，东界呈过渡形式与古城墟隆起西南部相连，南界以逆冲断裂带形式紧邻塘古孜巴斯坳陷，西以吐木休克断裂、巴东断裂为界与巴楚隆起分隔，西北与阿瓦提凹陷呈NW向倾斜过渡，北界以斜坡形式与满加尔坳陷相邻(贾承造，1997;刘克奇和金之钧，2004)，面积约2.75万km2。塔中隆起为一个比较完整的由多个次级构造带组成的大型断背斜构造，形成于中奥陶世末，定型于中泥盆世末期，石炭系及以上地层构成平缓单斜，表明后期被埋藏。多期构造演化和复杂的演化过程导致了该区存在多期地层不整合(李明杰等，2004)(图1)，其中中奥陶世末、奥陶纪末、中泥盆世末、二叠纪末、三叠纪末以及新生代的构造事件对该区的构造变形、地层不整合的形成以及油气的聚集与调整(或改造)等都具有重要的影响。

塔中隆起是塔里木盆地重要的油气聚集区带，主要产油气层系是奥陶系顶面风化壳和石炭系下部海相砂岩层(贾承造，1997)。在中奥陶世末，西昆仑库地洋向塔里木板块俯冲，区域应力场由拉伸状态变为侧向挤压状态，塔中隆起形成，在隆起顶部中下奥陶统遭受剥蚀，灰岩出露地表遭受长期的淋滤、风化形成塔中地区第一套有利的储集层系(张光亚，2000;刘克奇和金之钧，2004;刘训和王永，1995;汤良杰，1997)。中泥盆世末期的早海西运动使塔中隆起构造格局基本定型，以后进入构造相对稳定的发展时期(刘克奇和金之钧，2004)。油气藏主要形成时期为加里东晚期－海西早期、海西晚期－印支期、喜马拉雅期(康志宏等，2001;庞宏等，2010)。

图1 塔中隆起NS-520地震地质解释剖面(据李明杰等，2004略改)Fig.1 Interpreted geological seismic section of line NS-520 in the Tazhong uplift

本文选择塔中地区主要储集层系之一、构造改造强烈的中下奥陶统碳酸盐岩层为二维模拟对象，在认为水平应力为盆地构造变形的主导应力场的情况下，忽略上覆地层压力的影响，模拟其在中奥陶世末、中泥盆世末、三叠纪末和新近纪末等关键构造期的构造应力分布，并研究该区的油气聚集特征。

2 构造应力场的数值模拟

有限元方法是一种比较成熟的计算应力分布的数值方法，如果能提供合理的地质模型和边界条件，能比较真实地反映应力的分布状态。它是将地质体离散成有限个连续的单元，每个单元内赋予其实际的岩石力学参数，根据边界条件，用数值方法计算各单元内应力分布。二维有限元方法在计算中假设其垂向应力与应变忽略不计(汤良杰等，2007)。

2.1 地质模型的建立

数值模拟时以各时期断裂分布图作为模型的基本框架，再加上各期地层的岩相界线分区便构成了模拟的地质模型(张胜利等，2007)。

塔中隆起断裂发育(图2)，可划分为两期。中奥陶世末，塔中隆起南部的巴东南－塘北断裂带还没有形成，由于库地洋向塔里木板块下的的消减俯冲与挤压，台盆区整体处于北东向挤压构造背景，形成了塔中Ⅰ号断裂带、塔中Ⅱ号断裂带和塔中22井南断裂带等走向北西西－北西的逆冲断裂构造。巴东南－塘北断裂带、塔中1－8井断裂带、塔中3井断裂带等为走向北东、倾向北西的断裂系，形成于奥陶纪末期。因此，在进行应力场模拟时，不同时期进行不同的建模，中奥陶世末的构造应力场模拟仅考虑塔中隆起塔中Ⅰ号断裂带、塔中Ⅱ号断裂带、塔中22井南断裂带等北西西－北西向断裂。奥陶纪以后的各期模拟则考虑两期断裂。车尔臣－星星峡断裂带是该区东侧北东东向展布的多期活动的大型断裂带，在各期的模拟中都认为其存在，并主要发生走滑构造作用。

图2 塔中隆起中下奥陶统顶面断裂分布图Fig.2 Distribution of faults in the Middle-Lower Ordovician surface of the Tazhong uplift

根据钻井资料，模拟区早中奥陶世早期以白云岩、灰质白云岩为主，晚期以白云质灰岩、灰岩为主(王国司，2002)，岩性横向变化不大，总体上模拟层位为海相沉积的碳酸盐岩，对岩相地质单元进行均一化处理，选取统一的岩石力学参数。

2.2 数学模型的建立及参数的选取

模拟时选取ANSYS中具有二次形式的PLANE2单元模拟岩性均质连续体，通过岩心取样的岩石力学参数实验(刘聪等，2008)，我们得到有限元模拟的岩石物理参数为:弹性模量E=57.9 GPa，泊松比 μ =0.229。

本次模拟采用非连续模型的接触单元来模拟断层。连续模型与非连续接触模型的区别在于后者考虑了断层带与围岩之间不连续运动的接触关系，在结合处通过接触单元把二者联系起来，二者在接触位置可以发生错动;而对于连续模型，断层带与围岩在结合处具有公共节点，二者之间不会发生错动，研究结果表明，非连续接触模型比连续模型能够更好地模拟断层(雷显权和陈运平，2011)。应用非连续性接触模型模拟青藏高原断层活动的结果显示，非连续模型的运动场分布与GPS观测结果吻合程度大大高于连续体模型结果(郑勇等，2007)。

本文模型采用面面接触方式来处理断层和围岩的接触关系。将断层带与围岩的接触边界相互视为对方的接触面和目标面，从而组成接触对，然后，通过接触力学方法分析断层的运动状态和力学性质(ANSYS Inc，2001;郑勇等，2007;雷显权和陈运平，2011)。接触分析需要计算垂直于目标面的法向接触应力和平行于目标面的切向接触应力，采用Pinball算法进行计算。接触面和目标面之间的间隙(穿透量)用g来表示，当接触面穿过目标面时，就发生接触穿透，规定此时g为负值。罚函数法计算法向接触应力的公式如下:

其中fn为法向接触应力，Kn为法向接触刚度。

切向接触应力是由接触面在目标面上移动所产生的摩擦力引起的，若接触面沿目标面的切向位移的弹性分量为，则切向接触应力为:

其中ft为切向接触应力，Kt为切向接触刚度，F为静态/动态摩擦因子，τ为库伦滑动摩擦力，k为摩擦系数。

本次模拟采用扩展拉格郎日乘子法为接触算法。扩展拉格朗日乘子法通过罚函数的反复迭代得到精确的拉格朗日因子，与罚函数法相比，扩展拉格朗日乘子法常能得到更好的处理结果，对罚函数中参数的敏感性也更低(ANSYS Inc，2001)。

模拟中选取与PLANE2单元相适应的CONTA172和TARGE169面面高阶接触单元形成接触对来模拟断层。本次模拟处理断裂带8条，分别为塔中I号断裂带、塔中10井断裂带、塔中Ⅱ号断裂带、塔中22井南断裂带、巴东南－塘北断裂带、塔中1-8井断裂带、塔中3井断裂带和塔中5-48井断裂带，共生成8个接触对，并根据断层性质不同，赋予了不同的接触参数(表1)。扩展拉格郎日乘子法需设置四个主要参数表征接触单元，即:法向刚度Kn、切向刚度Kt、摩擦系数k及最大允许穿透值gmax。因区内多条断层呈现出挤压－走滑性质，为表现具走滑性质断层的切向滑动特征，相对于逆断层，模拟时通过减小摩擦系数、减小法向刚度、增大切向刚度的方式来体现。最大允许穿透值gmax的取值兼顾计算精度和耗时。本文参考地球物理学和地震学等资料，计算了多种参数情况下的模型结果，通过比较分析，采用如下表所示的断裂接触对参数设置。

表1 主要断裂接触对参数设置Table 1 Contact pair parameters of the main faults

2.3 网格划分及边界条件加载

在建立有限元数值模型后，对数值模型进行了三角形自动网格划分。在模型块体内部，网格尺度较大，三角形单元边长为10 km左右。在外框及接触单元边界，为保证计算精度，网格较为稠密，三角形单元边长甚至小于1 km。从块体内部到外框及接触单元边界，网格尺寸均匀过渡，无畸形网格，共划分出144545个单元，292711个结点。图3所示为新近纪末期网格划分后加载边界条件的有限元模型。

图3 新近纪末期加载边界条件的数学模型Fig.3 Mathematic model with boundary conditions of the Late Neogene

表2 塔中隆起中下奥陶统各时期主应力大小与方向①Table 2 Value and direction of the principal compressive stress in the Middle-Lower Ordovician of the Tazhong uplift

研究表明，塔里木盆地南缘构造活动对塔中古隆起构造的形成与发展起了重要的控制作用，寒武纪－奥陶纪，库地洋向塔里木板块俯冲，并于志留纪与塔里木大陆碰撞，形成库地缝合带，二叠纪古特提斯洋沿康西瓦－玛沁向北缘俯冲产生昆仑晚古生代岛弧，并于二叠纪末－三叠纪初受到甜水海地体碰撞，始新世印度板块与欧亚板块碰撞，并持续向北推挤(贾承造，1997)。受一系列塔里木盆地南缘构造活动的影响，模拟时外力载荷从模型的南部边界施加。边界外力的大小和方向主要根据前人研究塔里木盆地所得出的应力场(丁原辰，1996;王喜双等，1997;曾联波，2005)。在给模型南部边界施加不同时期载荷时(表2)，保持模型内部结构及车尔臣－星星峡东南边界相对位置不变，根据主应力方向，将南部边界旋转一定的角度，并从南部边界节点上施加垂直于南部边界的对应于模拟时期的最大主应力载荷(图4)。塔中隆起北部为满加尔凹陷，构造变形相对简单，地层平缓，模型北部设置为法向、切向位移全约束;模型东南以车尔臣－星星峡走滑断裂为界，设定为法向位移全约束，切向自由边界，以此模拟断裂的走滑特性。塔中隆起以西为巴楚隆起，巴楚隆起有着与塔中隆起近乎相同的变形程度，故模型西部设置为自由边界。

3 模拟结果分析与讨论

对上述加载不同边界条件的数学模型在ANSYS中求解，在后处理模块中可得到四个不同时期(中奥陶世末、中泥盆世末、三叠纪末、新近纪末)最大主压应力、最小主压应力、最大剪应力较清晰的应力分布云图，其中中奥陶世末的模拟结果见图5、6、7。

图4 中奥陶世末边界条件加载示意图Fig.4 Sketch map showing the assumed boundary conditions in the end of Middle Ordovician

模拟结果显示，构造应力场沿断裂呈条带状分布，在断裂末端、转折部位及断裂交叉处易形成应力场富集。

中奥陶世末最大主压应力模拟结果显示，塔中地区东、南、西三个方向上形成环绕塔中隆起的最大主压应力高值区，最大主压应力高值带呈北北东向或北东东向展布，隆起内部相对较低，塔中Ⅰ号、Ⅱ号、10井、22井南断裂带附近更低，形成沿断裂带北西－北西西走向低值分布条带，塔中Ⅰ号断裂带位于最大主压应力最低区。断裂带的西部倾末端都为最大主压应力集中点，且上盘比下盘最大主应力高。塔中隆起的东段与车尔臣－星星峡断裂相交部位最大主压应力分布集中，呈北东东向带状展布。断裂的走向发生变化的部位，如塔中5井区、塔中9井区最大主压应力相对较高。中泥盆世末－新近纪末，塔中隆起的断裂构造虽然没有活动或活动很弱，但在不同方向的挤压应力场作用下，其最大主压应力分布类似，东部以北东走向应力集中带为特征，隆起区为北西向低应力集中区。

图5 中奥陶世末构造应力场最大主压应力模拟结果Fig.5 Simulation results of the maximum compressive principal stress in the end of Middle Ordovician

图6 中奥陶世末构造应力场最小主压应力模拟结果Fig.6 Simulation results of the minimum compressive principal stress in the end of Middle Ordovician

图7 中奥陶世末构造应力场最大剪应力模拟结果Fig.7 Simulation results of the maximum shear stress in the end of Middle Ordovician

不同时期的最小主压应力等值线与断裂的展布方向一致，沿断裂带呈条带状高值分布，在断裂的上盘是张应力高值区，易形成张裂缝，在断裂走向发生变化和断裂的倾末端也是张应力高值区，是张裂缝发育的重要构造部位。东部的走滑断裂对最小主压应力的影响与最大主压应力相比明显减弱。

最大剪应力的分布与最大主压应力和最小主压应力的分布呈现类似规律，断裂的端点部位和走向发生变化的部位是剪应力高值分布区。在塔中地区西部形成北北东向剪应力高值区，在走滑明显的断裂附近剪应力偏高。中奥陶世末断裂的上盘最大剪应力高，下盘相对较低。中泥盆世末、三叠纪末和新近纪末的最大剪应力分布与中奥陶世有些差异，如在Ⅰ号断裂带下盘的最大剪应力较上盘高。

车尔臣走滑断裂对最大主压应力的影响作用最强，对最小主压应力影响作用最弱。中奥陶世末与后三期模拟结果相比，因为还没形成北东向断裂，在塔中22井南断裂以南，最大主压应力高于后三期，最小主压应力低于后三期。随不同时期边界外力方向的变化，各主应力高值区略有变化;由中奥陶世末到新近纪末边界外力大小递减，模拟对应期次的构造应力场大小亦随之递减。

在模型受边界挤压应力作用下，分析模型内部最大主压应力分布图，可以发现最大主压应力沿断裂带呈条带状低值分布，而且在不同时期塔中Ⅰ号断裂带都无一例外地处于最大主压应力最低区。该模拟结果与对该区进行的流体势能场分析结果一致。结合构造应力场与流体势能场分析结果，发现流体高势能区大多位于高压应力区，而低势能区常常与低压应力区对应，说明构造应力场对流体势能场具有直接的影响作用，高压应力和低压应力分别是形成高势能区和低势能区的必要条件。而油气运聚是由高压应力区运移到低压应力区并聚集成藏(王喜双等，1997)。因此，塔中隆起的断裂发育区始终是油气运移的有利指向区，塔中Ⅰ号断裂带始终处于最大主压应力最低区可能也是塔中Ⅰ号断裂带油气富集的原因之一。

致谢:在研究过程中曾得到中国铁道科学研究院铁道科学技术研究发展中心方兴助理研究员的帮助和支持，在此表示感谢!同时感谢审稿专家和编辑部老师提出的宝贵意见和建议!

丁原辰.1996.塔里木盆地北部油田古应力的AE法测量.地质力学学报，2(2):18－25.

贾承造.1997.中国塔里木盆地构造特征与油气.北京:石油工业出版社:2，99，277－279.

蒋有录，张乐，鲁雪松，李景明，王红军.2005.基于ANSYS的应力场模拟在库车坳陷克拉苏地区的初步应用.天然气工业，25(4):42－45.

康志宏，魏历灵，康艳芳.2001.塔里木盆地油气藏形成期分析.新疆石油地质，22(6):462－464.

雷显权，陈运平.2011.应用非连续模型研究断层对地壳应力的影响.中南大学学报(自然科学版)，42(8):2379－2386.

李明诚.1994.油气运移研究的现状和进展.石油勘探与开发，21(2):1－13.

李明杰，郑孟林，冯朝荣，张军勇.2004.塔中低凸起的结构特征及其演化.西安石油大学学报(自然科学版)，19(4):43－47.

李四光.1973.地质力学概论.北京:科学出版社:166.

刘聪，黄晓波，樊太亮，王增香，曾清波.2008.塔中地区奥陶系现今构造应力场模拟及裂缝预测.新疆石油地质，29(4):475－477.

刘翠荣.2002.川西坳陷喜山期构造应力场数值模拟及裂缝预测.天然气工业，22(3):10－13.

刘克奇，金之钧.2004.塔里木盆地塔中低凸起奥陶纪油气成藏体系.地球科学——中国地质大学学报，29(4):489－494.

刘训，王永.1995.塔里木板块及其周缘地区有关的构造运动简析.地球学报，16(3):246－260.

庞宏，庞雄奇，石秀平，向才富.2010.调整改造作用对塔中油气藏的影响.西南石油大学学报(自然科学版)，32(1):33－39.

桑广森，夏斌，张胜利，蔡周荣，梁正中，万志峰.2010.松辽盆地徐家围子三维构造应力场数值模拟研究.大地构造与成矿学，34(2):196－203.

单家增，张占文，陈绍生，许坤.2004.大民屯凹陷安福屯潜山带古构造应力场与裂缝发育特征的光弹物理模拟实验研究.石油勘探与开发，31(4):15－18.

石根华.1993.块体系统不连续变形数值分析新方法.北京:科学出版社:50－80.

苏生瑞.2002.断裂构造对地应力场的影响及其工程意义.岩石力学与工程学报，21(2):296－296.

孙晓庆.2008.古构造应力场有限元数值模拟的应用及展望.断块油气田，15(3):31－34.

汤良杰.1997.略论塔里木盆地主要构造运动.石油实验地质，19(2):108－114.

汤良杰，贾承造，王英民，曲国胜，曾联波，谭成轩，陈书平，刘豪.2007.塔里木叠合盆地构造解析和应力场分析//金之钧，王清晨.中国典型叠合盆地油气形成富集与分布预测丛书.北京:科学出版社:9.

佟彦明.2007.胶莱盆地莱阳期古构造应力场分析及模拟.大庆石油地质与开发，26(1):6－9.

万天丰.1988.古构造应力场.北京:地质出版社:2－4.

王国司.2002.塔里木盆地塔中地区奥陶系地层沉积特征.贵州地质，19(3):179－183.

王喜双，李晋超，王绍民，宋惠珍.1997.塔里木盆地构造应力场与油气聚集.石油学报，18(1):23－28.

王毅，宋岩，单家增.2005.构造应力在油气运聚成藏过程中的作用.石油与天然气地质，26(5):563－571.

杨奎锋，杨坤光，曾佐勋.2006.黄骅坳陷中区明下段末期构造应力场三维数值模拟.矿物岩石地球化学通报，25(1):82－86.

曾联波.2005.塔里木盆地库车山前构造带地应力分布特征.石油勘探与开发，32(3):59－60.

张光亚.2000.塔里木古生代克拉通盆地形成演化与油气.北京:地质出版社:53－57.

张胜利，夏斌，胡振华，张宴华.2007.丽水－椒江凹陷新生代构造应力场数值模拟与油气运聚关系探讨.大地构造与成矿学，31(2):180－185.

郑勇，陈颙，傅容珊，薛霆颙.2007.应用非连续性模型模拟断层活动对青藏高原应力应变场的影响.地球物理学报，50(5):1398－1408.

ANSYS Inc.2001.Theory Release 5.7.Canonsburg PA:ANSYS Inc:668－669，928.

Perrodon A.1992.Petroleum systems:Models and applications.Journal of Petroleum Geology，15(2):319－325.

Treagus S H.1981.A theory of stress and strain variations in viscous layers，and its geological implications.Tectonophysics，72(1－2):75－103.

Numerical Simulation of the Tectonic Stress Field in the Tazhong Area

QIU Dengfeng，ZHENG Menglin，ZHANG Yu and ZHANG Zhongpei
(Laboratory of Structural and Sedimentological Reservoir Geology，SINOPEC Petroleum Exploration＆Production Research Institute，Beijing100083，China)

The research of tectonic stress field is theoretical important and practical significant for petroleum migration and accumulation，reservoir characteristics description and crack formation analyses.The widely-used finite element analysis software(ANSYS)was applied to simulate tectonic stress field of the end of Middle Ordovician，Middle Devonian，Triassic and Neogene，based on the characteristics and evolution history of faults combined with well drilling data in the Tazhong uplift，with emphasis on the Middle-Lower Ordovician limestone.When faults are concerned，the discontinuous contact algorithm was selected to simulate the kinematic and dynamic characteristics of faults，which was also preferred by previous researches，and indeed discontinuous model is a better choice for fault modeling.On the base of structural analysis，the faults were divided into two stages，and two geological models were established.In the four tectonic events，different boundary conditions were imposed.The research indicates that the maximum compressive principal stress，which was affected markedly by the eastern strike-slip fault named Cherchen-xingxingxia，presents strip-shaped and low-value distribution along faults zone in uplift area.Tazhong 1stfault always had the minimum value of the maximum compressive principal stress in all of the tectonic stages，which reflects the Tazhong fault zone-especially the 1stfault-was a well-directed area for hydrocarbon migration.The minimum compressive principal stress was banded in high-value along faults zone.The hanging wall of fault was the high-value area of tensile stress and easy to form tensile cracks.The turn and the end of the faults were also the high-value area of tensile stress and important structure position of tensile cracks development.The effect on the minimum compressive principal stress and maximum shear stress was weakened from Cherchen-xingxingxia strike-slip fault.

Tazhong;tectonic stress field;finite element;numerical simulation;discontinuous

P554;TE121.2

A

1001-1552(2012)02-0168-008

2011－08－15;改回日期:2011－11－07

项目资助:国家科技重大专项“海相碳酸盐岩油气资源潜力、富集规律与战略选区”(2008ZX05005-001-001)、中国石化科技部“塔里木台盆区构造演化对油气成藏控制研究”(P06019)项目资助。

邱登峰(1983－)，男，助理工程师，构造地质及油气地质。Email:qiudf.syky@sinopec.com