CN1302386C - 低浓度表面活性剂与相态结合的三元复合驱计算机仿真方法 - Google Patents

Abstract

一种低浓度表面活性剂与相态结合的三元复合驱计算机仿真方法，包括一个完备的三元复合驱计算机仿真模型，与模型配套的求解计算方法和实现程序，以及仿真前后数据处理技术。建模时，方法采用了关键字数据文件输入、变长度数组生成、变厚度参数与死节点处理以及仿真模型自动生成技术；求解模型时采用了简化Newton法、参数Newton法和“分裂算法”等数值计算方法，避免了传统方法中经常出现的代码错误及奇异矩阵问题，有效地降低了非线性方程的阶数，将仿真速度提高了2～5倍；仿真给出的“三元复合驱段塞注入程序及井工作制度文件”，可有效指导油田后续工业化开采时的方案设计、实际开采、驱油效果分析和指标预测。

Description

低浓度表面活性剂与相态结合的三元复合驱计算机仿真方法

技术领域

本发明涉及一种用于石油开采的计算机仿真方法，特别是一种低浓度表面活性剂与相态相结合的三元复合驱计算机仿真方法。

背景技术

在未来相当长的时间内，石油、天然气等仍会在能源工业中扮演相当重要的作用。世界各国都在千方百计地勘探、寻找此类能源的新产地，或通过挖潜增加产量、延长老油气田寿命，尤其对于石油储量不太丰富的国家更是如此。我国以大庆油田为例，它经过30年(以年产5000万吨原油的产量)连续不断地开采，主力矿区目前基本上已进入高含水后期，但由于：(1)大庆油田的油藏是不同时期、不同规模、不同类型的河道砂岩交错叠置组成的，无论是成因还是横向分布状况都非常复杂，因而即使在传统的主采区，仍有存在着有开采价值的剩余油区块；(2)油层的纵向分布非均匀性严重，在此方向上的油层多达上百个，由于技术上的原因，有些层面尚未开采，尤其是其纵深层更是如此，也就是说在这些层面里也蕴藏着一些油藏等待开发；(3)油田分布面积较大，非均质性特别严重，根据地质学分析，在核心区周围还存在着一些“不连续油藏构造”。这说明大庆油田地下现在还埋藏着相当数量的剩余油，但其开采难度却越来越大，弄清油田剩余油的精细分布(特别是那些薄差油层、含水油层和位于传统核心区周边地区的离散分布油藏)非常重要，而能把剩余油开采出来，继续确保油田的高产、稳产及可持续发展，延缓油田的衰老，把地下油藏变成现实的能源财富则非常必要。然而过去许多油田无论是在石油勘探、还是开发方面，多采用“经验-预测分析-小试-中试”为主的“较为初级的”科研模式，其中的许多数据处理任务也多以“手工作业”为主。采用这种方法存在一些缺陷：(1)很难给出能准确、全面描述“地质结构非常复杂、非均质特性严重”这类油田实际的驱油体系的数学模型；(2)即使能给出简化的数学模型，但建模过程中原始数据准备的工作量非常巨大，求解模型也非常困难，不仅枯燥、耗时，甚至根本就无法进行。所以，用这种研究方法很难避免盲目性，研究的效费比也不高，因而很难适应象大庆油田这样的超大型企业持续发展的需要。由于上述两个原因，必须采用全新的、能把“人的智慧”与“计算机的数据存储与处理能力”结合起来、优势能够互补的研究方法，例如“计算机仿真”方法，来大大提高研究的质量，减少石油开采的盲目性，避免造成人力、物力、财力上的浪费。大庆油田从“八五”开始就开展了这方面的研究工作，尤其是三元复合驱数值模拟，并已显示出了这种方法的良好应用前景。但随着三元复合驱技术的发展，原有的“数值仿真软件”已经不适合大规模的矿场应用；1997年，申请人前后分别从加拿大Gcomp公司和美国德州大学先后引进了“ASP(Alkali-Surfactant-Polymer)碱—表活剂—聚合物复合驱模型”和“UTCHEM(University Texas Chemistry)三元复合驱油软件”，但是这两种模型及软件都是为学校教学服务的，仍然存在一些问题：(1)对稀体系、非相态驱油体系模型的描述能力较差，无法满足油田实际仿真的需要；(2)对仿真模型(即非线性微分方程组)的求解方法不稳定、计算速度慢，还经常因遇到“奇异矩阵求逆”及“g1、g2、s1、s2代码出错”等问题，经常致使仿真无法正常进行；(3)缺乏对仿真输入/输出数据进行有效处理(即前后数据处理)的数据加工程序，以致影响“模型的建立”及对仿真结果数据的利用；(4)实用性及商品化程度比较低，应用范围非常有限。针对ASP及UTCHEM所存在的缺陷，必须进行科技创新，发明一种技术更加先进的三元复合驱计算机仿真方法，该方法应具有“所采用的仿真数学模型功能齐全、模型的求解方法先进实用、所得到的仿真结果对油田的工业化开采有重要指导意义以及便于应用和普及”等特点，以便为油田的可持续发展提供强有力的技术支持。

发明的内容

本发明的目的在于提供一种既能描述稀体系非相态、又能描述浓体系相态的功能齐全的计算机仿真数学模型(非线性微分方程组)及具有强大数据处理能力(前处理)的建模技术，提供一个能对所述模型进行快速、高精度、可靠求解计算的“计算方法”，以及一种能有效利用所述仿真结果的数据处理(后处理)技术，以适应表面活性剂注入浓度比较低的油田对三元复合驱技术研究和实际应用的需要，为油田的可持续发展提供强有力的技术支持。

为实现上述目的，本发明采用如下技术方案：

一种低浓度表面活性剂与相态结合的三元复合驱计算机仿真方法，包括主、从计算机系统、计算机仿真软件包DQCHEM、计算机操作系统，主、从计算机系统包括一个位于前台的PC计算机和位于后台的工作站，以及与仿真配套的计算机外部设备，前台计算机中存放着系统仿真服务程序SSSP，通过人机对话对仿真任务进行组织及过程控制，后台工作站完成仿真过程中的各种数据处理及复杂计算任务，其特征是：所述的DQCHEM包括复杂油藏地质结构描述软件CGCDSW、三元复合驱体系模型描述及生成软件ASPEXTMDG、化学平衡反应方程组求解计算软件CHEMEQCNT、仿真前后数据处理软件FDEDPSW；所述的ASPEXTMDG既能描述相态、也能描述非相态驱油体系的理化平衡反应；用DQCHEM进行三元复合驱计算机仿真，包括如下主要步骤：

1)启动主、从计算机，系统仿真服务软件对系统初始化，待显示器上显出人机对话菜单后，操作人员键入具体仿真命令，由主机生成仿真控制字SWMCW，并设置初始迭代时间步长；

2)主机根据仿真控制字从数据源获得所论仿真区块的实际油藏地质特征、仿真节点数、所用驱油材料在驱油体系中所占组分、各类钻井特征及所用驱油段塞注入程序等基础数据，接着由FDEDPSW之关键字数据文件输入软件、变长度数组生成软件、变厚度参数与死节点处理软件对基础数据进行预处理，形成仿真输入数据流，送至工作站的CHEMEQCNT程序入口处；

3)先完成压力方程系数矩阵元素的计算，然后根据驱油体系中有无表活剂及其浓度的大小，对仿真驱油体系进行相态、非相态判断；

4)当注入到驱油体系中的表活剂后使临界胶束浓度CMC≤某一数值时，驱油体系为非相态，则按下述“关联度建模法”建立该体系的计算机仿真数学模型；当CMC＞某一数值时，则直接从5)开始进行计算机仿真计算：

(1)根据每一口井位的地质特性及其不同的地下坐标，先把整个驱油体系用网格划分为界面特性各不相同的k个“子驱油体系块”(k＝1、2、3……)，并用k个子化学平衡反应模型CHEMEQSM(k)对大系统中的每个子系统进行描述；

(2)根据在实验室条件下实测得到的“表活剂、碱、聚合物与界面张力活性关系图”，计算出每个子区块的界面张力等参数数据，作为原始输入数据，经SSSP处理形成数据流后，送CHEMEQCNT的程序入口，并对每个子化学平衡反应模型CHEMEQSM(k)的参数进行赋值；

(3)以k为“循环控制”和“关联控制”变量，按下述步骤依次对每个化学平衡反应模型CHEMEQSM(k)进行求解计算及解的关联重组，得出总体采油方案与井工制度；

5)在SSSP的监控下，先用“仿真输入数据流”对三元复合驱仿真模型ASPEXTMDG中的相对渗透率、粘度、浓度、饱和度、毛管力等变量赋值，并计算传导率、对流中的浓度项，然后用CHEMEQCNT依次进行以下求解计算：

(1)根据压力方程PRSSUREQ求解相压力、相流速和井的压力限制、井的注入和产出速率；

(2)根据物质平衡方程MASBLCEQ、浓度方程，求解出各种组分的总浓度、注入的PV数及采出程度、相比例和井筒产出液的组成；

(3)根据能量守恒方程ENGCSVEQ进行化学平衡计算，解出新的离子平衡环境；

(4)利用化学平衡模型计算有效含盐量、醇的分配系数和其它相关组分的浓度；

(5)进行“相态闪蒸计算”，求解相饱和度和相组成；

(6)利用上述结果计算界面张力、毛管数、相剩余饱和度、相对渗透率、相密度、相粘度、渗透率下降系数；

6)将本次的仿真结果存储于专用缓存区BUFFER的规定存储单元CELL(i)中；

7)判断当前仿真井的工作状态，即根据其是否完井，或是否已关井/开井，如已关井，则依次读入下一口井的类型及参数，即注入井、生产井、定压井、定产井、射孔方向(x，y，z)、井的工作状态，和井的运行参数和段塞成分，并计算下一次迭代所需压力方程系数矩阵的元素与迭代时间步长，经动态更新后，执行从4)开始的新一轮迭代计算；

8)以仿真推演时间变量T是否达到设定的最大值Tmax为判据，判断要否继续进行仿真，如需继续迭代，则更新时间变量后执行从2)开始的循环，直到迭代计算结果满足仿真要求；

9)用FDEDPSW对暂存于BUFFER之CELL(i)中的数据进行结果处理，生成“三元复合驱采油综合工艺文件”，然后以多种媒体的形式输出对比结果。

用所述的三元复合驱仿真模型ASPEXTMDG的求解计算软件CHEMEQCNT求解非线性微分方程时，可视情况分别采用如下方法编制相应的求解软件：

1)在用Newton法求解非线性方程组时，如此类方程组的F′(x)计算非常复杂但总是满秩矩阵，可用简化Newton法求解此类非线性方程组，即使用同一个矩阵F′(x)连续迭代m次然后再更新其矩阵元素；

2)当求解相态驱油体系模型时，可利用显式解法(IMPEC)求解浓度方程，用闪蒸法计算相饱和度和相组成，但在求解模型时若碰到F′(x*)奇异的情形时，就在参数ω的收敛区间0＜ω＜2内，适当选择参数ω的值，并且根据求解精度及迭代时间的要求，采用“变步长技术”选择适当的迭代时间步长，然后再用最小二乘法求解，避开“矩阵F′(x*)的奇异”问题，以使方程求解顺利而快速的进行；

