地球物理学进展  2016, Vol. 31 Issue (4): 1783-1788   PDF    
CPU-GPU异构并行计算在数字岩心线弹性静力学有限元模拟中的应用
朱伟1, 於文辉2     
1. 油气资源与勘探技术教育部重点实验室, 地球物理与石油资源学院, 长江大学, 武汉 430100
2. 地球物理与空间信息学院, 中国地质大学, 武汉 430074
摘要: 线弹性静力学有限元模拟是计算数字岩心弹性模量的有效方法之一,已被用于研究数字岩心的弹性模量与其微观结构、物质组成之间的关系.若要形成解决实际问题的能力必须计算足够多的数字岩心样本,而有限元模拟的计算量较大,因此并行计算对于该方法的成功应用非常重要.本文将数字岩心线弹性静力学有限元模拟算法分解为在CPU和GPU上执行的两个部分,由CPU负责协调控制,GPU负责大规模数值计算,实现CPU-GPU异构并行计算,获得计算效率提升.采用该并行算法计算孔隙数字岩心、裂缝数字岩心和裂缝—孔隙数字岩心的弹性模量,得到的弹性模量—孔隙度关系符合一般的岩石物理规律.CPU-GPU异构并行的线弹性静力学有限元模拟能够迅速计算大量数字岩心的弹性模量,提供相当于物理实验的“观测数据”,对岩石物理学研究具有重要的意义.
关键词数字岩心     弹性模量     GPU     并行计算    
Application of CPU-GPU heterogeneous parallelled computation in linear elastic statics finite element simulation on digital rock
ZHU Wei1 , YU Wen-hui2     
1. Key Laboratory of Exploration Technologies for Oil and Gas Resources (Yangtze University)-Ministry of Education, Geophysics and Oil Resource Institute, Yangtze University, Wuhan 430100, China
2. Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China
Abstract: Linear elastic statics finite element simulation is one of the effective methodologies to calculate elastic moduli of digital rock. It has been applied to study relationships between elastic moduli, microstructure and material composition of digital rocks. Parallelled computation is important to linear elastic statics finite element simulation, because computation amount of finite element simulation is large and a certain amount of digital rocks needed to resolve real problems. The linear elastic statics finite element simulation is carried out by CPU-GPU heterogeneous parallelled computation. One large part of code is executed on GPU and the other part of code is executed on CPU. Logical controlling code is executed on CPU, and the most part of numerical computation is on GPU. The computational efficiency increases compared to serial computation. Effective elastic moduli of porous digitial rocks, cracked digital rocks and cracked porous digital rocks are calculated. Relationship between effective elastic moduli and porosity is consistent with normal rock physics law. With the CPU-GPU heterogeneous parallelled computation we can calculate various macro elastic moduli of many digital rocks quickly, obtain “observation data” which is equivalent to laboratory experiments. CPU-GPU heterogeneous parallelled computation makes a great signifance to rock physics research.
Key words: digital rock     elastic moduli     GPU     parallelled computation    
0 引 言

数字岩石物理在数字岩心中模拟各种物理场并计算宏观物理参数,是联系岩石微观组构、微观物理场与宏观物理参数的桥梁.数字岩石物理至少具有以下四个特点:①数字岩心的微观组构是可控的;②微观物理场及其演化过程是可见的;③微观物理场与宏观物理参数间的关系是可定量表征的;④不同宏观物理参数的计算模型是可以相同的.数字岩石物理可以与物理实验、测井分析和岩石物理理论一起组成一个包含一系列理论、方法与技术的完整体系,对岩石物理学研究具有重要价值.