3)当需要求解非相态驱油体系模型时，可利用显式法(IMPEC)求解浓度方程、隐式法求解压力方程，或用等效拆分全隐式方法求解非相态驱油体系模型，其要领和具体过程是：假定对仿真区块划分网格时符合笛卡尔坐标离散条件，将总浓度的求解公式转化为近似求解“分相组分”浓度的相关表达式，求出n+r(0＜r＜1)时刻对流方程的解后，再利用n+r时刻的对流方程解，求出n+1时刻扩散方程解。

所述的FDEDPSW包括输入数据前处理软件FDPSW，计算结果后处理及显示软件EDPSW，其中FDPSW包括关键字数据文件输入软件、变长度数组生成软件、变厚度参数与死节点处理软件，EDPSW包括一维曲线图、二维色彩图、三维立体图生成软件；

1)用所述的FDPSW对输入数据进行预处理的方法是：

(1)先将输入数据文件按类型分成“题头注释和油藏数据、输出选件信息、油藏性质、一般物理性质数据、地质化学反应剂的物化参数、生物降解数据、井数据”七个向量数组，每个向量数组设一关键字，然后对关键字进行排序；

(2)根据对关键字的排序结果，依次把与关键字对应的、字长各不相同的向量数组里的数据，自动写入到“容量动态可变的数据暂存区”中，并标出各向量数组在所述存储区中的其起始地址；

(3)将带有“关键字”和“地址指针”标志且暂存于“数据暂存区”中原始基础数据，按地址指针进行连接，生成一个带有所述标志的一维大数组，然后以数据流的形式把它送至后台工作站的CHEMEQCNT程序入口处；

(4)将带有“关键字”和“地址指针”标志的一维大数组，按标志分解成与基础数据对应的向量数组，对相关数学模型的参数赋值后，用CHEMEQCNT正式进行仿真计算；

2)用所述的EDPSW对计算结果进行处理的方法是：

(1)用VC⁺⁺中的类CcurveView、CtwoView、CcubeView，将暂存于BUFFER中的数据分别绘制成一维曲线图、二维色彩图、三维立体图，以便给出全井及单井的注入量、产液量、产水量、产油量等与时间的对应关系，或把二维截面上的数据显示为平面彩色图和等值线图，把三维数据转换为可随意旋转的三维立体图形，以便更直观的显示仿真结果；

(2)将暂存于BUFFER中的“仿真结果数据”和“三元复合驱采油实验区所用采油工艺”数据进行拟合对比，生成“三元复合驱段塞注入程序及井工作制度文件”；

(3)保存经EDPSW处理后的出的结果，然后或在终端设备上打印、显示处理结果，或通过网络化仿真服务软件，将处理结果送局域网上，供注册用户观看、查阅。

附图说明

图1是本发明的总体组成示意方框图；

图2是本发明的仿真过程流程图；

图3是杏二西区块仿真结果与实际驱采结果间的拟合对比曲线；

具体实施例

鉴于本发明涉及的技术层面非常广，既要交代与实现等计算机仿真技术直接相关的“仿真数学模型及其建立、计算方法、程序设计与实现”等，还要介绍不同“驱油体系”的“驱油机理”，可以说是多个学科的交叉。注意到大庆油田的油藏地质结构非常复杂的实际，必须用一系列非线性微分方程(化学平衡反应方程)描述所论“驱油体系”，才能切合油田的实际。这样一来，寻求速度快、精度高、运行可靠的“非线性微分方程(或方程组)”的求解方法，就成为本发明必须解决的技术关键。在介绍实施例时只讲具体的方法、步骤，就不容易深刻理解本发明，而如果在讲方法、步骤的同时混杂着也讲仿真模型、驱油机理、计算方法甚至具体程序，那就使人很难抓住要领。因而集中介绍一下本发明所用的仿真数学模型、驱油机理、计算方法、程序设计与实现，然后再结合实施例，集中介绍三元复合驱计算机仿真的具体方法。

1、本发明的组成及各部分的作用

本发明的技术创新点是适合油田实际的“三元复合驱计算机仿真模型建立”及“非线性微分方程的求解方法和计算机程序的研制”。为了解决这些关键技术，必须靠硬件和软件相互协调才能实现，图1给出了本发明的总体组成方框图。从图中可以看出：硬件部分包括主、从计算机系统，配套的计算机外部设备(多功能显示器、人机对话工具组)等，软件部分指的是“三元复合驱计算机仿真软件包DQCHEM”，它包括复杂油藏地质结构描述软件CGCDSW、三元复合驱仿真模型自动生成软件ASPEXTMDG、化学平衡反应方程组即“三元复合驱仿真模型”求解计算软件CHEMEQCNT，以及仿真服务与仿真前后数据处理软件SSSP & FDEDPSW。所述的主、从计算机系统包括一个位于前台的PC微型计算机和位于后台的计算机工作站WS，前台计算机中存放着系统仿真服务程序SSSP，通过人机对话对仿真任务进行组织及过程控制，后台工作站存放着仿真模型生成、模型求解以及仿真数据处理软件等，完成仿真过程中的各种数据处理及模型求解等复杂计算任务；主机通过调用SSSP及其网络化仿真服务子程，产生“仿真工作模式控制字SWMCW”，然后在SSSP组织及控制协调下，最终由后台工作站完成仿真过程中的化学平衡反应方程的求解计算与各种数据处理任务。

2、关于系统仿真服务程序SSSP

系统仿真服务软件SSSP在“仿真任务的组织及过程控制”方面具有重要作用，SSSP的作用是：(1)通过对人机对话工具状态的识别生成任务控制字；(2)根据任务控制字选择输入数据的来向和输出数据的流向；(3)对整个仿真过程进行控制协调。为了完成这些任务，它必须具有多重中断响应控制子程序、人机对话界面生成子程序、数据处理作业控制子程序、仿真过程监控子程序、网络化仿真服务子程序等，才能从技术上保证仿真工作的正常进行。虽然多重中断响应控制子程序、人机对话界面生成子程序、仿真过程监控子程序很重要，但人们比较熟悉，故不赘述，下面重点介绍其中的网络化仿真服务子程序和数据处理作业子程序。

2.1网络化仿真服务子程

该程序是用以控制输入/输出数据的流向、进而确定是单站独立仿真，还是组网仿真的，其核心部件是两个独立的、受前台计算机控制的“单刀双掷数据选择软件开关”。根据仿真数据输入/输出方式的不同，可给出四种仿真工作模式控制字(SWMCW＝00、01、11、10)，然后“数据选择软件开关”根据工作模式控制字的状态，通过切换决定输入/输出数据的来/去向，从而得出四种仿真服务模式：

“00”输入数据由本地提供——输出数据提供给本地；

“01”输入数据由本地提供——输出数据提供给局域网；

“10”输入数据由局域网提供——输出数据提供给本地；

“11”输入数据由局域网提供——输出数据提供给局域网；

2.2数据处理作业子程序

该程序是为仿真准备初始数据，对仿真所得结果进行数据整理，以减少仿真过程中的出错率，提高仿真效率、进而拓展仿真结果的利用率而专门开发的一个仿真前后数据处理软件FDEDPSW(也可称为数据加工软件)。FDEDPSW包括输入数据前处理软件FDPSW，计算结果后处理及显示软件EDPSW，其中FDPSW包括关键字数据文件输入软件、变长度数组生成软件、变厚度参数与死节点处理软件，EDPSW包括一维曲线图、二维色彩图、三维立体图生成软件。

在本发明以前的计算机仿真时，大多数方法(包括UTCHEM)输入给“化学反应方程的求解软件”的信息也是数据流，但由于这些数据流的生成和输入方式都比较落后(大多采用“固定格式的顺序输入法”；由于这种方法对“数据行”和“行中数据之间的顺序”要求非常严格，而且用户准备这些参数的工作量非常巨大，稍一疏忽就容易出错)，因此在DQCHEM模型中，我们采取了关键字的输入格式。由于化学驱油仿真不仅要用到数个专业的知识(非其它类型的油藏仿真所能比)，而且要和各种类型的输入数据(参数数目众多)打交道，因此关键字输入的难点在于设计。在我们研发的DQCHEM模型中，设计“关键字输入文件”的方法是：将原始数据先分成七组(题头注释和油藏数据、输出选件信息、油藏性质、一般物理性质数据、地质化学反应剂的物化参数、生物降解数据、井数据)，每组数据都设一“关键字”(经统计，这里所用的关键字多达400多个)，这就大大方便了用户的使用。但是，仅用关键字还有缺点，即它的所有数组长度都是固定的，一旦我们计算具有不同网格结构和不同组分的驱油体系的化学反应时，固定长度的数组就放不下这些数据，所以就必须重新定义数据的字长、修改数据处理源程序、以及重新对程序进行编译链接才能适应新情况，使用起来很不方便。为此，我们决定研发新的仿真模型DQCHEM，设计新的数据格式。开始采用的是公用块结构(模型使用了100多个公用块)，这种方法的缺点是：1)主程序可能会遗漏某些公共块；2)出现数据块重名现象，以致无法对模型进行顺利求解；3)某些公用块太短，只有一两个元素，使用很不规范。我们在建立DQCHEM模型时，一是借用“一维大数组”的设计思想，并用数组指针的办法为每个小数组分配存储位置，从而精简了一些多余的公用块，有效地避免了公用快的遗漏和重名；二是采用Fortran新版本所提供的直接分配数组空间的功能，为关键字数据文件动态设置“子数组”的边界值(需要多大的数组就开多大的数据长度)，最佳地使用计算机存储资源。实践证明：这种办法直观，可读性好，不易出错。

在仿真时还经常会碰到对厚度不同的油藏地质结构区块(非均匀性非常剧烈的不连续油藏地质区块)的仿真问题，这自然要求数据输入软件有对这类结构的描述能力。本数值仿真模型DQCHEM允许用两种格式输入油层垂向厚度：①如果层间厚度相等，则只需输入一个厚度常数；②层间厚度不等，需要输入N个厚度常数，N≤NZ。对第二种情况，尽管层间厚度可以不等，但对每一层来说，如各网格的厚度都一样，这显然不能准确地描述实际的油藏构造。我们增加第三种厚度数据输入格式——变厚度输入格式。对于一个垂向被划分为NX*NY*NZ个网格节点的油藏，X方向及Y方向尺寸按UTCHEM原始方法输入数据，Z方向输入NX*NY*NZ个厚度数据，即每一个网格节点就可以有对应的厚度值；为处理“死节点”问题，我们在输入数据格式中还特地增加了相应的关键字，定义一个与网格节点数等长的整型数组NKN(NX*NY*NZ)。在模型的初始化阶段，对该数组的所有元素赋初值为1。当要仿真的实际油藏需要处理死节点时，就从数据文件中读取死节点标记：即读入一个长度为NX*NY*NZ的整型数组，赋给数组NKN。数组元素为1表示正常节点，为0表示死节点。数据文件中，代表死节点的数据由油藏描述人员给出，或者由油田的数据文件编制人员根据油层的实际情况给出。以上就是处理非均匀油藏地质结构的输入数据时的编程思想，即“变厚度参数与死节点处理软件”的编制方法。

根据上述设计思路，用FDPSW对输入数据进行预处理的方法是：

1)先判断输入油藏地质结构特征数据时有无变厚度参数和死节点问题，如有，就先按上述方法处理“厚度参数和死节点”问题，接着再按下述方法处理；