数字岩心线弹性静力学有限元模拟(下文简称为LESS)是数字岩石物理的重要内容之一,其研究始于20世纪90年代,Garboczi和Day(1995)Garboczi(1998)关于复合材料力学性质的模拟是这方面较早的工作.基于三维数字岩心的LESS计算量较大,Bohn和Garboczi(2003)发展了基于集群的并行算法,提高模拟效率.Arns等(2002)张晋言和孙建孟(2012)姜黎明等(2012)孙建孟等(2014)利用LESS研究孔隙流体性质、饱和度变化对数字岩心弹性模量的影响,数值结果与实验测量值、Gassmann理论预测值等一致.Knackstedt等(2003)Grechka等(2006)刘学锋(2010)Zhang和ToksÖz(2012)赵建鹏等(2014)利用LESS研究岩石微观结构变化对数字岩心弹性模量的影响.Knackstedt等(2003)指出在数字岩心建立过程中对孔隙结构和骨架的控制可以减少干扰数据分析的因素.Andr等(2013)将LESS与其他两种弹性模量模拟方法进行了比较,发现三种方法的结果差异较小.刘向君等(2014)联合使用Avizo数字岩心处理软件和Comsol有限元软件计算数字岩心的弹性模量,丰富了现有数字岩石物理的研究手段.朱伟和单蕊(2014)认为数字岩石物理存在间接(理论研究)和直接(实际应用)两种研究路线.从实际应用的角度,Dvorkin等(2011)认为岩石的非均匀性尺度很大,毫米级数字岩心的模拟值与厘米级真实岩心的实测值不能直接对比,但在一些情况下宏观物理参数之间的相关性在较大的尺度范围内具有不变性.

综述前人研究,本文得出:①利用LESS计算数字岩心的弹性模量具有广泛的应用前景,可有效补充岩石物理学研究的方法与技术;②若要形成解决问题的能力必须计算足够多的数字岩心样本,以形成宏观物理参数的变化规律.因此,并行计算对于LESS的成功应用特别重要,但传统的并行计算方式成本较大,难以被广泛应用.图形处理单元(GPU)是实现低成本大规模并行计算的有效途径.利用GPU并行计算提高效率在地球物理学研究中已经存在若干成功的案列(陈召曦等,2012刘守伟等,2013唐祥德等,2015).因此,本文将LESS采用CPU-GPU异构并行实现,计算孔隙数字岩心、裂缝数字岩心和裂缝-孔隙数字岩心的弹性模量,分析弹性模量与孔隙度的关系,讨论CPU-GPU异构并行计算对于岩石物理学研究的意义.

1 数字岩心线弹性静力学有限元模拟的基本原理

数字岩心的基本构成单元是像素.在数字岩心中每一个像素都被赋予一个整数值——像数值,以区分像素所归属的物质类型.根据弹性力学理论,可利用有限元方法计算数字岩心的弹性模量,其基本原理是在数字岩心的表面加载应变(位移)边界条件,计算初始弹性势能,通过共轭梯度法求解最小弹性势能对应的应力-应变状态,计算平均应变张量和平均应力张量,反推平均弹性模量,即数字岩心的弹性模量.

弹性体的弹性势能En与其弹性模量和所受应变有关,其积分形式为

(1)

式中,εα,εβ表示弹性体V内部微小单元的应变,Cαβ表示其弹性模量.

若以像素为有限元单元体,则单元体的弹性势能可以表示为

(2)

式中,u表示单元结点位移向量,D=NTLTCLN表示单元刚度矩阵,N表示形函数矩阵,L表示形函数空间偏导算子矩阵,C表示介质刚度矩阵.若在单元上采用Simpson数值积分,单元体弹性势能可写为求和形式为

(3)

式中,Dxyz表示积分点刚度矩阵.

考虑周期性边界条件,边界单元体的势能函数具有以下形式

(4)

式中,向量δ与单元所处位置有关.数字岩心的总体势能E由各个单元势能方程总装而成,公式为

(5)

式中,U是所有结点位移向量,A由单元刚度矩阵D根据结点位置总装而成,BδTD总装而成,H是常数.

在给定的边界条件下,处于平衡状态的数字岩心的弹性势能取极值,公式为

(6)

通常,A为大型稀疏矩阵.

数字岩心的弹性模量由平均应变张量和平均应力张量反算得到.当数字岩心表现为弹性各向同性时,其体积模量K和剪切模量G可以表示为

(7)

式中,σ11,σ22,σ33,σ12,σ23,σ31表示平均应力张量的分量,ε11,ε22,ε33,ε12,ε23,ε31表示平均应变张量的分量.

关于LESS的更详细说明,读者可以参考Garboczi(1998)刘学锋(2010)的工作.

2 CPU-GPU异构并行计算

CPU和GPU在结构和功能上采用了不同的设计理念.CPU芯片的设计目标是缩短单个线程的执行时间,而GPU芯片的设计目标是提高线程的总吞吐量(Kirk and Hwu,2010).因此,CPU芯片的主要部分是逻辑控制单元和缓存,GPU芯片的主要部分是数值运算单元. GPU的并行方式是单程序多线程并行,因此GPU并行计算特别适合矩阵或数组运算.一块GPU芯片上的计算核心一般只有数十至数百个,而一个内核函数生成的线程数量可达成千上万个,因此GPU中执行相同指令段的线程也有可能分批执行.

线程组织对GPU并行计算非常重要.一方面,流处理器上能够分配的线程块的数量和大小受到硬件资源的限制.另一方面,线程的组织会影响显存数据的读取与存储速度,线程块中的线程数量为32的倍数比较有利. 采用CPU-GPU异构并行计算利用CPU和GPU各自的优势,提高算法执行效率.本文将有限元模拟中需要大量循环迭代的环节用GPU并行实现,其他执行逻辑控制和只需少量循环迭代的环节仍然由CPU串行执行.本文模拟所用微机的CPU主频为3.2 GHz,内存容量为8 GB,GPU卡包含384个计算核心,显存为4 GB.通过计算一个大小为1503像素的数字岩心,比较共轭梯度法迭代过程的时间,GPU并行相比串行的执行效率可以提升8倍以上.

3 并行计算实例

本文建立了若干孔隙数字岩心、裂缝数字岩心和裂缝-孔隙数字岩心,并采用CPU-GPU异构并行模拟计算数字岩心的弹性模量.CPU-GPU异构并行计算可以迅速模拟大量数字岩心的弹性模量,进而分析弹性模量与孔隙度等特征参数的关系.这无论是对岩石物理学自身的发展,还是对油气勘探实践都具有十分重要意义.

3.1 孔隙数字岩心

孔隙数字岩心是指只含孔隙,不含裂缝的数字岩心.本文以砂岩样品的CT扫描图像为基础建立一个数字岩心,其骨架设为均匀介质.为使数字岩心能采样较大的岩心空间,同时又能减少像素数量,本文对数字岩心进行重采样.重采样后的数字岩心的像素总量减小至原来的1/8,其大小为1503像素,孔隙度为0.212.图 1中展示的即是所建三维数字岩心的图像.

为了增加数字岩心的样本数量,本文对数字岩心的骨架进行膨胀处理,简单模拟某种胶结物可能的生长方式,产生多个数字岩心.胶结物生长遵循如下述原则.当一个孔隙像素具有N个相邻骨架点时,它的像素值被设为N×I.其中,I大于原始数字岩心中最大的像素值.然后,将像素值大于给定阈值的孔隙像素转换为胶结物像素.这种方法可以十分方便地控制胶结物含量.以图 1中的数字岩心为基础,本文建立了3个含胶结物的数字岩心.图 2展示了其中1个数字岩心的局部图像.

图 1 三维数字岩心 Figure 1 Three dimensional digital rock
图 2 含胶结物数字岩心的局部图像 Figure 2 Local image of digital rock with cement

当采用LESS计算数字岩心的弹性模量时,孔隙流体的体积模量和剪切模量分别取0.00001 GPa和0.0 GPa;骨架矿物的体积模量和剪切模量分别取37.0 GPa和44.0 GPa;胶结物的体积模量和剪切模量定义3组,分别是第一组(29.6 GPa、35.2 GPa),第二组(37.0 GPa、44.0 GPa)和第三组(44.4 GPa、52.8 GPa).模拟结果表明这些数字岩心均存在微弱的各向异性,但本文近似按照各向同性情况处理.数字岩心的弹性模量与孔隙度的关系如图 3所示.

图 3 数字岩心的弹性模量与孔隙度关系 Figure 3 Relationship between elastic moduli and porosity of digital rocks

图 3中将同一胶结物对应的数字岩心的体积模量和剪切模量分别用虚线相连.

由于定义了3种弹性模量不相同的胶结物,故在同一孔隙度位置存在3组体积模量和剪切模量数据,即6个数据点.胶结物的弹性模量越大,数字岩心的弹性模量越大.当胶结物含量较高时,在同一孔隙度位置数字岩心的弹性模量数值差异较大;当胶结物含量较低时,弹性模量的数值比较接近.这体现了胶结物的含量与弹性模量对数字岩心弹性模量的影响.