2)将输入数据文件按类型分成“题头注释和油藏数据、输出选件信息、油藏性质、一般物理性质数据、地质化学反应剂的物化参数、生物降解数据、井数据”七个向量数组，每个向量数组设一关键字，然后对关键字进行排序；

3)根据对关键字的排序结果，依次把与关键字对应的、字长各不相同的向量数组里的数据，自动写入到“容量动态可变的数据暂存区”中，并标出各向量数组在所述存储区中的其起始地址；

4)将带有“关键字”和“地址指针”标志且暂存于“数据暂存区”中原始基础数据，按地址指针进行连接，生成一个带有所述标志的一维大数组，然后以数据流的形式把它送至后台工作站的CHEMEQCNT程序入口处；

5)将带有“关键字”和“地址指针”标志的一维大数组，按标志分解成与基础数据对应的向量数组，对相关数学模型的参数赋值后，用CHEMEQCNT正式进行仿真计算；

上述方法就是所述的“关键字输入法”、“变长度数组生成法”、“变厚度参数与死节点处理方法”，实际上也就是编程实现的具体技术。

用EDPSW对CHEMEQCNT的计算结果进行处理的方法是：

1)用VC⁺⁺中的类CcurveView、CtwoView、CcubeView，将暂存于BUFFER中的数据分别绘制成一维曲线图、二维色彩图、三维立体图；

2)将暂存于BUFFER中的“仿真结果数据”和“三元复合驱采油实验区所用采油工艺”数据进行拟合对比，生成“三元复合驱段塞注入程序及井工作制度文件”；

3)保存经EDPSW处理后得出的结果(如果需要，也可把一维曲线图转换成EXCEL电子表格的形式，直观地给出全井及单井的注入量、产液量、产水量、产油量等与时间的对应关系及驱油段塞注入程序；也可从二维色彩图的数据中提取出等值线图数据，把三维数据转换为可随意旋转的三维立体图形)，然后或在终端设备上打印、直观显示处理结果，或通过网络化仿真服务软件，将处理结果送局域网上，供注册用户观看、查阅，以便共享资源。

2.3仿真过程监控子程序

在SSSP中设有仿真过程监控子程序，一方面用来协助仿真计算的正常进行，再者，当仿真过程中出现问题时，向操作人员发出提示信息。虽然本程序很重要，但由于它的功能比较简单，也不是要重点保护的内容，只是为了内容的完整在此稍提一下，以便全面地了解整个仿真过程。

3、三元复合驱计算机仿真模型介绍与驱油机理描述

要进行计算机仿真，必须有仿真的数学模型。本发明给出的“三元复合驱计算机仿真数学模型DQCHEM”，对相态类型、驱油机理(表面活性剂驱、碱驱、聚合物驱等)的描述都比较全面；在对相态类型描述方面，模型既考虑到浓体系相态、也考虑到稀体系非相态，即它除考虑水相、油相、微乳液相外，对对碱、聚合物和表面活性剂的功能及它们之间的关系，也都给以了充分描述；此外它对示踪剂、凝胶等功能的描述也较为成功，因而本模型既适用于由“高浓度表面活性剂”所生成的相态“驱油体系”、也适用于由“低浓度表面活性剂”所生成的非相态“驱油体系”之情况，真正做到了相态与非相态的结合，因而适应范围更广；该模型对驱油机理的描述也非常全面，无论“驱油体系”处在什么状态(相态、非相态)，都可把驱油机理描述成：在驱油过程中降低油水两相之间的界面张力，进而降低毛细管力对油的捕集作用，使被驱替过的油层剩余油饱和度降低，渗透能力提高。

3.1三元复合驱计算机仿真所用的基本仿真数学模型

本发明所研究的对象，是一个涉及到相互间可发生复杂物化反应的、由“地下水、地下原油、注入水、注入表面活性剂、岩石、碱等多种物质”所构成的化学反应体系，这一化学平衡体系不但是一个严重的非线性系统，同时组分浓度之间的变化差异十分巨大。本发明所用的基本仿真数学模型DQCHEM，就是以“流体力学、物理化学、聚合物流变学”之“三大方程(即物质平衡方程MASBLCEQ、能量守恒方程ENGCSVEQ、压力方程PRSSUREQ)”为依据，由“非线性微分方程组”确定的油田实际油藏地质结构、所用驱油材料以及由它们所形成的“液态驱油体系”和体系中原油的渗出特性(油水两相之间的界面张力、毛细管力对油的捕集力、油层剩余油饱和度等)等参量之间的相互制约关系。

(1)物质平衡方程MASBLCEQ

根据达西定律，K组分质量连续性可根据单位孔隙体积

中k组分的总体积表达式为：

$\frac{\∂}{\∂t}(\mathrm{\φ}{\tilde{C}}_k\·{\mathrm{\ρ}}_k)+\▿\left[\underset{l=1}{\overset{n_p}{\mathrm{\Σ}}}{\mathrm{\ρ}}_k(C_{\mathrm{kl}}\·\stackrel{\‾}{U}-{\tilde{D}}_{\mathrm{kl}})\right]=R_k...(A.1-1)$

式中，单位孔隙体积k组分的总孔隙体积为包括吸附相的所有相之和。

${\tilde{C}}_k=[1-\underset{k=1}{\overset{n_{\mathrm{cv}}}{\mathrm{\Σ}}}\·{\hat{C}}_k]\underset{k=1}{\overset{n_{\mathrm{cv}}}{\mathrm{\Σ}}}{S}_l\·{\hat{C}}_{k}l+C_k,k=1,2,\·\·\·,\mathrm{nc}...(A.1-2)$

其中：n_cv，占体积组分总数，这些组分为水、油、表活剂、空气；_np为相数；k组分吸附浓度；ρ_k，纯k组分在参考压力(1atm)下的密度，S_l为相态饱和度。

这里，输入参数为n_cv、n_p、ρ_k、φ、C_kl、S_l为输入变量， U、D_kl为渗流力学中的达西定律及数学离散表达式，

为中间变量(见公式B.1.1-1)。在进行数值仿真计算时，首先在数据流中给出实际区块的上述输入参数，再根据达西定律及离散数学表达式，通过中间变量的结果，计算出达到一定迭代精度的组分浓度的数值结果。

(2)能量守恒方程ENGCSVEQ

假设能量仅为温度的函数且油相或液相中能量仅以对流和热传导方式进行，则能量守恒方程可推导为：

$\frac{\∂}{\∂t}[(1-\mathrm{\φ}){\mathrm{\ρ}}_s\·C_{\mathrm{vs}}+\mathrm{\φ}\underset{l=1}{\overset{n_p}{\mathrm{\Σ}}}{\mathrm{\ρ}}_l\·S_l\·C_{\mathrm{vl}}]T+\▿[\underset{l=1}{\overset{n_p}{\mathrm{\Σ}}}{\mathrm{\ρ}}_l\·C_{\mathrm{pl}}\·U_{l}T-{\mathrm{\λ}}_T\▿T]=\mathrm{qH}-Q...(A.1-3)$

式中：T，油藏温度；C_vs、C_vl为固体和相1在恒定体积下的热容量；λ_T为热传导率；qH为单位体积焓源项；Q₁为向上下盖层及固体中的热损失；

(3)压力方程PRSSUREQ

根据参考相(即水相)压力，压力方程可表述为

$\mathrm{\φ}{C}_t\frac{\∂P_l}{\∂t}+\▿\stackrel{\‾}{\stackrel{\‾}{K}}\·{\mathrm{\λ}}_{\mathrm{rTc}}\·\▿P_l=-\▿\underset{l=1}{\overset{n_p}{\mathrm{\Σ}}}\stackrel{\‾}{\stackrel{\‾}{K}}\·{\mathrm{\λ}}_{\mathrm{rlc}}\·\▿h+\▿\·\underset{l=1}{\overset{n_p}{\mathrm{\Σ}}}\stackrel{\‾}{\stackrel{\‾}{K}}\·{\mathrm{\λ}}_{\mathrm{rlc}}\·\▿P_{\mathrm{cll}}+\underset{k=1}{\overset{n_{\mathrm{cv}}}{\mathrm{\Σ}}}{Q}_k...(A.1-4)$

式中： ${\mathrm{\λ}}_{\mathrm{rlc}}=\frac{k_{\mathrm{rl}}}{{\mathrm{\μ}}_l}\·\underset{k=1}{\overset{n_{\mathrm{cv}}}{\mathrm{\Σ}}}{\mathrm{\ρ}}_k\·C_{\mathrm{kl}},{\mathrm{\λ}}_{\mathrm{rTc}}=\underset{l=1}{\overset{n_p}{\mathrm{\Σ}}}{\mathrm{\λ}}_{\mathrm{rlc}};$ C_t，总压缩系数；

弥散张量。

借助这三大方程，原则上就可对“驱油体系”中的三个液相(水、油、微乳液，其中微乳液相一般取决于相环境中各组分相对量和有效电解质浓度或含盐度)进行计算机仿真了。

3.2对驱油的机理的数学描述

由“表活剂”生成的“液态驱油体系”中有三个液相(水相、油相、微乳液相)，其他情况(如碱驱)可视为“表活剂驱”的的特例，所以这里我们主要用数学描述(即解析式)对“表面活性剂的驱油机理”及“驱油过程中相关参量的计算方法”进行说明。

3.2.1表面活性剂驱的物化机理

表面活性剂的主要驱油(物化)机理是在驱油过程中降低油水两相之间的界面张力，降低毛细管力对油的捕集作用，使被表面活性剂溶液驱替过的油层剩余油饱和度降低，表面活性剂驱油作用的过程为：表面活性剂、水和油所组成的体系首先根据体系环境的有效含盐量变化情况形成不同类型的相态，不同类型的相态组成决定了体系中相之间的界面张力，然后由界面张力确定了反应粘滞力作用效果的毛管数，毛管数的高低影响各相剩余饱和度的大小，最终使相之间的相对渗透能力发生了变化，进而导致所有相饱和度同时降低，流动能力相对提高。当驱油体系中表面活性剂浓度比较低时，体系就不会出现相的变化，整个体系也仅由油水两相组成，我们称这种体系为非相态稀体系；非相态稀体系三元复合驱的驱油机理是通过表面活性剂、碱、油、和水之间的协同效应来实现的。

1)吸附的表征方法

通常使用langmuir型吸附等温式来描述表面活性剂的吸附现象。表活剂吸附浓度可由下式给出：

${\hat{C}}_k=\min [{\tilde{C}}_k\·\frac{{\mathrm{\α}}_k({\stackrel{\‾}{C}}_k-{\hat{C}}_k)}{1+b_k({\tilde{C}}_k-{\hat{C}}_k)}],k=\mathrm{3,4}...(B1-1)$

吸附计算中的浓度是被水归一化的浓度，上述取最小值是为了保证吸附浓度随有效盐度呈线性增加，随渗透率增加而降低。

2)相态的表征