图 3可知,几何结构相同的数字岩心,当胶结物的弹性模量变化时,数字岩心的弹性模量也将发生变化.在Xu-White理论中,将含泥质砂岩中的孔隙分为砂岩孔隙和泥岩孔隙,通过改变两类孔隙的扁率来拟合实际数据.如果采用Xu-White理论拟合模拟结果,那么当数字岩心的几何结构相同时,反演的孔隙扁率也会发生变化.这说明孔隙扁率不仅仅是一个描述孔隙形状的几何参数.

3.2 裂缝数字岩心

裂缝数字岩心是指在均匀介质中嵌入裂缝所形成的模型.在地层中裂缝系统的尺度变化非常大.大型裂缝的尺度远远超过岩心的尺寸,而微裂缝往往又在CT的成像分辨率之外.尺度差异较大的两组裂缝系统不能直接包含在同一数字岩心中.前人一般采用简单的几何形状来表示裂缝,如Grechka等(2006)将嵌入包裹体的均匀介质作为数字岩心计算弹性模量;Orlowsky等(2003)在均匀介质中设置线状裂缝建立二维数字岩心进行波场模拟.仿照这种思路,本文首先生成硬币状裂缝面,然后将它绕一个轴旋转,最后嵌入到三维均匀介质中,形成裂缝数字岩心.裂缝面的直径和空间方向的变化量都是在规定的范围内随机选取的.利用这种方法本文建立了4个裂缝数字岩心.图 4展示了其中1个裂缝数字岩心的三维图像.受GPU卡计算性能的限制,本文所建的裂缝数字岩心的分辨率较低,致使裂缝的张开度较大,裂缝孔隙度较高.

图 4 裂缝数字岩心 Figure 4 Cracked digital rock without pore

当采用LESS计算数字岩心的弹性模量时,流体的体积模量和剪切模量分别取0.00001 GPa和0.0 GPa;骨架矿物的体积模量和剪切模量分别取37.0 GPa和44.0 GPa.裂缝数字岩心的弹性模量与孔隙度之间的关系如图 5所示.在图 5所示的孔隙度范围之内,4个裂缝数字岩心的弹性模量差异不明显,这说明裂缝数字岩心的建模方法比较稳定.

图 5 裂缝数字岩心的弹性模量与孔隙度关系 Figure 5 Relationship between elastic moduli and porosity of cracked digital rocks
3.3 裂缝-孔隙数字岩心

裂缝-孔隙数字岩心是指同时包含裂缝和孔隙的数字岩心.本文将3.1节中图 1所示的数字岩心分别与3.2节中的4个裂缝数字岩心进行布尔运算,即孔隙ⓧ裂缝→孔隙,孔隙ⓧ骨架→孔隙,骨架ⓧ裂缝→裂缝,骨架ⓧ骨架→骨架,生成4个裂缝-孔隙数字岩心.图 6是其中的1个裂缝-孔隙数字岩心的三维图像. 当采用LESS计算数字岩心的弹性模量时,流体的体积模量和剪切模分别取0.00001 GPa和0.0 GPa,骨架的体积模量和剪切模分别取37.0 GPa和44.0 GPa.裂缝-孔隙数字岩心的弹性模量与孔隙度的关系如图 7所示.与3.2节中的裂缝数字岩心相比,裂缝-孔隙数字岩心的弹性模量相对变化较大.这说明裂缝系统和孔隙系统的融合增加了孔隙结构的变化.

图 6 裂缝-孔隙数字岩心 Figure 6 Cracked digital rock with pore
图 7 裂缝-孔隙数字岩心的弹性模量与孔隙度关系 Figure 7 Relationship between elastic moduli and porosity of cracked digital rocks with pore
4 讨论与结论

4.1 本文将LESS以CPU-GPU异构并行的方式实现,模拟若干数字岩心的弹性模量,并简单分析弹性模量与孔隙度关系.

4.2 本文认为相对低成本的CPU-GPU异构并行计算对于LESS以及数字岩石物理的应用具有重要的意义,具体如下:

(1) 多数岩石物理模型都是在各向同性假设前提下建立的,因此一般要求数字岩心为各向同性.实际上,本文建立的几类数字岩心均存在一定程度的各向异性.如果实际岩心确实具有各向异性,那就要采用适应各向异性的方法.如果各向异性是由于数字岩心包含的微观结构类型不充分引起的(即采样的岩心空间不够大),那就需增大数字岩心的尺寸.

(2) 岩石的微观结构异常复杂,只有当数字岩心具有足够高的分辨率和足够大的尺寸时,计算的宏观物理参数才能合理,否则将对经验关系、理论模型的评价与建立产生不利影响.扩大数字岩心的采样空间一般会增大数字岩心的尺寸.

(3) 数字岩心的尺寸一般小于实验室标准岩样,孤立的模拟结果的意义有限.若要形成解决实际问题的能力,必须计算足够多的样本,建立不同宏观物理量之间的相关性.

因此,LESS以及数字岩石物理对计算效率性价比的提高具有迫切性,CPU-GPU并行的LESS可以迅速计算大量数字岩心的宏观弹性模量,提供相当于物理实验的“观测数据”,这对岩石物理学研究具有重要的意义.

致谢 感谢匿名审稿专家的支持和帮助.
参考文献
[1] Andrä H, Combaret N, Dvorkin J, et al.2013. Digital rock physics benchmarks-Part Ⅱ: Computing effective properties[J]. Computers&Geosciences, 50 (SI) . DOI:10.1016/j.cageo.2012.09.008
[2] Arns C H, Knackstedt M A, Pinczewski W V, et al.2002. Computation of linear elastic properties from microtomographic images: Methodology and agreement between theory and experiment[J]. Geophysics, 67 (5) : 1396–1405. DOI:10.1190/1.1512785
[3] Bohn R B, Garboczi E J. 2003. User manual for finite element and finite difference programs: A parallel version of NISTIR 6269[R]. National Institute of Standards and Technology Internal Report 6997.
[4] Chen Z X, Meng X H, Guo L H, et al.2012. Three-dimensional fast forward modeling and the inversion strategy for large scale gravity and gravimetry data based on GPU[J]. Chinese Journal of Geophysics (in Chinese), 55 (12) : 4069–4077. DOI:10.6038/j.issn.0001-5733.2012.12.019
[5] Dvorkin J, Derzhi N, Diaz E, et al.2011. Relevance of computational rock physics[J]. Geophysics, 76 (5) . DOI:10.1190/GE02010-0352.1
[6] Garboczi E J, Day A R.1995. An algorithm for computing the effective linear elastic properties of heterogeneous materials: Three-dimensional results for composites with equal phase Poisson ratios[J]. Journal of the Mechanics and Physics of Solids, 43 (9) : 1349–1362. DOI:10.1016/0022-5096(95)00050-S
[7] Garboczi E J. 1998. Finite element and finite difference programs for computing the linear electric and elastic properties of digital images of random materials[R]. National Institute of Standards and Technology Internal Report 6269.
[8] Grechka V, Vasconcelos I, Kachanov M.2006. The influence of crack shape on the effective elasticity of fractured rocks[J]. Geophysics, 71 (5) . DOI:10.1190/1.2240112
[9] Jiang L M, Sun J M, Liu X F, et al.2012. Numerical study of the effect of natural gas saturation on the reservoir rocks' elastic parameters[J]. Well Logging Technology (in Chinese), 36 (3) : 239–243.
[10] Kirk D B, Hwu W W.2010. Programming Massively Parallel Processors: A Hands-on Approach[M]. Singapore: Elsevier Inc .
[11] Knackstedt M A, Arns C H, Pinczewski W V.2003. Velocity-porosity relationships, 1: Accurate velocity model for clean consolidated sandstones[J]. Geophysics, 68 (6) : 1822–1834. DOI:10.1190/1.1635035
[12] Liu S W, Wang H Z, Chen S C, et al.2013. Implementation strategy of 3D reverse time migration on GPU/CPU clusters[J]. Chinese Journal of Geophysics (in Chinese), 56 (10) : 3487–3496. DOI:10.6038/cjg20131023
[13] Liu X F. 2010. Numerical simulation of elastic and electrical properties of rock based on digital cores [PhD Thesis] (in Chinese). Dongying: China University of Petroleum (East China).
[14] Liu X J, Zhu H L, Liang L X.2014. Digital rock physics of sandstone based on micro-CT technology[J]. Chinese Journal of Geophysics (in Chinese), 57 (4) : 1133–1140. DOI:10.6038/cjg20140411
[15] Orlowsky B, Saenger E H, Guéguen Y, et al.2003. Effects of parallel crack distributions on effective elastic properties - A numerical study[J]. International Journal of Fracture, 124 (3-4) . DOI:10.1023/B:FRAC.dqwlxjz-31-4-178322563.29991.80
[16] Sun J M, Yan G L, Jiang L M, et al.2014. Research of influence laws of fluid properties on elastic parameters of fractured low permeability reservoir rocks based on digital core[J]. Journal of China University of Petroleum (Natural Science Edition) (in Chinese), 38 (3) : 39–44. DOI:10.3969/j.issn.1637-5005.2014.03.006
[17] Tang X D, Liu H, Zhang H.2015. Frequency-space domain acoustic modeling based on an average-derivative and GPU implementation[J]. Chinese Journal of Geophysics (in Chinese), 58 (4) : 1341–1354. DOI:10.6038/cjg20150421
[18] Zhang J Y, Sun J M.2012. Rock elastic properties determined by using digital rock and effective medium model[J]. Journal of Oil and Gas Technology (in Chinese), 34 (2) : 65–70.
[19] Zhang Y, Toksöz M N.2012. Impact of the cracks lost in the imaging process on computing linear elastic properties from 3D microtomographic images of Berea sandstone[J]. Geophysics, 77 (2) : R95–R104. DOI:10.1190/GE02011-0126.1
[20] Zhao J P, Sun J M, Jiang L M, et al.2014. Effects of cementation on elastic property permeability of reservoir rocks[J]. Earth Science-Journal of China University of Geosciences (in Chinese), 39 (6) : 769–774. DOI:10.3799/dqkx.2014.072
[21] Zhu W, Shan R.2014. Progress of digital rock physics[J]. Oil Geophysical Prospecting (in Chinese), 49 (6) : 1138–1146.
[22] 陈召曦, 孟小红, 郭良辉, 等.2012. 基于GPU并行的重力?重力梯度三维正演快速计算及反演策略[J]. 地球物理学报, 55 (12) : 4069–4077. DOI:10.6038/j.issn.0001-5733.2012.12.019
[23] 姜黎明, 孙建孟, 刘学锋, 等.2012. 天然气饱和度对岩石弹性参数影响的数值研究[J]. 测井技术, 36 (3) : 239–243.
[24] 刘守伟, 王华忠, 陈生昌, 等.2013. 三维逆时偏移GPU/CPU机群实现方案研究[J]. 地球物理学报, 56 (10) : 3487–3496. DOI:10.6038/cjg20131023
[25] 刘学锋. 2010. 基于数字岩心的岩石声电特性微观数值模拟研究[博士论文]. 东营: 中国石油大学(华东).
[26] 刘向君, 朱洪林, 梁利喜.2014. 基于微CT技术的砂岩数字岩石物理实验[J]. 地球物理学报, 57 (4) : 1133–1140. DOI:10.6038/cjg20140411
[27] 孙建孟, 闫国亮, 姜黎明, 等.2014. 基于数字岩心研究流体性质对裂缝性低渗透储层弹性模量的影响规律[J]. 中国石油大学学报(自然科学版), 38 (3) : 39–44. DOI:10.3969/j.issn.1637-5005.2014.03.006
[28] 唐祥德, 刘洪, 张衡.2015. 基于加权平均导数的频率-空间域正演模拟及GPU实现[J]. 地球物理学报, 58 (4) : 1341–1354. DOI:10.6038/cjg20150421
[29] 张晋言, 孙建孟.2012. 应用数字岩心和有效介质模型研究岩石弹性性质[J]. 石油天然气学报, 34 (2) : 65–70.
[30] 赵建鹏, 孙建孟, 姜黎明, 等.2014. 岩石颗粒胶结方式对储层岩石弹性及渗流性质的影响[J]. 地球科学(中国地质大学学报), 39 (6) : 769–774. DOI:10.3799/dqkx.2014.072
[31] 朱伟, 单蕊.2014. 虚拟岩石物理研究进展[J]. 石油地球物理勘探, 49 (6) : 1138–1146.