表活剂/油/水的相态理论主要考虑五种占体积组分(油、水、表活剂和两种醇)在溶液中形成三种拟组分的问题。如果没有醇，则仅仿真三种组分，通常被做成三角图。含盐度和二价离子浓度强烈地影响相态，在低盐度时，富油相基本上是纯油，微乳液相包括水和电解质、表活剂及一些溶剂油，这种相环境类型被称为Winsor I型或II(-)型。另外，如果表活剂浓度低于CMC，则两相为液相，包含所有表活剂、电解质及溶剂油和纯富油相，上述情况为低盐度时；对于高盐度，存在一富水相和包含了大部分表活剂、原油及一些溶解水的微乳液，这种相环境类型被称为WinsorII型或II(+)型；而II(-)和II(+)之间的相则被叫做第三相(WinsorIII)，这些相包含富集油相、富集水相和微乳液相。

当绘制双结点曲线和褶点线后，表活剂/油/水相态可表达为含盐度的函数，图2-1也给出了随着含盐量的增加相态变历的类型。当有两种相态或三种相态存在的情况下，计算方法分别如下：

(1)两相相态

对于II(-)型，可由双结点曲线方程给出：

$\frac{C_{3l}}{C_{2l}}=A\·{\left[\frac{C_{3l}}{C_{1l}}\right]}^B,1=\mathrm{1,2,3}...(B1-2)$

式中A、B为经验参数，C_kl为l相中k组分的浓度。

对于II(-)和II(+)的平衡相，可由褶点线方程给出：

$\frac{C_{3l}}{C_{2l}}=E\·{\left[\frac{C_{3l}}{C_{1l}}\right]}^F,1=\mathrm{1,2}...(B1-3)$

式中l＝1，2分别为II(-)和II(+)型相态，在缺少褶点线数据时，F＝-1/B；对于对称双结点曲线B＝-1时，F＝1；因此

$E=\frac{C_{1p}}{C_{2p}}=\frac{1-C_{2p}-C_{3p}}{C_{2p}}...(B1-4)$

式中C_kp为褶点处组分k的浓度，对II(-)和II(+)型的相环境均为输入值。

(2)三相相态

III型三相区相组成计算中，可简单假设富集油相和液相为单一线组成，微乳液相组成可由恒定不变点M的坐标确定，计算如下：

$\left\{\begin{array}{c}\frac{C_{3m}}{C_{2m}}=A{\left(\frac{C_{3m}}{C_{1m}}\right)}^B\\ C_{1m}+C_{2m}+C_{3m}=1\\ C_{2m}=\frac{C_{\mathrm{SE}}-C_{\mathrm{SEL}}}{C_{\mathrm{SEU}}-C_{\mathrm{SEL}}}\end{array}\right....(B1-5)$

式中C_SE为体系的含盐量，C_SEL是三相区开始形成的有效含盐量的下限，C_SEU为三相区消失时对应的有效含盐量上限。

3)相饱和度

存在表活剂的油层相饱和度可通过相浓度、总组分浓度和约束条件确定，因此，此时的相环境和相组成均已知道。总组分浓度和饱和度约束条件为：

$\left\{\begin{array}{c}{C}_k=\underset{l=1}{\overset{3}{\mathrm{\Σ}}}{S}_l\·C_{\mathrm{kl}}\\ \underset{l=1}{\overset{3}{\mathrm{\Σ}}}{S}_l=1\end{array}\right.,l=\mathrm{1,2,3}...(B1-6)$

渗流层中相饱和度(存在三相)可通过总组分浓度和饱和度约束条件确定：

$S_2=\frac{C_2-C_1}{1-C_{2l}}\·S_l=\frac{C_l}{1-C_{1l}},S_3=1-S_1-S_2...(B1-7)$

4)界面张力

在DQCHEM模型中，计算微乳液/水和微乳液/油之间界面张力的方法有Healy和Huh两种，这里我们给出Healy方法。

当相组成确定后，微乳液和富集相间的界面张力作为溶解参数的函数进行计算：

$\left\{\begin{array}{c}{\log }_{10\mathrm{\&sigma;l}3}={\log }_{l0}{F}_l+G_{l2}+\frac{G_{l1}}{1+G_{l3}\&CenterDot;R_{13}},R_{l3}>1\\ {\log }_{l0\mathrm{\&sigma;l}3}={\log }_{l0}{F}_l+(1-R_{l3})\&CenterDot;{\log }_{l0}{\mathrm{\&sigma;}}_{\mathrm{ow}}+R_{l3}\&CenterDot;[G_{l2}+\frac{G_{l1}}{1+G_{l3}}],R_{l3}<1\end{array}\right....(B1-8)$

式中：G_l1、G_l2、G_l3(l＝1，2)均为输入参数，R_l3为溶解比； $R_{l3}=\frac{G_{l3}}{G_{33}},$ 修正系数； $F_l=\frac{1-e-\sqrt{\mathrm{Co}{n}_l}}{e-\sqrt{2}}$ (l＝1，2)； ${\mathrm{Con}}_l=\underset{k=1}{\overset{3}{\mathrm{\&Sigma;}}}{(C_{k1}-C_{k3})}^2;$ 在没有表活剂或表活剂浓度小于CMC时，IFT＝σ_ow。

5)毛管压力

毛管压力与饱和度之间的关系与界面张力、渗透率和孔隙度密切相关，吸入过程与排出过程所反应的毛管压力特征是不一样的。通常把水驱和注入表活剂过程假定为吸入过程，因此考虑下列两种形式的毛管压力。

(1)渗流水层

存在三相时渗流水层毛管压力公式中隐含的假设是润湿性按水、有机物和空气方向逐渐降低的，且水相总是存在，即有：

${\left[\frac{P_b}{P_{c1l}}\right]}^{{\mathrm{\&lambda;}}_i}=1-S_{\mathrm{nl}},1=\mathrm{2,4}...(B1-9)$

式中λ_i为孔隙介质中孔隙大小分布测定值，压力P_b、P_c1l由岩石渗透率和孔隙度标定。

(2)饱和层

对于三相微乳液体系，饱和层的毛管压力可根据表活剂相态计算：

$\left\{\begin{array}{c}{\left[\frac{P_{b1}}{P_{c13}}\right]}^{{\mathrm{\&lambda;}}_i}=1-S_{\mathrm{nl}}\\ [{\frac{P_{b2}}{P_{c23}}]}^{{\mathrm{\&lambda;}}_i}=1-\frac{S_2-S_{2r}}{(S_1-S_{1r})+(S_3-S_{3r})}\\ S_{\mathrm{nl}}=\frac{S_1-S_{1r}}{1-S_{1r}-S_{2r}-S_{3r}}\\ P_{b1}=C_{\mathrm{pci}}\&CenterDot;\frac{\mathrm{\&sigma;l}3}{\mathrm{\&sigma;l}2}\&CenterDot;\sqrt{\frac{\mathrm{\&phi;}}{K}}\end{array}\right....(B1-10)$

6)捕集数

注表活剂提高采收率的一个重要机理就是捕集到的原油通过注表活剂使界面张力降低而产生移动，除了界面张力外，浮力也是影响油滴流动的重要因素，可以用Bond数表示。捕集并运动的非润湿相的Bond数和毛管数通常处理两组独立的无因次数组，一组为重力/毛管力(Bond数)，另一组为粘滞力/毛管力(毛管数)，两组结合起来发展为一种新的无因次数—捕集数。

计算公式如下：

毛管数： $N_{\mathrm{cl}}=\frac{\stackrel{\&OverBar;}{K}\&CenterDot;\&dtri;{\mathrm{\&Phi;}}_{l^{\&prime;}}}{\mathrm{\&sigma;}1l^{\&prime;}l}l=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n_p...(B1-11)$

l和l’为被驱替及驱替液体。

Bond数： $N_{\mathrm{Bl}}=\frac{\mathrm{kg}(\mathrm{\&rho;l}-{\mathrm{\&rho;l}}^{\&prime;})}{\mathrm{\&sigma;}{1l}^{\&prime;}}l=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n_p...(B1-12)$

7)相对渗透率

毛管数的变化会引起各相剩余饱和度的改变，毛管数与剩余饱和度之间的关系为：

$S_{\mathrm{lr}}={S_{\mathrm{lr}}}^H+\frac{{S_{\mathrm{lr}}}^L-{S_{\mathrm{lr}}}^H}{1+T_l{N}_{\mathrm{cl}}}l=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n_p...(B1-13)$

式中T_l是模型输入参数，S_lr ^L、S_lr ^H分别是低毛管数和理想极限高毛管数下l相的剩余饱和度，为模型输入参数。

每相剩余饱和度的变化必然引起相对渗透率曲线发生变化，变化后的l相相对渗透率曲线端点值k_r1和指数值n₁根据剩余饱和度值由低毛管数相对渗透率曲线和高毛管数相对渗透率曲线的端点值和指数值进行线形插值计算：

$\left\{\begin{array}{c}{K}_{\mathrm{rl}}={K_{\mathrm{rl}}}^L+\frac{{S_{\mathrm{lr}}}^L-S_{\mathrm{lr}}}{{S_{\mathrm{lr}}}^L-{S_{\mathrm{lr}}}^H}({K_{\mathrm{rl}}}^H-{K_{\mathrm{rl}}}^L)l=1,\&CenterDot;\&CenterDot;\&CenterDot;,n_p\\ n_l={n_l}^L+\frac{{S_{\mathrm{lr}}}^L-S_{\mathrm{lr}}}{{S_{\mathrm{lr}}}^L-{S_{\mathrm{lr}}}^H}({n_l}^H-{n_l}^L)l=1,\&CenterDot;\&CenterDot;\&CenterDot;,n_p\end{array}\right....(B1-14)$

8)相粘度液相粘度可根据纯组分粘度和油、水、表活剂的相浓度进行计算：

${\mathrm{\&mu;}}_l=C_{1l}{\mathrm{\&mu;}}_p{e}^{{\mathrm{\&alpha;}}_1(C_{1l}+C_{2l})}+C_{2l}{\mathrm{\&mu;}}_o{e}^{{\mathrm{\&alpha;}}_2(C_{1l}+C_{3l})}+C_{3l}{\mathrm{\&alpha;}}_3{e}^{{\mathrm{\&alpha;}}_4{C}_{1l}+{\mathrm{\&alpha;}}_5{C}_{2l}}...(B1-15)$

式中α_i为输入参数。当存在中间相聚合物时，μ_w可由μ_p代替。

3.2.2碱的物化机理

DQCHEM模型对碱在化学驱过程中的主要驱油机理包括：

1)碱与原油中的酸性组分反应产生表面活性剂；

2)碱与地层水和岩石矿物的反应，引起碱耗和离子环境的变化；

3)反应生成物和高PH值影响相态特性、界面张力、表面活性剂的吸附。

设反应体系中包含有n_i个流体组分，n_k个固体组分n_j个介质吸附的阳离子和n_m个与胶束缔合的阳离子，这些组分都是由N个基本元素组成的，可由下面的元素质量守恒方程提供的N个方程及液体相和粘土表面上的电中性方程解出：

$\left\{\begin{array}{c}{C}_n^l=\underset{j=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{h}_{\mathrm{nj}}{C}_j+\underset{k=1}{\overset{n_k}{\mathrm{\&Sigma;}}}\mathrm{gnk}{\hat{C}}_k+\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{f}_{\mathrm{ni}}{\stackrel{\&OverBar;}{C}}_i+\underset{m=1}{\overset{n_m}{\mathrm{\&Sigma;}}}{e}_{\mathrm{nm}}{\stackrel{\stackrel{\&OverBar;}{\&OverBar;}}{C}}_m\\ 0=\underset{j=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{z}_j{C}_j+\underset{m=1}{\overset{n_m}{\mathrm{\&Sigma;}}}{\stackrel{\stackrel{\&OverBar;}{\&OverBar;}}{z}}_m{\stackrel{\stackrel{\&OverBar;}{\&OverBar;}}{C}}_m\\ \mathrm{Qv}=\underset{i=1}{\overset{n_i}{\mathrm{\&Sigma;}}}{\stackrel{\&OverBar;}{z}}_i{\stackrel{\&OverBar;}{C}}_i\end{array}\right....(B2-1)$

式中Qv为所有阳离子(包括H⁺)的交换能力。

事实上对碱驱过程中的多数驱油机理描述在胶束/聚合物过程中也存在。因此，可将碱驱看成是表面活性剂驱的一种特殊情况，只是这里的表面活性剂不是人工加入、而是向地下加碱后，经物化反应在地下产生的。只要生成“表活剂”，以后的情况就和“表活剂驱”的过程相同了，故不再赘叙。

3.2.3聚合物的物化机理

聚合物溶液的高粘度能够改善油水相间的流度比，抑制注入液的突进，达到扩大波及体积，提高采收率的目的。模型对聚合物的驱油机理分以下几个方面进行描述：

1)聚合物溶液粘度

实验结果表明，在一定的剪切速率下聚合物溶液的粘度μ_p ⁰是聚合物浓度和含盐量的函数，一般情况下可用下述函数表达：

${\mathrm{\&mu;}}_p^0={\mathrm{\&mu;}}_w(1+(A_{p1}{C}_{4l}+A_{p2}{C_{4l}}^2+A_{p3}{C_{4l}}^3)C_{\mathrm{SEP}}^{\mathrm{SP}})l=\mathrm{1,3}...(B3-1)$

式中C_4l表示水相(l＝1)和微乳液相(l＝3)中聚合物的浓度；μ_w是水相粘度，A_p1、A_p2、A_p3为实验室参数；CSEP是含盐量浓度；SP为实验室参数，它决定着粘度μ_p ⁰依赖于含盐量的程度。

2)聚合物溶液流变特征

一般来说，高分子聚合物溶液都具有某种流变特征，即认为其粘度依赖于剪切速率，利用Meter方程表达这种依赖关系：

${\mathrm{\&mu;}}_p={\mathrm{\&mu;}}_w+\frac{{{\mathrm{\&mu;}}_p}^0-{\mathrm{\&mu;}}_w}{1+{\left(\frac{\mathrm{\&gamma;}}{{\mathrm{\&gamma;}}^{1/2}}\right)}^{p_{\mathrm{\&alpha;}}^{-1}}}...(B3-2)$

式中γ为多孔介质中多相流中l相的等效剪切速率，γ_1/2是 ${\mathrm{\&mu;}}_p=\frac{1}{2}({{\mathrm{\&mu;}}_p}^0+{\mathrm{\&mu;}}_w)$ 对应的剪切率，P_α是经验参数；μ_p称为聚合物溶液在多孔介质中流动的视粘度。

3)渗透率下降系数和残余阻力系数

聚合物溶液在多孔介质中渗流时，由于聚合物在岩石表面的吸附必然引起流度下降和流动阻力增加。为了表达这种现象，定义渗透率下降系数R_k和残余阻力系数R_rf分别为：

$R_k=\frac{K_w}{K_p},$

DQCHEM模型用下列公式仿真渗透率下降系数

$\left\{\begin{array}{c}{R}_k=1+\frac{(R_{k\max }-1)b_{\mathrm{rk}}{C}_{41}}{1+b_{\mathrm{rk}}{C}_{41}}\\ R_{k\max }={[1-\frac{c_{\mathrm{rk}}{\left(A_{\mathrm{pl}}{C}_{\mathrm{SEP}}^{\mathrm{SP}}\right)}^{1/3}}{{\left[\frac{\sqrt{K_x{K}_y}}{\mathrm{\&Phi;}}\right]}^{1/2}}]}^{-4}\end{array}\right....(B3-4)$

式中b_rk、c_rk为输入参数。

4)不可及孔隙体积

实验发现，孔隙介质中聚合物要比溶液中的示踪剂流动的快，这是由于聚合物的高分子结构决定了聚合物能够流经的孔隙体积小的原因。聚合物不能进入的这部分孔隙体积称为不可及孔隙体积，在DQCHEM模型中表示为：

$\mathrm{XIPV}=\frac{\mathrm{\&phi;}-{\mathrm{\&phi;}}_p}{\mathrm{\&phi;}}...(B3-5)$

式中φ为盐水测的孔隙度，φ_p为聚合物溶液测的孔隙度。

5)聚合物吸附的表征方法

渗透率多孔介质中聚合物分子的吸附滞留是由吸附于固相表面和捕集在小孔隙中造成的。这种吸附滞留与表面活性剂吸附滞留一样将使聚合物速度降低，并消耗聚合物段塞。因此，聚合物吸附可作为渗透率、含盐度及聚合物浓度的函数。

聚合物的有效盐度可表示为：

$C_{\mathrm{sep}}=\frac{C_{51}+({\mathrm{\&beta;}}_p-1)C_{61}}{C_{11}}...(B3-6)$

式中的C₅₁、C₆₁、C₁₁分别为阴离子、Ca²⁺、水在液相中的浓度；其中β由实验室测定。

3.2.4其它相关特性

注入到地下的电解质与地下流体及岩矿通常要发生阳离子交换。这种交换影响到溶液中离子的传输，并对最佳盐度和表活剂相态产生明显影响，阳离子的类型和浓度同样也影响水力传导能力，仿真计算时还要考虑这种影响。一般阳离子都以自由离子、粘土表面吸附和表活剂胶束上吸附三种形式存在，因此，Ca²⁺、Na⁺离子在粘土和表活剂上的交换守恒方程可描述为：

$\frac{{\left(C_{12}^s\right)}^2}{C_6^s}={\mathrm{\&beta;}}^S\&CenterDot;C_3^m\&CenterDot;\frac{{\left(c_{12}^f\right)}^2}{C_c^f}...(B4-1)$

$\frac{{\left(C_{12}^c\right)}^2}{C_6^c}={\mathrm{\&beta;}}^c\&CenterDot;Q_v\&CenterDot;\frac{{\left(C_{12}^f\right)}^2}{C_c^f}...(B4-2)$

而相密度可由相比重来表示：γ＝g·ρ，相比重γ_l是相压力和相组成的函数：

$\left\{\begin{array}{c}{\mathrm{\&gamma;}}_l=C_{1l}{\mathrm{\&gamma;}}_l+C_{2l}{\mathrm{\&gamma;}}_{2l}+C_{3l}{\mathrm{\&gamma;}}_{3l}+0.02533C_{6l}+C_{8l}{\mathrm{\&gamma;}}_{8l}\\ {\mathrm{\&gamma;}}_{\mathrm{kl}}={\mathrm{\&gamma;}}_{\mathrm{kR}}[1+C_k^0(P_l-P_{R0})]\end{array}\right.l=1,\&CenterDot;\&CenterDot;\&CenterDot;,n_p...(B4-3)$

式中f、c、s分别表示自由阳离子、粘土上吸附阳离子和胶束上吸附阳离子，β^c、β^S、C₃ ^m为粘土、表活剂上的交换常数及表活剂浓度，Q_v为阳离子的交换容量；γ_kR为组分K在参考压力下的比重，为输入参数，P_R0为参考压力。

3.2.5非相态稀体系的驱油机理描述

通常的三元复合驱过程中，由于经济效益因素的影响，所采用的表面活性剂浓度比较低。如果表面活性剂、水和油组成的体系中表面活性剂浓度比较低，体系就不会出现相的变化，整个体系也仅由油水两相组成，我们称这种体系为非相态稀体系。非相态稀体系三元复合驱的驱油机理通过表面活性剂、碱、油、和水之间的协同效应实现。

1)界面张力

表面活性剂、碱、油和水之间的协同效应通过界面张力活性函数描述

σ＝σ(C_S，C_A)                  (B5-1)

式中，σ是油水相间的界面张力，C_S是表面活性剂浓度，C_A是碱浓度；界面张力活性函数由实测获得；

2)毛管数

毛管数是反应由于粘性力的作用而使相剩余饱和度发生改变的一个无因次变量。毛管数的定义如下：

$N_{\mathrm{cl}}=\frac{|\stackrel{\&RightArrow;}{k}\&CenterDot;\stackrel{\&RightArrow;}{\&dtri;}{\mathrm{\&Phi;}}_{l^{\&prime;}}|}{{\mathrm{\&sigma;}}_{l{l}^{\&prime;}}}l=\mathrm{1,2}...(B5-2)$

式中，l和l′分别代表被驱替和驱替相，势梯度 表达式为：

Φ_l′＝P_l′-gρ_l′h                   (B5-3)

3)相对渗透率曲线

毛管数与相剩余饱和度之间的关系为：

$S_{\mathrm{lr}}=S_{\mathrm{lr}}^H+\frac{S_{\mathrm{lr}}^L-S_{\mathrm{lr}}^H}{1+T_l{N}_{\mathrm{cl}}}l\backslash \mathrm{1,2}...(B5-4)$

式中，T_l是常数，S_lr ^L，S_lr ^H分别是低毛管数和理想极限高毛管数下l相的剩余饱和度。

相剩余饱和度的变化必然引起相对渗透率曲线发生改变，变化后的l相对渗透率。曲线端点值k_rl和指数值 n_l可根据剩余饱和度值由低毛管数相对渗透率曲线和高毛管数相对渗透率曲线的端点值及指数值(k_lr ^L，k_rl ^H，n_l ^L，n_l ^H)进行计算。相对渗透率曲线端点值和指数值计算公式：

$k_{\mathrm{rl}}=k_{\mathrm{rl}}^L+\frac{S_{l^{\&prime;}r}^L-S_{l^{\&prime;}r}}{S_{l^{\&prime;}r}^L-S_{l^{\&prime;}r}^H}(k_{\mathrm{rl}}^H-k_{\mathrm{rl}}^l)l=\mathrm{1,2}...(B5-5)$

$n_l=n_l^L+\frac{S_{l^{\&prime;}r}^L-S_{l^{\&prime;}r}}{S_{l^{\&prime;}r}^L-S_{l^{\&prime;}r}^H}(n_l^H-n_l^L)l=l,2...(B5-6)$

4、三元复合驱仿真模型的求解与计算

如上所述，本发明的数学模型可用“高阶非线性微分方程组”来描述，因而，仿真的问题最终就成了“非线性微分方程组”的求解问题。搞过计算机数字仿真的科技工作者都知道，“非线性微分方程组”的求解问题，从来都是一项具有很大挑战性的工作，人们总是期望在所论的边界条件下，通过选择适当的算法、调整迭代参数，以便能快速、稳定地得到该“非线性微分方程组”的一组工程解，换言之，保证“迭代的收敛性、计算的稳定性与快速性和解的高精度”，一直是工程技术人员努力的目标。但是，求解这样的方程组，往往需要进行数百万、乃至上千万次的迭代循环计算，在迭代过程中如算法或迭代步长不得当，仿真计算根本就无法正常进行下去，即使能得出计算解，也不会有任何使用价值(例如单从仿真计算的意义上讲，迭代的时间步长Δt＝10^-7可能是有意义的，但对采油工程而言10^-7天就毫无实际意义)。这些问题在求解高阶微分方程组时尤其会经常遇到，因而算法问题关系到计算机仿真的成败。为了避免在工程计算时出现此类问题，我们花费了相当大的精力，开发研制出了与模型DQCHEM的求解相配套的算法及具体计算机求解程序CHEMEQCNT，取得了较好的仿真效果。概括起来，这些方法包括：计算相态体系时的显式(IMPEC)解法；计算非相态体系时利用隐式方法求解压力方程，用显式(IMPEC)法求解浓度方程(也叫等效拆分全隐式的方法)；计算相饱和度和多组分体系平衡常数时用相态闪蒸计算法，它以相平衡理论为基础，通常用以表示地层中随着压力的下降、粘度增大、密度减小、体积系数相应变化的过程。所谓IMPEC是implict pressure explicit concentration的缩写。在迭代计算中，显式求解浓度方程、隐式求解压力方程，就是已知第i步的参数，求解第i+1步的变量；隐式求解浓度方程和压力方程，就是已知第i-1步和i+1步的参数，求解第i步的变量；等效拆分全隐式法是将浓度方程等效拆分为两个相孤立的流动弥散方程，这样，假设认为相间的相互转化过程不影响组分总质量，就可以只用孤立相的过程求得总的质量，所求变量的解法采用全隐式解法。只要用某种算法语言(我们用的是Fortran77、F90、Visual Basic6.0、Visual C⁺⁺6.0)依法编成相应的计算机程序，然后由计算机反复调用、执行这些程序，就可完成模型的求解计算和计算机仿真工作。下面分别介绍求解仿真模型的具体算法及其应用。

4.1(相态)仿真模型的一般求解计算方法

由于程序CHEMEQCNT面临的是“高阶非线性微分方程组”。为了求解这类方程组，我们先从普通解法入手，然后再详细介绍本发明所用的方法，这样就会更深的理解本发明。

不失一般性，设有非线性微分方程F(x)＝0(任何形式的方程总可通过移项等变成这种形式)，其中X为N维向量，初始条件为X∈X₀；业内一般人士都知道：Newton法是解非线性方程组的最常用的方法，所以人们求解非线性方程组时首先都要选用“Newton法”；但是，用Newton法求解非线性方程组时存在着如下明显不足：

(1)用Newton法求解非线性微分方程组时要进行迭代计算，而每次迭代都要计算F′(x)；由于F′(x)是由F对向量X的n个偏导数值构造的矩阵，而在化学平衡方程组中，矩阵的每个元素表达式都很复杂，因而计算量非常巨大；

(2)仿真实践证明，在许多情况下，对初始条件X∈X₀的要求有较严格的限制，而在实际应用中给出确保收敛的迭代初值是十分困难的，给出收敛到“期望”的初值就更加困难；

(3)如果迭代过程中某一步在x_k处P(x_k)出现奇异或几乎奇异，则Newton法的计算将无法进行下去，特别是如果P(x*)在F(X)＝0解x*处奇异(类似一元代数方程产生重根)，计算就变得十分困难、复杂；

(4)而对于具有强非线性的化学平衡方程组而言，其工程解往往就在奇异点附近，因此，在求解计算过程中经常还会出现所谓“g1”，“g2”，“s1”错误(自定义错误代码)，致使迭代计算无法进行下去而停止。

为了克服Newton法的上述缺点，本发明研发了求解“非线性微分方程组”的一些改进算法，并把它们用在对“杏二西三元复合驱工业化采油试验区”的计算机仿真计算中，取得了好的效果。下面分别介绍这些方法：

(1)简化Newton法。方法的核心内容是：每进行m次Newton法迭代，更换一次F′(x^*)。对于一维问题，其几何意义是，过一点作一切线，接着作m-1次平行于该切线的割线，而不像Newton法那样过每个迭代点都作切线。用这种方法虽然所用迭代步数会有所增加，但计算精度与Newton法相同，关键是要选择一个优化的“m”，它的选择原则是：在迭代计算中，既能使所求变量的相对误差(精度)不降低，又能使迭代的次数最少。事实上m是与所论系统的规模有关的，系统的规模越大，m就应越大。实践证明，采用简化Newton法以后，在精度与Newton法基本相同的前提下，大大提高了非线性微分方程的求解效率。

(2)参数Newton法

对于一般的非线性方程组，其参数Newton法为

x^k+1＝x^k-ωF(x^k)^-1F(x^k)，k＝0，1，2，…(C-1)

其中ω是一个适当选取的参数，使‖F(x_k+1)‖≤‖F(x_k)‖-；实际计算中ω的选择是每步变化的，可以证明：对于充分靠近x*的x₀，当0＜ω＜2时参数Newton法收敛，因而可在此范围内选择ω的大小。

由于用Newton法进行仿真计算时还经常会碰到F(x*)在F(X)＝0的解x*处奇异(或近似奇异)的情形，这就有一个设法处理好奇异(或近似奇异)矩阵的问题。先把一般Newton迭代公式改写为：

x_k+1＝x_k-F(x_k)⁺F(x_k)，(k＝0，1，2，…)     (C-2)

其中[F′(X^K)]^-1为F′(X_K)的Penrose-Moore逆。

对于k＝0，1，2，…直到收敛，求解最小二乘问题，即

F(x_k)Δx_k＝-F(x_k)，                  (C-3)

x_k+1＝x_k+Δx_k

由于当F′(x*)为奇异时，F(x*)不一定为零，因而可求出满足

[F′(x*)]⁺F(x*)＝0的解x*，亦即F(x)＝0的广义解。

显然，先以ω为参量，将一般Newton迭代改为含参量公式，并进行迭代计算，当F′(x*)在F(X)＝0的解x*处奇异(或近似奇异)时，再用最小二乘法，可以得到方程F(x)＝0的广义解。这样不仅可以使方程求解顺利进行，而且大大加快了计算的速度。

4.2非相态描述模型的求解计算

由于本发明所采用的模型在保留原有相态描述功能的同时，增加了非相态驱油机理描述功能，它是根据实验室测定的表活剂、碱和界面张力活性关系图计算相间界面张力，实现模型非相态描述功能的，为了求解这样的模型，特地在计算方法中增加了隐式解法，从而使现在的模型真正做到既能对浓体系、也可以对稀体系进行计算机仿真，大大扩展了本发明的应用范围。下面简要叙述在非相态情况下在算法中增加的隐式程度。

1)非相态时的组分浓度(连续性)方程

以组分浓度的形式表示组分k的质量守恒方程为：

$\frac{\&PartialD;}{\&PartialD;t}\left(\mathrm{\&phi;}{\mathrm{\&rho;}}_k{\tilde{C}}_k\right)+\stackrel{\&OverBar;}{\&dtri;}\&CenterDot;\left[\underset{l=1}{\overset{N_p}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_k(C_{\mathrm{kl}}\stackrel{\&OverBar;}{u_l}-\stackrel{\&OverBar;}{D_{\mathrm{kl}}})\right]=R_k...(C-4)$

式中φ为孔隙度，ρ_k＝1+C_pk(P_w-P_st)，C_pk为k组分的压缩系数；P_st参考压力。扩散项表达式为Fickian形式：

${\stackrel{\&RightArrow;}{D}}_{\mathrm{kl}}=\mathrm{\&phi;}{S}_l{\stackrel{\&RightArrow;}{\stackrel{\&RightArrow;}{K}}}_{\mathrm{kl}}\&CenterDot;\stackrel{\&RightArrow;}{\&dtri;}{C}_{\mathrm{kl}}...(C-5)$

$\mathrm{\&phi;}{S}_l{\stackrel{\&RightArrow;}{\stackrel{\&RightArrow;}{K}}}_{\mathrm{kl}}=\left[\begin{array}{ccc}{K}_{\mathrm{xx}}& K_{\mathrm{xy}}& K_{\mathrm{xz}}\\ K_{\mathrm{yx}}& K_{\mathrm{yy}}& K_{\mathrm{yz}}\\ K_{\mathrm{zx}}& K_{\mathrm{zy}}& K_{\mathrm{zz}}\end{array}\right],(C-6)$

其分量可表达为

$K_{\mathrm{qr}}=\mathrm{\&phi;}{S}_l\frac{D_{\mathrm{kl}}}{\mathrm{\&tau;}}{\mathrm{\&delta;}}_{\mathrm{qr}}+\frac{({\mathrm{\&alpha;}}_L-{\mathrm{\&alpha;}}_T)u_q{u}_r}{U_l}+{\mathrm{\&alpha;}}_T{U}_l{\mathrm{\&delta;}}_{\mathrm{qr}}.(C-7)$

式中： ${\mathrm{\&delta;}}_{\mathrm{qr}}=\left\{\begin{array}{cc}1& (q=r)\\ 0& (q\&NotEqual;r)\end{array}\right.$ (q，r＝x，y，z)为Kronecker delta函数； $U_l=\sqrt{u_x^2+u_y^2+u_z^2}$ 为相速率。

利用如下所示的分裂算法，将总浓度的求解公式转化为近似求解“分相组分”浓度的相关表达式，就可进行后续的其他处理了。

2)分裂算法

从形式上看，分相组分浓度的渗流驱动方程是一个对流扩散方程，并且流场的物理过程决定了该方程以强对流为主。物理上，对流过程和扩散过程特性完全不同；对流过程和扩散过程的物理特性差别决定了我们可以采用不同的方法分别求解它们，即所谓的“分裂算法”。“分裂算法”的实质是“时间导数分步走”：先得出对流方程解，但它并非n+1时刻的解，而是中间过程n+r(0＜r＜1)时刻的对流方程解，再利用n+r时的对流方程解求出扩散方程的解，才是n+1时刻的扩散方程解，因而“分裂算法”也叫隐式法(这里我们假定求解区域可以按笛卡尔坐标离散为NX*NY*NZ个网格，即对仿真区块划分网格时符合笛卡尔坐标离散条件)。

4.3迭代时间步长的选择与更新在模型(由“非线性微分方程组”来描述)的求解计算中，一般都要经过成千上万次迭代才能完成，中间少不了对“迭代时间步长”这一参数的选择与更新。所谓“迭代时间步长t”，是指在仿真计算中多长时间进行一次迭代计算，通常这个时间步长在用户输入的数据流中给出，并且这个数值很小。在对本模型求解算法的改造中，采用了所谓“变步长技术”，因而大大提高了运算速度。迭代时间步长的更新原则是必须使迭代计算收敛，具体更新方法是：如果迭代虽然收敛，但耗时过长，就适当增大迭代时间步长，如果迭代不收敛或计算误差相对较大，即缩小时间步长，其具体大小的控制一般由数值模拟工作者的经验而定。

5、具体实施例

为了从理论与实践的结合上把本发明讲清楚，并验证DQCHEM模型的实际应用效果，我们应用“低浓度表活剂与相态结合的三元复合驱数值仿真方法”，对大庆油田“杏二西三元复合驱试验区”进行了计算机仿真及跟踪拟合。下面先把所论区块的基本情况、即原始基础数据(包括对仿真试验区的油藏地质情况的描述、油井分布方案、所用驱油材料、所采用的三元复合驱段塞注入方式与井工作制度描述等)介绍一下，然后再对仿真模型的建立、仿真数据流的生成、仿真模型(即非线性微分方程组)的求解，以及仿真结果的处理等具体仿真方法(过程)进行介绍。

5.1对杏二西三元复合驱试验区基本情况描述

杏二区三元复合驱油矿场试验位于杏树岗背斜构造西翼、北起203号断层、南至2区3排、西起204号断层、东至211号断层的比较封闭的纯油区内。试验区呈东高西低、南高北低的趋势，高差为10-20米，油层埋藏深度为800-1200米，平均为980.3米。整个试验区的面积为0.3Km²，孔隙体积43.5×10⁴m³，地质储量24.01×10⁴t，采用四注九采五点法面积井网，共有油水井13口，注入井4口，采出井9口(包括1口中心井，八口平衡井)；注采井距200米，生产井距280米。根据杏二西三元复合驱试验区实际井位分配，数值仿真模型DQCHEM对该区域的油藏地质结构的解释是：纵向上分3个小层，仿真井数为18口，仿真节点数为31×26×3＝2418个。上述驱试验区的油藏地质基础数据见表3-1；针对该试验区的油藏地质特点，分别制定了三种驱油预案，即水驱、聚合物驱、三元复合驱，进行计算机仿真。其中三元复合驱的注入材料与注入方式为：

0.275PV×1500聚合物的前置段塞；

1.0％碱+0.35％表活剂+1200ml聚合物(0.3PV)主段塞；

1.0％碱十0.15％表活剂+1200ml聚合物(0.15PV)副段塞；

0.2PV×800ml的后续保护段塞；

聚合物驱的聚合物用量为570PV.mg/L，采用单一整体段塞注入。

实际仿真时拟采用的驱油方案、所注入的驱油材料等数据描述、拟采用的驱油段塞注入程序、此驱油体系的相态组成及其参数(包括醇分配、界面张力等计算参数)分别示于表3-2、表3-3，部分相态参数示于表3-4；

表3-1：试验区基础数据表

表3-2：杏二西三元复合驱方案及实际注入情况(聚合物浓度mg/L)

实际注入化学剂量	注入量(m3)	33314	16350	162955	45576	21784	43503	21921	127562	462965
											注入(PV)	0.076	0.0376	0.351	0.105	0.05	0.1	0.05	0.29	1.05
	碱(m3)			4411.1	1413.7					5824.8
											聚合物(t)		21	332.9	78.8	23.9	34.4	17.73		508.77
	表活剂(t)			818.1	88.2					906.3

表3-3：组分组成表(相号、相名称：1、2、3；水相、油相、微乳液相)

组分号	组分	是否占有体积	浓度单位
				1	水	占	体积分数
2	油	占	体积分数
				3	表活剂	占	体积分数
4	聚合物	否	重量百分数
				5	阴离子	否	Meq/ml
6	阳离子	否	Meq/ml
				7	醇#1	占	体积分数
8	醇#2	占	体积分数
				9	示踪剂#1	否	重量百分数

表3-4：部分相态参数表

参数名称	意义
		C2PLC	II(+)型区域中褶点的油浓度
C2PLRC	II(-)型区域中褶点的油浓度
		EPSME	临界胶束浓度(CMC)——形成胶束的最小表面活性剂浓度
HBNS70	有效矿化度为0时，醇1对表活剂双结曲线影响的斜率
		HBNC70	有效含盐量为0时，醇1对表活剂双结曲线影响的截距
HBNS71	最佳矿化度为0时，醇1对表活剂双结曲线影响的斜率
		HBNC71	最佳矿化度为0时，醇1对表活剂双结曲线影响的截距
HBNC72	两倍最佳矿化度为0时，醇1对双结曲线影响的斜率
		HBNT70	在矿化度0条件下，温度对双结曲线影响的斜率
HBNT71	在两倍矿化度0条件下，温度对双结曲线影响的斜率
		HBNT72	斜率参数BT，在三相条件下与温度相关的参数

CSEL7	醇1和钙离子为0时，III型区域形成的最低含盐量浓度
		CSEU7	醇1和钙离子为0时，III型区域形成的最高含盐量浓度
BETA6	钙离子对有效含盐量影响的曲线斜率参数
		………	…………………………………………
T22	油相残余油饱和度相关的毛管数参数

根据对所述驱试验区基本情况的描述，在对杏二西三元复合驱矿场试验区块数值仿真计算中，物质平衡方程MASBLCEQ中的参量的含义与取值分别为：

n_cv为所注入三元复合体系的总组分数，数值为13，其中包括水组分、油组分、表活剂组分、聚合物组分、钙离子组分、镁离子组分等；

n_p为参与计算的相态数。当区块注入表活剂浓度小于临界胶束浓度(CMC)时为非相态计算，取n_p＝2；而当注入表活剂浓度大于临界胶束浓度时为相态计算，取n_p＝3；

ρ_k为区块原油密度，数值为0.86g/cm³。

φ为区块孔隙度，数值为0.26(无量纲)。

C_kl为各相中注入体系的组分浓度，k为组分序号，l为相的序号；

S_l为各相的剩余油饱和度；

5.2数值计算输入数据流生成

注意到，DQCHEM和CHEMEQCNT需要的输入是“仿真区块的油藏地质特征精细描述数据、驱油材组成详细体系参数”等原始数据，而其输出是能够度量驱油效果、确定驱油段塞注入程序的“流体各相界面间的张力、毛细管对油的扑集作用力、油相剩余饱和度、渗透率下降系数的大小”等物理量。由于这里输入的原始基础数据由本地提供，输出数据也是提供给本地，因而“网络化仿真服务子程SSSP”之仿真工作模式控制字WMC＝“00”；当SSSP完成系统初始化、原始数据准备好之后，就从本地数据库读入这些原始数据，并启动仿真前后数据处理软件FDEDPSW、进而关键字数据文件输入软件、变长度数组生成软件、变厚度参数与死节点处理软件，生成仿真所需格式的输入数据流，并把它们送到位于后台工作站的化学平衡反应方程组求解计算软件CHEMEQCNT之入口处，就可在SSSP的协调控制下正式开展计算机仿真计算(对模型的求解计算)。

5.3仿真过程说明

综上所述，我们结合图2给出的主程序流程，对用DQCHEM对“大庆杏二西三元复合驱试验区的计算机仿真”的过程做如下的归纳：

1)启动主、从计算机，系统仿真服务软件SSSP对系统初始化，待显示器上显出人机对话菜单后，操作人员键入具体仿真命令，由主机生成包括SWMCW在内的仿真控制字；

2)设置初始时间步长，主机根据仿真控制字(SWMCW＝00)从本地数据源获得包括仿真区块的实际油藏地质特征参数、所用驱油材料及其数量、各类钻井的类型(即是生产井还是注入井，是定压井还是定产井)及其描述参数(如否完井、射孔方向的坐标x，y，z)、注入的段塞组分与参数等在内的原始基础数据，接着用FDEDPSW之关键字数据文件输入软件、变长度数组生成软件、变厚度参数与死节点处理软件，对基础数据进行预处理，形成仿真输入数据流，送至工作站的CHEMEQCNT程序入口处；

3)根据驱油体系环境中有无表活剂注入，以及表活剂浓度的大小，用临界胶束浓度的临界值CMC对仿真驱油体系进行相态、非相态判断；当CMC≤.99就先进行如下处理，然后再执行从4)开始的程序，若无表活剂或CMC＞.99，就直接从4)开始进行仿真模型DQCHEM的求解计算：

(1)据每一口井位的地质特性及其不同的地下坐标，先把整个驱油体系用网格划分为界面特性各不相同的k个“子驱油体系块”(k＝1、2、3……)，并用k个化学平衡反应方程CHEMEQM(k)对大系统中的每个子系统进行描述；

(2)根据在实验室条件下实测得到的“表活剂、碱、聚合物与界面张力活性关系图”，计算出每个子区块的界面张力、毛管数、相对渗透率、渗透率下降系数等参数数据，送CHEMEQCNT的程序入口，并对每个子模型CHEMEQM(k)的参数进行赋值；

(3)以k为“循环控制”和“关联控制”变量，按下述步骤依次对每个化学平衡反应方程、即子模型CHEMEQM(k)进行求解计算及解的关联重组，得出总体采油方案与井工制度；

4)在SSSP的监控下，首先给相对渗透率、粘度、浓度、饱和度、毛管力等变量赋值，计算传导率、对流中的浓度项，完成压力方程PRSSUREQ系数矩阵元素的计算，然后用CHEMEQCNT对ASPEXTMDG依次进行以下求解计算：

(1)根据压力方程PRSSUREQ求解相压力、相流速和井的压力限制、井的注入和产出速率；

(2)根据物质平衡方程MASBLCEQ、浓度方程，求解出各种组分的总浓度、注入的PV数及采出程度、相比例和井筒产出液的组成；

(3)根据能量守恒方程ENGCSVEQ进行化学平衡计算，解出新的离子平衡环境；

(4)利用化学平衡模型计算有效含盐量、醇的分配系数和其它相关组分的浓度；

(5)进行“相态闪蒸计算”，求解相饱和度和相组成；

(5)利用上述结果计算界面张力、毛管数、相剩余饱和度、相对渗透率、相密度、相粘度、渗透率下降系数；

(6)对下一次迭代所需压力方程系数矩阵的元素与时间步长进行计算、动态更新后，执行从3)开始的新一轮迭代计算，直到迭代计算结果满足仿真要求；

5)将本次仿真所得数据存储于仿真结果专用缓存区BUFFER中；

6)操作人员根据人机对话菜单，重新键入仿真命令，由主机生成新的仿真控制字，仿真执行从2)开始的循环；

7)用FDEDPSW对暂存于BUFFER数据进行处理，生成“三元复合驱段塞注入程序及井工作制度文件”，并和“三元复合驱采油实验区所用采油工艺”数据进行拟合对比，然后以多种媒体形式给出对比结果，最后由仿真控制字(SWMCW控制输出仿真所得结果，供井队指导采油、科技人员阅读选用。

附图3给出了“三元复合驱采油仿真计算结果”与“实际驱采结果”之间的数据对比和二者之间的曲线拟合情形。从拟合曲线可以看出：(1)用仿真结果指导“三次采油”，可把含水率由99.9％下降至49.2％(降幅接近50％)，到2001年2月，中心井的综合含水为98.3％，累积增油11392t，提高采收率19.42％，截止到2001年6月，全区综合含水已达98％，累积增油58499t，采收率提高19.42％，说明用本发明给出的计算机仿真结果不仅对实际矿场采收有指导意义，而且效果良好；(2)如把仿真结果用于对试验区的工业开采前景进行预测，可以看出：当试验区含水为98％时，三元复合驱最终采收率可提高20％；

三元复合驱计算机仿真方法的技术效果析

由于本发明在“数学模型、数值计算方法和相关软件”三方面进行了科技创新，因而取得非常明显的技术效果。

1)技术方案先进，仿真模型适应范围广。本发明针对大庆油田的油藏地质构造的复杂性(非均匀性十分严重)的实际，给出了“网格化子模型建模法”(即先把整个驱油体系用网格划分为界面特性各不相同的k个“子驱油体系块”，然后再对每个子区块分别建立仿真模型)、“物质间的关联度法”的建模方法，针对大庆油田在低表活剂情况下较难形成相态的实际，给出了既能描述浓体系相态、也能描述稀体系非相态，既能对层间非均质、又适于层内非均质的地质模型进行描述，既适于单口井、又适于多口井的仿真，既适用于单一驱油材料、又适于多元复合驱油材料的计算机仿真模型和仿真方法，因而完全可满足各类大型油田进行油气勘探开发研究等方面的需要，应用范围非常广；

2)计算方法与仿真程序先进，确保了“三元复合驱计算机仿真”的顺利进行和“快速、精确、可靠、高效”。针对所建仿真模型为高阶非线性微分方程组因而难以求解的实际，本发明给出了包括简化Newton法、参数Newton法和专门用以求解非相态特性的“分裂算法”(即等效拆分全隐式方法)等先进算法，有效避免了用已知方法时经常出现的奇异解错误及“g1、g2、s1、s2”代码错误等，并把非线性微分方程的阶数从传统方法的14阶降到3-4阶，大大降低了对仿真模型的求解难度；在数据处理方面，开发了关键字数据文件输入软件、变长度数组生成软件、变厚度参数与死节点处理软件等，既能保证仿真的快速、准确、顺利地进行，又可方便地生成优化的“三元复合驱采油工艺及井工制文件”。实际测试表明：同已有的仿真方法相比，本发明给出的仿真方法把仿真效率提高了3-5倍；

3)仿真结果对生产实践的指导作用明显。用本发明给出的“仿真结果(后)处理与显示软件”处理仿真结果数据，可以得出与给定区块对应的驱油方案(包括井位的确定、各类井的段塞工艺、所用驱油材料及其组分等)，用该方案指导“三元复合驱工业化采油试验区”的采油实践证明：仿真结果与实采动态产油指标之间的曲线拟合良好，表明开采剩余油时的实际增油效果理想。因而，用本发明给出的方法对指导大庆油田今后的三元复合驱工业化推广的方案设计及开采、指标预测和三元复合驱效果分析有重要意义。

4)本发明所提供的计算机仿真模型、方法全面，计算技术先进，与之配套的仿真软件不仅可用于对各种类型油田的“三元复合驱采油”进行仿真，也可用于相关专业的教学，该仿真软件的通用化程度也得到大大提高。

总之，本发明在“数学模型、数值计算方法和相关软件”三个方面进行了创新，所用仿真模型功能全，适用范围广，所用数值计算方法和数据处理技术先进，大大提高了仿真的效率、精度、可靠性与仿真结果的应用范围，因而不仅对指导大庆油田今后的三元复合驱工业化推广的方案设计及开采、指标预测和三元复合驱效果分析有重要指导意义，因而它也必将为油田的勘探、开发带来巨大的经济效益和社会效益。

Claims (3)

1、一种低浓度表面活性剂与相态结合的三元复合驱计算机仿真方法，包括主、从计算机系统、计算机仿真软件包DQCHEM、计算机操作系统，主、从计算机系统包括一个位于前台的PC计算机和位于后台的工作站，以及与仿真配套的计算机外部设备，前台计算机中存放着系统仿真服务程序SSSP，通过人机对话对仿真任务进行组织及过程控制，后台工作站完成仿真过程中的各种数据处理及复杂计算任务，其特征是：所述的DQCHEM包括三元复合驱体系模型描述及生成软件ASPEXTMDG、化学平衡反应方程组求解计算软件CHEMEQCNT、仿真前后数据处理软件FDEDPSW；所述的ASPEXTMDG既能描述相态、也能描述非相态驱油体系的理化平衡反应；用DQCHEM进行三元复合驱计算机仿真，包括如下主要步骤：

1)主、从计算机启动，系统仿真服务软件对系统初始化，待显示器上显出人机对话菜单后，根据操作人员键入的具体仿真命令，主机生成仿真控制字SWMCW，并设置初始迭代时间步长；

2)主机根据仿真控制字从数据源获得所论仿真区块的实际油藏地质特征、仿真节点数、所用驱油材料在驱油体系中所占组分、各类钻井特征及所用驱油段塞注入程序等基础数据，接着由FDEDPSW之关键字数据文件输入软件、变长度数组生成软件、变厚度参数与死节点处理软件对基础数据进行预处理，形成仿真输入数据流，送至工作站的CHEMEQCNT程序入口处；

所述的FDEDPSW包括输入数据前处理软件FDPSW，计算结果后处理及显示软件EDPSW，其中FDPSW包括关键字数据文件输入软件、变长度数组生成软件，EDPSW包括一维曲线图、二维色彩图、三维立体图生成软件；

所述FDPSW对输入数据进行预处理过程是：

(1)先将输入数据文件按类型分成“题头注释和油藏数据、输出选件信息、油藏性质、一般物理性质数据、地质化学反应剂的物化参数、生物降解数据、井数据”七个向量数组，每个向量数组设一关键字，然后对关键字进行排序；

(2)根据对关键字的排序结果，依次把与关键字对应的、字长各不相同的向量数组里的数据，自动写入到“容量动态可变的数据暂存区”中，并标出各向量数组在所述存储区中的其起始地址；

(3)将带有“关键字”和“地址指针”标志且暂存于“数据暂存区”中原始基础数据，按地址指针进行连接，生成一个带有所述标志的一维大数组，然后以数据流的形式把它送至后台工作站的CHEMEQCNT程序入口处；

(4)将带有“关键字”和“地址指针”标志的一维大数组，按标志分解成与基础数据对应的向量数组，对相关数学模型的参数赋值后，用CHEMEQCNT正式进行仿真计算；

3)先完成压力方程系数矩阵元素的计算，然后根据驱油体系中有无表活剂及其浓度的大小，对仿真驱油体系进行相态、非相态判断；

4)当注入到驱油体系中的表活剂后使临界胶束浓度CMC≤某一数值时，驱油体系为非相态，则按下述“关联度建模法”建立该体系的计算机仿真数学模型；当CMC＞某一数值时，则直接从5)开始进行计算机仿真计算：

(1)根据每一口井位的地质特性及其不同的地下坐标，先把整个驱油体系用网格划分为界面特性各不相同的k个“子驱油体系块”，k＝1、2、3……，并用k个子化学平衡反应模型CHEMEQSM(k)对大系统中的每个子系统进行描述；

(2)根据在实验室条件下实测得到的“表活剂、碱、聚合物与界面张力活性关系图”，计算出每个子区块的界面张力，作为原始输入数据，经SSSP的FDPSW处理形成数据流后，送CHEMEQCNT的程序入口，并对每个子化学平衡反应模型CHEMEQSM(k)的参数进行赋值；

(3)以k为“循环控制”和“关联控制”变量，按下述步骤依次对每个化学平衡反应模型CHEMEQSM(k)进行求解计算及解的关联重组，得出总体采油方案与井工制度；

5)在SSSP的监控下，先用“仿真输入数据流”对三元复合驱仿真模型ASPEXTMDG中的相对渗透率、粘度、浓度、饱和度、毛管力变量赋值，并计算传导率、对流中的浓度项，然后用CHEMEQCNT依次进行以下求解计算：

(1)根据压力方程PRSSUREQ求解相压力、相流速和井的压力限制、井的注入和产出速率；

(2)根据物质平衡方程MASBLCEQ、浓度方程，求解出各种组分的总浓度、注入的PV数及采出程度、相比例和井筒产出液的组成；

(3)根据能量守恒方程ENGCSVEQ进行化学平衡计算，解出新的离子平衡环境；

(4)利用化学平衡模型计算有效含盐量、醇的分配系数和其它相关组分的浓度；

(5)进行“相态闪蒸计算”，求解相饱和度和相组成；

(6)利用上述结果计算界面张力、毛管数、相剩余饱和度、相对渗透率、相密度、相粘度、渗透率下降系数；

6)将本次的仿真结果存储于专用缓存区BUFFER的规定存储单元CELL(i)中；

7)判断当前仿真井的工作状态，即根据其是否完井，或是否已关井/开井，如已关井，则依次读入下一口井的类型及参数，即注入井、生产井、定压井、定产井、射孔方向(x，y，z)、井的工作状态，和井的运行参数和段塞成分，并计算下一次迭代所需压力方程系数矩阵的元素与迭代时间步长，经动态更新后，执行从4)开始的新一轮迭代计算；

8)以仿真推演时间变量T是否达到设定的最大值为判据，判断要否继续进行仿真，如需继续迭代，则更新时间变量后执行从2)开始的循环，直到迭代计算结果满足仿真要求；

9)用FDEDPSW的EDPSW对暂存于BUFFER之CELL(i)中的数据进行结果处理，生成“三元复合驱采油综合工艺文件”，然后以多种媒体的形式输出对比结果；

所述用EDPSW对计算结果进行处理的过程是：

(1)用VC⁺⁺中的类CcurveView、CtwoView、CcubeView，将暂存于BUFFER中的数据分别绘制成一维曲线图、二维色彩图、三维立体图，以便给出全井及单井的注入量、产液量、产水量、产油量等与时间的对应关系，或把二维截面上的数据显示为平面彩色图和等值线图，把三维数据转换为可随意旋转的三维立体图形，以便更直观的显示仿真结果；

(2)将暂存于BUFFER中的“仿真结果数据”和“三元复合驱采油实验区所用采油工艺”数据进行拟合对比，生成“三元复合驱段塞注入程序及井工作制度文件”；

(3)保存经EDPSW处理后的出的结果，然后或在终端设备上打印、显示处理结果，或通过网络化仿真服务软件，将处理结果送局域网上，供注册用户观看、查阅。

2、如权利要求1所述的一种低浓度表面活性剂与相态结合的三元复合驱计算机仿真方法，其特征是：用所述的三元复合驱仿真模型的求解计算软件CHEMEQCNT求解非线性微分方程时，可视情况分别采用如下方法编制相应的求解软件：

1)在用Newton法求解非线性方程组时，如此类方程组的F′(x)计算非常复杂但总是满秩矩阵，可用简化Newton法求解此类非线性方程组，即使用同一个矩阵F′(x)连续迭代m次然后再更新其矩阵元素，这里m是一个优化参量；

2)当求解相态驱油体系模型时，可利用显式解法IMPEC求解浓度方程，用闪蒸法计算相饱和度和相组成，但在求解模型时若碰到F′(x*)奇异的情形时，就在参数ω的收敛区间0＜ω＜2内，适当选择参数ω的值，并且用线性变换放大时间步长，然后再用最小二乘法求解，就可把“矩阵F′(x*)的奇异”问题避开，以使方程求解顺利而快速的进行，并且保证仿真计算的精度不降低；

3)当需要求解非相态驱油体系模型时，可利用显式法IMPEC求解浓度方程、隐式法求解压力方程，即用等效拆分全隐式方法求解非相态驱油体系模型，其要领和具体过程是：假定对仿真区块划分网格时符合笛卡尔坐标离散条件，将总浓度的求解公式转化为近似求解“分相组分”浓度的相关表达式，先求出n+r(0＜r＜1)时刻对流方程的解，再利用n+r时刻的对流方程解，求出n+1时刻扩散方程解。

3、如权利要求1所述的一种低浓度表面活性剂与相态结合的三元复合驱计算机仿真方法，其特征是：所述计算机仿真方法的编程语言环境为Fortran77或F90、Visual Basic6.0、Visual C⁺⁺6.0，仿真程序在Windows操作系统下运行，或在ALPHA、SUN、及SGI工作站的UNIX操作系统下运行。

Images