Dynamic Optimization of a Polymer Flooding Process Based on Implicit Discrete Maximum Principle

Abstract

Polymer flooding is one of the most important technologies for enhanced oil recovery (EOR). In this paper, an optimal control model of distributed parameter systems (DPSs) for polymer injection strategies is established, which involves the performance index as maximum of the profit, the governing equations as the fluid flow equations of polymer flooding, and some inequality constraints as polymer concentration and injection amount limitation. The optimal control model is discretized by full implicit finite-difference method. To cope with the discrete optimal control problem (OCP), the necessary conditions for optimality are obtained through application of the calculus of variations and Pontryagin’s discrete maximum principle. A modified gradient method with new adjoint construction is proposed for the computation of optimal injection strategies. The numerical results of an example illustrate the effectiveness of the proposed method.

1. Introduction

It is of increasing necessity to produce oil fields more efficiently and economically because of the ever-increasing demand for petroleum worldwide. Since most of the significant oil fields are mature fields and the number of new discoveries per year is decreasing, the use of EOR processes is becoming more and more imperative. At present, polymer flooding technology is the best method for chemically EOR [1]. It could reduce the water-oil mobility ratio and improve sweep efficiency [2–5].

Due to the high cost of chemicals, it is essential to optimize polymer injection strategies to provide the greatest oil recovery at the lowest cost. The optimization procedure involves maximizing the objective function (cumulative oil production or profit) from a polymer flooding reservoir by adjusting the injection concentration. One way of solving this problem is direct optimization with the reservoir simulator. Numerical models are used to evaluate the complex interactions of variables affecting development decisions, such as reservoir and fluid properties and economic factors. Even with these models, the current practice is still the conventional trial and error approach. In each trial, the polymer concentration of an injection well is selected based on the intuition of the reservoir engineer. This one-well-at-a-time approach may lead to suboptimal decisions because engineering and geologic variables affecting reservoir performance are often nonlinearly correlated. And the problem definitely compounds when multiple producers and injectors are involved in a field development case. The use of the optimal control method offers a way out.

The optimal control method has been researched in EOR techniques in recent years. Ramirez et al. [6] firstly applied the theory of optimal control to determine the best possible injection strategies for EOR processes. Their study was motivated by the high operation costs associated with EOR projects. The objective of their study was to develop an optimization method to minimize injection costs while maximizing the amount of oil recovered. The performance of their algorithm was subsequently examined for surfactant injection as an EOR process in a one-dimensional core flooding problem [7]. The control for the process was the surfactant concentration of the injected fluid. They observed a significant improvement in the ratio of the value of the oil recovered to the cost of the surfactant injected from 1.5 to about 3.4. Optimal control was also applied to steam flooding by Liu et al. [8]. They developed an approach using optimal control theory to determine operating strategies to maximize the economic attractiveness of steam flooding process. Their objective was to maximize a performance index which is defined as the difference between oil revenue and the cost of injected steam. Their optimization method also obtained significant improvement under optimal operation. Ye et al. [9] were involved in the study of optimal control of gas-cycling in condensate reservoirs. It was shown that both the oil recovery and the total profit of a condensate reservoir can be enhanced obviously through optimization of gas production rate, gas injection rate, and the mole fractions of each component in injection gas. Brouwer and Jansen [10], Sarma et al. [11], Zhang et al. [12], and Jansen [13] used the optimal control theory as an optimization algorithm for adjusting the valve setting in smart wells of water flooding. The water flooding scheme that maximized the profit was numerically obtained by combining reservoir simulation with control theory practices of implicit differentiation. They were able to achieve improved sweep efficiency and delayed water breakthrough by dynamic control of the valve setting.

For the previous work on optimal control of polymer flooding, Guo et al. [14] and Lei et al. [15] applied the iterative dynamic programming algorithm to solve the OCP of polymer flooding. However, the optimal control model used in their study is so simple that it is not adapt for practical oilfield development. As a result of the complicated nature of reservoir models with nonlinear constraints, it is very tedious and troublesome to cope with a large number of grid points for the state variables and control variables. To avoid these difficulties, Li et al. [16] and Lei et al. [17, 18] used the genetic algorithms to determine the optimal injection strategies of polymer flooding and the reservoir model equations were treated as a “black box.” The genetic algorithms are capable of finding the global optimum on theoretical sense, but as Sarma et al. [11] point out, they require a tens or hundreds of thousand reservoir simulation runs of very large model and are not able to guarantee monotonic maximization of the objective function.

In this paper, an optimal control model of DPS for polymer flooding is established which maximizes the profit by adjusting the injection concentration. Then, the determination of polymer injection strategies turns to solve this OCP of DPS. The model is discretized by the full implicit finite difference method. Necessary conditions for optimality are obtained by Pontryagin's discrete maximum principle. A new gradient-based numerical algorithm which makes it relatively easy to create adjoint equations is presented for solving the OCP. Finally, an example of polymer flooding project involving a heterogeneous reservoir case is investigated and the results show the efficiency of the proposed method.

2. Mathematical Formulation of Optimal Control

2.1. Performance Index

Let Ω∈𝑅² denote the domain of reservoir with boundary 𝜕Ω, let 𝑛 be the unit outward normal on 𝜕Ω, and let (𝑥,𝑦)∈Ω be the coordinate of a point in the reservoir. Given a fixed final time 𝑡_𝑓, we set Ψ=Ω×[0,𝑡_𝑓], Γ=𝜕Ω×[0,𝑡_𝑓], and suppose that there exist 𝑁_𝑤 injection wells and 𝑁_𝑜 production wells in the oilfield. The injection and production wells are located at 𝐿_𝑤={(𝑥_𝑤𝑖,𝑦_𝑤𝑖)∣𝑖=1,2,…,𝑁_𝑤} and 𝐿_𝑜={(𝑥_𝑜𝑗,𝑦_𝑜𝑗)∣𝑗=1,2,…,𝑁_𝑜}, respectively. This descriptive statement of the cost functional must be translated into a mathematical form to use quantitative optimization techniques. The oil value can be formulated as ^𝑡𝑓₀_Ω𝜉_𝑜1−𝑓_𝑤𝑞_out𝑑𝜎𝑑𝑡,(2.1) where 𝜉_𝑜 is the cost of oil per unit mass (10⁴$/m³), 𝑓_𝑤(𝑥,𝑦,𝑡) is the fractional flow of water, and 𝑞_out(𝑥,𝑦,𝑡) is the flow velocity of production fluid (m/day). We define 𝑞_out(𝑥,𝑦,𝑡)≥0 at (𝑥,𝑦)∈𝐿_𝑜 and 𝑞_out(𝑥,𝑦,𝑡)≡0 at (𝑥,𝑦)∉𝐿_𝑜.

The polymer cost is expressed mathematically as ^𝑡𝑓₀_Ω𝜉_𝑝𝑞_in𝑐_pin𝑑𝜎𝑑𝑡,(2.2) where 𝜉_𝑝 is the cost of oil per unit volume (10⁴$/m³),𝑐_pin(𝑥,𝑦,𝑡) is the polymer concentration of the injection fluid (g/L), and 𝑞_in(𝑥,𝑦,𝑡) is the flow velocity of injection fluid (m/day). We define 𝑞_in(𝑥,𝑦,𝑡)≥0 at (𝑥,𝑦)∈𝐿_𝑤 and 𝑞_in(𝑥,𝑦,𝑡)≡0 at (𝑥,𝑦)∉𝐿_𝑤.

The objective functional is, therefore, max𝐽=^𝑡𝑓₀_Ω𝜉_𝑜1−𝑓_𝑤𝑞_out−𝜉_𝑝𝑞_in𝑐_pin𝑑𝜎𝑑𝑡.(2.3)

2.2. Governing Equations

The maximization of the cost functional 𝐽 given by (2.3) is not totally free but is constrained by the system process dynamics. The governing equations of the polymer flooding process must, therefore, be developed to describe the flow of both the aqueous and oil phases through the porous media of a reservoir formation. The equations used in this paper allow for the adsorption of polymer onto the solid matrix in addition to the convective and dispersive mechanisms of mass transfer. Let 𝑝(𝑥,𝑦,𝑡),𝑆_𝑤(𝑥,𝑦,𝑡), and𝑐_𝑝(𝑥,𝑦,𝑡) denote the pressure, water saturation, and polymer concentration of the oil reservoir, respectively, at a point (𝑥,𝑦)∈Ω and a time 𝑡∈[0,𝑡_𝑓], then 𝑝(𝑥,𝑦,𝑡),𝑆_𝑤(𝑥,𝑦,𝑡), and𝑐_𝑝(𝑥,𝑦,𝑡) satisfy the following partial differential equations (PDEs).(1)The flow equation for oil phase: 𝜕𝜕𝑥𝑘^𝑝𝑟_𝑜𝜕𝑝𝜕𝑥+𝜕𝜕𝑦𝑘^𝑝𝑟_𝑜𝜕𝑝𝜕𝑦−1−𝑓^𝑤𝑞_out=ℎ𝜕𝜕𝑡𝜙1−𝑆^𝑤𝐵_𝑜,(𝑥,𝑦,𝑡)∈Ψ.(2.4)(2)The flow equation for water phase: 𝜕𝜕𝑥𝑘^𝑝𝑟_𝑤𝜕𝑝𝜕𝑥+𝜕𝜕𝑦𝑘^𝑝𝑟_𝑤𝜕𝑝𝜕𝑦+𝑞ⁱⁿ−𝑓_𝑤𝑞_out=ℎ𝜕𝜕𝑡𝜙𝐵𝑆^𝑤_𝑤,(𝑥,𝑦,𝑡)∈Ψ.(2.5)(3)The flow equation for polymer component: 𝜕𝜕𝑥𝑘^𝑑𝑟_𝑑𝜕𝑐_𝑝𝜕𝑥+𝜕𝜕𝑥𝑘^𝑝𝑟_𝑐𝜕𝑝𝜕𝑥+𝜕𝜕𝑦𝑘^𝑑𝑟_𝑑𝜕𝑐_𝑝𝜕𝑦+𝜕𝜕𝑦𝑘^𝑝𝑟_𝑐𝜕𝑝𝜕𝑦+𝑞ⁱⁿ𝑐_pin−𝑓_𝑤𝑞_out𝑐_𝑝=ℎ𝜕𝜕𝑡𝜙^𝑝𝑆_𝑤𝑐_𝑝𝐵_𝑤+𝜌_𝑟(1−𝜙)𝐶_rp,(𝑥,𝑦,𝑡)∈Ψ.(2.6)(4)The boundary conditions and initial conditions: 𝜕𝑝𝜕𝑛|||_𝜕Ω=0,𝜕𝑆^𝑤𝜕𝑛||||_𝜕Ω=0,𝜕𝑐^𝑝𝜕𝑛||||_𝜕Ω𝑝(𝑥,𝑦,0)=𝑝₀=0,(𝑥,𝑦,𝑡)∈Γ,(2.7)(𝑥,𝑦),𝑆_𝑤(𝑥,𝑦,0)=𝑆^0𝑤(𝑥,𝑦),𝑐_𝑝(𝑥,𝑦,0)=𝑐^0𝑝(𝑥,𝑦),(𝑥,𝑦)∈Ω.(2.8) where the corresponding parameters are defined as 𝑟𝑘_𝑝=𝐾ℎ,𝑘_𝑑=𝐷ℎ,_𝑜=𝑘^𝑟𝑜𝐵_𝑜𝜇_𝑜,𝑟_𝑤=𝑘^𝑟𝑤𝐵_𝑤𝑅_𝑘𝜇_𝑤,𝑟_𝑐=𝑘^𝑟𝑤𝑐_𝑝𝐵_𝑤𝑅_𝑘𝜇_𝑝,𝑟_𝑑=𝜙^𝑝𝑆_𝑤𝐵_𝑤.(2.9)The constant coefficient 𝐾(𝑥,𝑦) is the absolute permeability (𝜇m²), 𝐷(𝑥,𝑦) is the diffusion coefficient of polymer (m²/s), ℎ(𝑥,𝑦) is the thickness of the reservoir bed (m), 𝜌_𝑟 is the rock density (kg/m³), and𝜇_𝑜is the oil viscosity (mPa · s).

The oil volume factor 𝐵_𝑜, the water volume factor 𝐵_𝑤, the rock porosity 𝜙, and the effective porosity to polymer 𝜙_𝑝 are expressed as functions of the reservoir pressure 𝑝: 𝐵_𝑜=𝐵^𝑜𝑟𝐵1+𝐶_𝑜𝑝−𝑝_𝑟_𝑤,=𝐵^𝑤𝑟𝜙=𝜙_𝑟1+𝐶_𝑤𝑝−𝑝_𝑟,𝜙1+𝐶_𝑅𝑝−𝑝_𝑟_𝑝,=𝑓_𝑎𝜙,(2.10) where𝑝_𝑟is the reference pressure (MPa), 𝜙_𝑟, 𝐵_𝑜𝑟, and 𝐵_𝑤𝑟 denote the porosity, the oil and water volume factor under the condition of the reference pressure, respectively, 𝑓_𝑎 is the effective pore volume coefficient, 𝐶_𝑜, 𝐶_𝑤, and 𝐶_𝑅 denote the compressibility factors of oil, water, and rock, respectively.

Functions relating values of the oil and water relative permeabilities 𝑘_𝑟𝑜 and 𝑘_𝑟𝑤 to the water saturation 𝑆_𝑤 are 𝑘_𝑟𝑤=𝑘_{𝑟𝑤𝑟𝑜}𝑆^𝑤−𝑆_𝑤𝑐1−𝑆_𝑤𝑐−𝑆_𝑜𝑟^𝑛𝑤𝑘_𝑟𝑜,=𝑘_{𝑟𝑜𝑐𝑤}1−𝑆^𝑤−𝑆_𝑜𝑟1−𝑆_𝑤𝑐−𝑆_𝑜𝑟^𝑛𝑜,(2.11) where 𝑘_{𝑟𝑤𝑟𝑜} is the oil relative permeability at the irreducible water saturation 𝑆_𝑤𝑐, 𝑘_{𝑟𝑤𝑐𝑤} is the water relative permeability at the residual oil saturation 𝑆_𝑜𝑟, 𝑛_𝑤, and 𝑛_𝑜 are constant coefficients.

The polymer solution viscosity 𝜇_𝑝 (mPa⋅s), the permeability reduction factor 𝑅_𝑘, and the amount adsorbed per unit mass of the rock 𝐶_𝑟𝑝 (mg/g) which depend on the polymer concentration 𝑐_𝑝 are given by 𝜇_𝑝𝑅=𝜇_𝑤1+𝑎𝑝₁𝑐_𝑝+𝑎𝑝₂𝑐^2𝑝+𝑎𝑝^3𝑝_𝑘,=1+𝑅^𝑘max−1⋅𝑏_𝑟𝑘⋅𝑐_𝑝1+𝑏_𝑟𝑘⋅𝑐_𝑝𝐶_𝑟𝑝,=𝑎𝑐^𝑝1+𝑏𝑐_𝑝,(2.12) where 𝜇_𝑤 is the viscosity of the aqueous phase with no polymer (mPa · s), 𝑎𝑝₁, 𝑎𝑝₂, 𝑎𝑝₃, 𝑅_𝑘max, 𝑏_𝑟𝑘, 𝑎 (cm³/g), and 𝑏(L/g) are constant coefficients.

The fractional flow of water 𝑓_𝑤 is given by 𝑓_𝑤=𝑟^𝑤𝑟_𝑜+𝑟_𝑤.(2.13)

2.3. Constraints

Since the polymer injection amount is limited, and the negative and overhigh injection concentrations are not allowed, the constraints in polymer flooding are expressed mathematically as follows.(1)The polymer injection amount constraint: ^𝑡𝑓₀_Ω𝑞_in𝑐_pin𝑑𝜎𝑑𝑡≤𝑚_𝑝max.(2.14)(2)The injection polymer concentration constraint: 0≤𝑐_pin≤𝑐_max,(2.15) where 𝑚_𝑝max is the maximum polymer injection amount and 𝑐_max is the maximum injection polymer concentration.

3. Full Implicit Finite Difference Method

The governing equations given by (2.4)–(2.8) are nonlinear PDEs. Several finite difference approximations for the numerical simulation of such DPS are possible. We adopt a full-implicit finite-difference scheme for the calculation of the governing equations. The scheme is described below.

For two space variables, we now consider the grid system with which we divide up the reservoir region in the 𝑥-𝑦 plane. The integer 𝑖 (𝑖=1,2,…,𝑛_𝑥) is used as the index in the 𝑥-direction, and the integer 𝑗 (𝑗=1,2,…,𝑛_𝑦) for the index in the 𝑦-direction. Thus, 𝑥_𝑖 is the 𝑖th value of 𝑥, and 𝑦_𝑗 is the 𝑗th value of 𝑦. Double indexing is used to identify functions within the two-dimensional region. Let 𝐮(𝑥,𝑦,𝑡)=[𝑝,𝑆_𝑤,𝑐_𝑝]^𝑇 denote the system state vector. The discrete state vector in grid point (𝑥_𝑖,𝑦_𝑗) is described by 𝐮_𝑖,𝑗=𝐮𝑥_𝑖,𝑦_𝑗.(3.1)

The reservoir domain is divided into 𝑛_𝑥×𝑛_𝑦 blocks, and the point (𝑥_𝑖,𝑦_𝑗) is considered to be at the center of block (𝑖,𝑗). There are 𝑛_𝑥 blocks in the 𝑥-direction and 𝑛_𝑦 blocks in the 𝑦-direction. Further details concerning this grid are given in Figure 1. We identify the coordinate 𝑥_𝑖−(1/2) with the left side of the block (𝑖,𝑗), and 𝑥_𝑖+(1/2) with the right side of the block. Similarly, 𝑦_𝑗−(1/2) is identified with the bottom of the block, and 𝑦_𝑗+(1/2) with the top.

Figure 1 

Details of block-centred grid.

Using the predefined grid system, the derivatives in (2.4)–(2.6) are replaced by finite differences. Three formulas are useful in this context: 𝜕𝐮𝜕𝑡=𝐮^𝑛+1−𝐮^𝑛Δ𝑡^𝑛𝜕𝐮,𝜕𝑥=𝐮^{𝑖+(1/2),𝑗}−𝐮_{𝑖−(1/2),𝑗}𝜕𝐮Δ𝑥,𝜕𝑦=𝐮^{𝑖,𝑗+(1/2)}−𝐮_{𝑖,𝑗−(1/2)}Δ𝑦,(3.2) where 𝑛=0,1,…,𝑁−1 is the index of reservoir simulation time step, 𝑁 is the number of simulation time steps, Δ𝑡^𝑛 represents the size of the 𝑛th time step in days, 𝐮^𝑛+1 represents the state vector at time 𝑡^𝑛+1,𝑡^𝑛+1 is the total simulation time in days at the end of the 𝑛th time step (𝑡^𝑁=𝑡_𝑓), and Δ𝑥 and Δ𝑦 (m) denote the space step lengths in the 𝑥-direction and 𝑦-direction, respectively.

Apply the formulas (3.2) to discretize the governing equations (2.4)–(2.6), and multiply the grid area Δ𝑥Δ𝑦 (m²) to the two sides of the equations. For the full implicit finite difference method, the second spatial derivative items in (2.4)–(2.6) are evaluated at time 𝑡^𝑛+1 instead of at 𝑡^𝑛. Then, the full-implicit finite-difference equations of the polymer flooding model has the general form: 𝐠^𝑛=𝐅^𝑛+1𝐮̃^𝑛+1+𝐖^𝑛+1𝐮̃^𝑛+1,𝐯^𝑛−𝐀^𝑛+1𝐮̃^𝑛+1−𝐀^𝑛(𝐮̃^𝑛)=𝟎,(3.3) where 𝐠=[𝐠^𝑇1,1,…,𝐠^{𝑇𝑖,𝑗},…,𝐠^𝑇𝑛_{𝑥,𝑛𝑦}]^𝑇∈𝑅^{𝑛𝑥×𝑛𝑦×3} is the discrete governing equations of polymer flooding model, 𝐮̃=[𝐮^𝑇1,1,…,𝐮^{𝑇𝑖,𝑗},…,𝐮^𝑇𝑛_{𝑥,𝑛𝑦}]^𝑇∈𝑅^{𝑛𝑥×𝑛𝑦×3} is the discrete state vector, 𝐯=[𝑣_𝑖,𝑗∣𝑣_𝑖,𝑗=𝑐_pin𝑖,𝑗,(𝑖,𝑗)∈𝜅_𝑤]^𝑇∈𝑅^𝑁𝑤 is the control vector, 𝜅_𝑤 is the set of grid coordinates where injection wells located, 𝐅 refers to the flux terms, 𝐖 refers to the source terms, and 𝐀 refers to the accumulation terms.

The discrete governing equations𝐠_𝑖,𝑗=[𝑔_{𝑜𝑖,𝑗},𝑔_{𝑤𝑖,𝑗},𝑔_{𝑐𝑖,𝑗}]^𝑇 for the grid point (𝑖,𝑗) at time 𝑡^𝑛 are given by 𝑔𝑔^{𝑛𝑜𝑖,𝑗}=𝐹^{𝑛+1𝑜𝑖,𝑗}+𝑊^{𝑛+1𝑜𝑖,𝑗}−𝐴^{𝑛+1𝑜𝑖,𝑗}+𝐴^{𝑛𝑜𝑖,𝑗}=0,_{𝑛𝑤𝑖,𝑗}𝑔=𝐹^{𝑛+1𝑤𝑖,𝑗}+𝑊^{𝑛+1𝑤𝑖,𝑗}−𝐴^{𝑛+1𝑤𝑖,𝑗}+𝐴^{𝑛𝑤𝑖,𝑗}=0,_{𝑛𝑐𝑖,𝑗}=𝐹^{𝑛+1𝑐𝑖,𝑗}+𝑊^{𝑛+1𝑐𝑖,𝑗}−𝐴^{𝑛+1𝑐𝑖,𝑗}+𝐴^{𝑛𝑐𝑖,𝑗}=0.(3.4) Equations (3.4) are discrete flow equations for oil phase, water phase, and polymer component, respectively. The full-implicit finite-difference schemes for the flux terms 𝐅_𝑖,𝑗=[𝐹_{𝑜𝑖,𝑗},𝐹_{𝑤𝑖,𝑗},𝐹_{𝑐𝑖,𝑗}]^𝑇 are given by 𝐹^{𝑛+1𝑜𝑖,𝑗}=Δ𝑦𝑘^{𝑝𝑖+(1/2),𝑗}𝑟^𝑜𝑛+1_{𝑖+(1/2),𝑗}Δ𝑥𝑝^{𝑛+1𝑖+1,𝑗}−𝑝^{𝑛+1𝑖,𝑗}+𝑘^{𝑝𝑖−(1/2),𝑗}𝑟^𝑜𝑛+1_{𝑖−(1/2),𝑗}Δ𝑥𝑝^{𝑛+1𝑖−1,𝑗}+Δ𝑥𝑘_{𝑝𝑖,𝑗+(1/2)}−𝑝^{𝑛+1𝑖,𝑗}𝑟^𝑜𝑛+1_{𝑖,𝑗+(1/2)}Δ𝑦𝑝^{𝑛+1𝑖,𝑗+1}−𝑝^{𝑛+1𝑖,𝑗}+𝑘^{𝑝𝑖,𝑗−(1/2)}𝑟^𝑜𝑛+1_{𝑖,𝑗−(1/2)}Δ𝑦𝑝^{𝑛+1𝑖,𝑗−1}𝐹_{𝑛+1𝑤𝑖,𝑗}−𝑝^{𝑛+1𝑖,𝑗},=Δ𝑦𝑘^{𝑝𝑖+(1/2),𝑗}𝑟^𝑤𝑛+1_{𝑖+(1/2),𝑗}Δ𝑥𝑝^{𝑛+1𝑖+1,𝑗}−𝑝^{𝑛+1𝑖,𝑗}+𝑘^{𝑝𝑖−(1/2),𝑗}𝑟^𝑤𝑛+1_{𝑖−(1/2),𝑗}Δ𝑥𝑝^{𝑛+1𝑖−1,𝑗}+Δ𝑥𝑘_{𝑝𝑖,𝑗+(1/2)}−𝑝^{𝑛+1𝑖,𝑗}𝑟^𝑤𝑛+1_{𝑖,𝑗+(1/2)}Δ𝑦𝑝^{𝑛+1𝑖,𝑗+1}−𝑝^{𝑛+1𝑖,𝑗}+𝑘^{𝑝𝑖,𝑗−(1/2)}𝑟^𝑤𝑛+1_{𝑖,𝑗−(1/2)}Δ𝑦𝑝^{𝑛+1𝑖,𝑗−1}𝐹_{𝑛+1𝑐𝑖,𝑗}−𝑝^{𝑛+1𝑖,𝑗},=Δ𝑦⎡⎢⎢⎣𝑘^{𝑑𝑖+(1/2),𝑗}𝑟^𝑑𝑛+1_{𝑖+(1/2),𝑗}Δ𝑥𝑐^{𝑛+1𝑝𝑖+1,𝑗}−𝑐^{𝑛+1𝑝𝑖,𝑗}+𝑘^{𝑑𝑖−(1/2),𝑗}𝑟^𝑑𝑛+1_{𝑖−(1/2),𝑗}Δ𝑥𝑐^{𝑛+1𝑝𝑖−1,𝑗}+Δ𝑥⎡⎢⎢⎣𝑘_{𝑑𝑖,𝑗+(1/2)}−𝑐^{𝑛+1𝑝𝑖,𝑗}⎤⎥⎥⎦𝑟^𝑑𝑛+1_{𝑖,𝑗+(1/2)}Δ𝑦𝑐^{𝑛+1𝑝𝑖,𝑗+1}−𝑐^{𝑛+1𝑝𝑖,𝑗}+𝑘^{𝑑𝑖,𝑗−(1/2)}𝑟^𝑑𝑛+1_{𝑖,𝑗−(1/2)}Δ𝑦𝑐^{𝑛+1𝑝𝑖,𝑗−1}+Δ𝑦𝑘_{𝑝𝑖+(1/2),𝑗}−𝑐^{𝑛+1𝑝𝑖,𝑗}⎤⎥⎥⎦𝑟^𝑐𝑛+1_{𝑖+(1/2),𝑗}Δ𝑥𝑝^{𝑛+1𝑖+1,𝑗}−𝑝^{𝑛+1𝑖,𝑗}+𝑘^{𝑝𝑖−(1/2),𝑗}𝑟^𝑐𝑛+1_{𝑖−(1/2),𝑗}Δ𝑥𝑝^{𝑛+1𝑖−1,𝑗}+Δ𝑥𝑘_{𝑝𝑖,𝑗+(1/2)}−𝑝^{𝑛+1𝑖,𝑗}𝑟^𝑐𝑛+1_{𝑖,𝑗+(1/2)}Δ𝑦𝑝^{𝑛+1𝑖,𝑗+1}−𝑝^{𝑛+1𝑖,𝑗}+𝑘^{𝑝𝑖,𝑗−(1/2)}𝑟^𝑐𝑛+1_{𝑖,𝑗−(1/2)}Δ𝑦𝑝^{𝑛+1𝑖,𝑗−1}−𝑝^{𝑛+1𝑖,𝑗}.(3.5) The discrete source terms 𝐖_𝑖,𝑗=[𝑊_{𝑜𝑖,𝑗},𝑊_{𝑤𝑖,𝑗},𝑊_{𝑐𝑖,𝑗}]^𝑇 are expressed as 𝑊𝑊^{𝑛+1𝑜𝑖,𝑗}_{𝑛+1𝑤𝑖,𝑗}=−1−𝑓^{𝑛+1𝑤𝑖,𝑗}𝑄^{𝑛out𝑖,𝑗},𝑊_{𝑛+1𝑐𝑖,𝑗}=𝑄^𝑛_in𝑖,𝑗−𝑓^{𝑛+1𝑤𝑖,𝑗}𝑄^𝑛_out𝑖,𝑗,=𝑄^𝑛_in𝑖,𝑗𝑐^𝑛_pin𝑖,𝑗−𝑓^{𝑛+1𝑤𝑖,𝑗}𝑄^𝑛_out𝑖,𝑗𝑐^{𝑛+1𝑝𝑖,𝑗},(3.6) where 𝑄_out is the fluid production rate (𝑄_out=Δ𝑥Δ𝑦𝑞_out, m³/d, 𝑄_out𝑖,𝑗≥0 at (𝑖,𝑗)∈𝜅_𝑜, and 𝑄_out𝑖,𝑗≡0 at (𝑖,𝑗)∉𝜅_𝑜), 𝜅_𝑜 is the set of grid coordinates where production wells located, 𝑄_in is the fluid injection rate (𝑄_in=Δ𝑥Δ𝑦𝑞_in, m³/d, 𝑄_in𝑖,𝑗≥0 at (𝑖,𝑗)∈𝜅_𝑤, and 𝑄_in𝑖,𝑗≡0 at (𝑖,𝑗)∉𝜅_𝑤). The full-implicit finite-difference schemes for the accumulation terms 𝐀_𝑖,𝑗=[𝐴_{𝑜𝑖,𝑗},𝐴_{𝑤𝑖,𝑗},𝐴_{𝑐𝑖,𝑗}]^𝑇 are 𝐴^{𝑛+1𝑜𝑖,𝑗}=𝑉^𝑖,𝑗Δ𝑡^𝑛𝜙^{𝑛+1𝑖,𝑗}1−𝑆^{𝑛+1𝑤𝑖,𝑗}𝐵^{𝑛+1𝑜𝑖,𝑗},𝐴^{𝑛𝑜𝑖,𝑗}=𝑉^𝑖,𝑗Δ𝑡^𝑛𝜙^{𝑛𝑖,𝑗}1−𝑆^{𝑛𝑤𝑖,𝑗}𝐵^{𝑛𝑜𝑖,𝑗}𝐴_{𝑛+1𝑤𝑖,𝑗},=𝑉^𝑖,𝑗Δ𝑡^𝑛𝜙^{𝑛+1𝑖,𝑗}𝑆^{𝑛+1𝑤𝑖,𝑗}𝐵^{𝑛+1𝑤𝑖,𝑗},𝐴^{𝑛𝑤𝑖,𝑗}=𝑉^𝑖,𝑗Δ𝑡^𝑛𝜙^{𝑛𝑖,𝑗}𝑆^{𝑛𝑤𝑖,𝑗}𝐵^{𝑛𝑤𝑖,𝑗}𝐴_{𝑛+1𝑐𝑖,𝑗},=𝑉^𝑖,𝑗Δ𝑡^𝑛𝜙^{𝑛+1𝑝𝑖,𝑗}𝑆^{𝑛+1𝑤𝑖,𝑗}𝑐^{𝑛+1𝑝𝑖,𝑗}𝐵^{𝑛+1𝑤𝑖,𝑗}+𝜌_𝑟1−𝜙^{𝑛+1𝑖,𝑗}𝐶^{𝑛+1𝑟𝑝𝑖,𝑗}𝐴_{𝑛𝑐𝑖,𝑗},=𝑉^𝑖,𝑗Δ𝑡^𝑛𝜙^{𝑛𝑝𝑖,𝑗}𝑆^{𝑛𝑤𝑖,𝑗}𝑐^{𝑛𝑝𝑖,𝑗}𝐵^{𝑛𝑤𝑖,𝑗}−𝜌_𝑟1−𝜙^{𝑛𝑖,𝑗}𝐶^{𝑛𝑟𝑝𝑖,𝑗},(3.7) where 𝑉 is the grid volume (𝑉_𝑖,𝑗=ℎ_𝑖,𝑗Δ𝑥Δ𝑦, m³).

The boundary condition (2.7) is equivalent to 𝜕𝐮𝜕𝑥=0,𝜕𝜕𝐮𝑦=0,∀(𝑥,𝑦,𝑡)∈𝜕Ω×0,𝑡^𝑓.(3.8) Therefore, by using finite difference and (3.8), we have 𝐮𝐮_𝑖,00,𝑗=𝐮_𝑖,1,𝐮_{𝑖,𝑛𝑦+1}=𝐮_{𝑖,𝑛𝑦},𝑖=1,2,…,𝑛_𝑥,=𝐮_1,𝑗,𝐮_{𝑛𝑥+1,𝑗}=𝐮_{𝑛𝑥,𝑗},𝑗=1,2,…,𝑛_𝑦.(3.9)

The initial condition (2.8) is discretized as 𝐮^0𝑖,𝑗=𝐮𝑥_𝑖,𝑦_𝑗,0,𝑖=1,…,𝑛_𝑥,𝑗=1,…,𝑛_𝑦.(3.10)

In solving the governing equations (2.4)–(2.8) by the full-implicit finite-difference methods, a problem that was originally described by PDEs defined over continuous time and spatial domains is transformed to problem that is described by a set of nonlinear discrete algebraic equations that are implicit form.

4. Necessary Conditions of OCP for Polymer Flooding

By using the full-implicit finite-difference method, the optimal control model of polymer flooding is discretized as the general form: max𝐽=^𝑁−1_𝑛=0𝐽^𝑛=^𝑁−1_𝑛=0Δ𝑡^𝑛⎡⎢⎢⎣𝜉_𝑜(𝑖,_{𝑗)∈𝜅𝑜}1−𝑓^{𝑛+1𝑤𝑖,𝑗}𝑄^{𝑛out𝑖,𝑗}−𝜉_𝑝(𝑖,_{𝑗)∈𝜅𝑤}𝑐^𝑛_pin𝑖,𝑗𝑄^𝑛_in𝑖,𝑗s.t.𝐠_𝑛⎤⎥⎥⎦,(4.1)𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,𝑛=𝟎,𝑛𝐮̃=0,…,𝑁−1,(4.2)_𝑛||_𝑛=0=𝐮̃⁰,(4.3)_𝑁−1_𝑛=0Δ𝑡^𝑛_(𝑖,_{𝑗)∈𝜅𝑤}𝟎≤𝑐𝐯^𝑛_{pin𝑖,𝑗𝑛}𝑄^𝑛_in𝑖,𝑗≤𝑚_𝑝max,(4.4)≤𝐯_max,(4.5) where (4.1) is the discrete performance index and (4.4) is the discrete constraint of polymer injection amount.

We desire to find a set of necessary conditions for the state vector, 𝐮̃, and the control, 𝐯, to be extremals of the functional 𝐽 (4.1) subject to the governing (4.2)-(4.3) and the constraints (4.4)-(4.5). A convenient way to cope with such a discrete OCP is through the use of Pontryagin's maximum principle. The first step is to form an augmented functional by adjoining the governing equations to the performance index 𝐽: 𝐽_𝐴=^𝑁−1_𝑛=0𝐽^𝑛𝐮̃^𝑛+1,𝐯^𝑛,𝑛+𝝀^𝑛+1^𝑇𝐠^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,𝑛,(4.6) where 𝝀=[𝝀^𝑇1,1,…,𝝀^{𝑇𝑖,𝑗},…,𝝀^𝑇𝑛_{𝑥,𝑛𝑦}]^𝑇 is the adjoint vector and 𝝀_𝑖,𝑗=[𝜆_{𝑜𝑖,𝑗},𝜆_{𝑤𝑖,𝑗},𝜆_{𝑐𝑖,𝑗}]^𝑇.

We can simplify the augmented performance function of (4.6) by introducing the Hamiltonian, 𝐻, as 𝐻^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,𝝀^𝑛+1,𝑛=𝐽^𝑛𝐮̃^𝑛+1,𝐯^𝑛,𝑛+𝝀^𝑛+1^𝑇𝐠^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,𝑛.(4.7) Therefore, the augmented performance function can be written as 𝐽_𝐴=^𝑁−1_𝑛=0𝐻^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,𝝀^𝑛+1,𝑛.(4.8)

Following the standard procedure of the calculus of variations, the increment of 𝐽_𝐴, denoted by Δ𝐽_𝐴, is formed by introducing variations𝛿𝐮̃^𝑛+1,𝛿𝐮̃^𝑛,𝛿𝐯^𝑛, and𝛿𝝀^𝑛+1 giving Δ𝐽_𝐴=𝐽_𝐴𝐮̃^𝑛+1+𝛿𝐮̃^𝑛+1,𝐮̃^𝑛+𝛿𝐮̃^𝑛,𝐯^𝑛+𝛿𝐯^𝑛,𝝀^𝑛+1+𝛿𝝀^𝑛+1−𝐽_𝐴𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,𝝀^𝑛+1.(4.9) This formulation assumes that the final time, 𝑡_𝑓, is fixed.

Expanding (4.9) in a Taylor series and retaining only the linear terms gives the variation of the functional, 𝛿𝐽_𝐴, 𝛿𝐽_𝐴=^𝑁−1_𝑛=0𝜕𝜕𝐮𝐻̃^𝑛𝑛+1^𝑇𝛿𝐮̃^𝑛+1+𝜕𝐻^𝑛𝜕𝐮̃^𝑛^𝑇𝛿𝐮̃^𝑛+𝜕𝐻^𝑛𝜕𝝀^𝑛+1^𝑇𝛿𝝀^𝑛+1+𝜕𝐻^𝑛𝜕𝐯^𝑛^𝑇𝛿𝐯^𝑛.(4.10) The variations 𝛿𝐮̃^𝑛+1 and 𝛿𝐮̃^𝑛 are not independent but can be related by a discrete version of the integration by parts as ^𝑁⁻¹_𝑛=0𝜕𝜕𝐮𝐻̃^𝑛𝑛+1^𝑇𝛿𝐮̃^𝑛+1=^𝑁−1_𝑛=0𝜕𝐻𝜕𝐮̃^𝑛−1𝑛^𝑇𝛿𝐮̃^𝑛−𝜕𝐻⁻¹𝜕𝐮̃⁰^𝑇𝛿𝐮̃⁰+𝜕𝐻^𝑁−1𝜕𝐮̃^𝑁^𝑇𝛿𝐮̃^𝑁.(4.11)

By substituting (4.11) into (4.10), the first variation 𝛿𝐽_𝐴 becomes 𝛿𝐽_𝐴=^𝑁−1_𝑛=0𝜕𝐻𝜕𝐮̃^𝑛−1𝑛+𝜕𝐻^𝑛𝜕𝐮̃^𝑛^𝑇𝛿𝐮̃^𝑛−𝜕𝐻⁻¹𝜕𝐮̃⁰^𝑇𝛿𝐮̃⁰+𝜕𝐻^𝑁−1𝜕𝐮̃^𝑁^𝑇𝛿𝐮̃^𝑁+^𝑁−1_𝑛=0𝜕𝜕𝝀𝐻^𝑛𝑛+1^𝑇𝛿𝝀^𝑛+1+^𝑁−1_𝑛=0𝜕𝜕𝐻𝐯^𝑛𝑛^𝑇𝛿𝐯^𝑛+1.(4.12)

The following necessary conditions for optimality are obtained when we apply Pontryagin's discrete maximum principle.

4.1. Adjoint Equations

Since the variation 𝛿𝐮̃^𝑛 is free and not zero, the coefficient terms involving the 𝛿𝐮̃^𝑛 variation in (4.12) are set to zero. This results in the following adjoint equations: 𝜕𝐻^𝑛−1𝜕𝐮̃^𝑛+𝜕𝐻^𝑛𝜕𝐮̃^𝑛=0.(4.13) Substituting the Hamiltonian (4.7) into (4.13), the adjoint equations reduce for the polymer flooding problem under consideration as given by 𝜕𝜕𝐠𝐮̃^𝑛−1𝑛^𝑇𝝀^𝑛+𝜕𝐠^𝑛𝜕𝐮̃^𝑛^𝑇𝝀^𝑛+1+𝜕𝐽^𝑛−1𝜕𝐮̃^𝑛=0.(4.14)

4.2. Governing Equations

Since the variation 𝛿𝝀^𝑛+1 is free, we have 𝜕𝐻^𝑛𝜕𝝀^𝑛+1=𝐠^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,𝑛=0.(4.15)

4.3. Transversality Conditions

Since the initial state is specified, the variation 𝛿𝐮̃⁰ of (4.12) is zero. However, the final state is not specified; therefore, the variation 𝛿𝐮̃^𝑁 is free and nonzero. This means the following relation must be zero: 𝜕𝐻^𝑁−1𝜕𝐮̃^𝑁=0.(4.16) By substituting the Hamiltonian (4.7) into (4.16), the transversality conditions of adjoint equations for the polymer flooding OCP are expressed as 𝜕𝐽^𝑁−1𝜕𝐮̃^𝑁+𝝀^𝑁^𝑇𝜕𝐠^𝑁−1𝜕𝐮̃^𝑁=0.(4.17) Equation (4.17) shows that the adjoint variables are known at the final time 𝑡^𝑁. Since the state variables are known at the initial time and the adjoint variables are known at the final time, the discrete OCP is a split two-point boundary-value problem.

4.4. Optimal Control

With all the previous terms vanishing, the variation of the functional 𝛿𝐽_𝐴 becomes 𝛿𝐽_𝐴=^𝑁−1_𝑛=0𝜕𝜕𝐻𝐯^𝑛𝑛^𝑇𝛿𝐯^𝑛.(4.18)

When the state and control regions are not bounded, the variation of the functional must vanish at an extremal (the fundamental theorem of the calculus of variations), therefore, 𝜕𝐻^𝑛𝜕𝐯^𝑛=𝜕𝐽^𝑛𝜕𝐯^𝑛+𝜕𝐠^𝑛𝜕𝐯^𝑛^𝑇𝝀^𝑛+1=0.(4.19)

Since there are bounds and amount constraints on the control 𝐯^𝑛, we use the discrete maximum principle to assert the following necessary conditions for optimality: 𝐻^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝝀^𝑛+1,𝐯^𝑛∗=max_𝐯𝑛𝐻^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝝀^𝑛+1,𝐯^𝑛,(4.20) where 𝐯^𝑛∗ denotes the optimal control vector.

5. Numerical Solution

A gradient-based iterative numerical technique is proposed to determine the optimal injection strategies of polymer flooding. The computational procedure is based on adjusting estimates of control vector 𝐯^𝑛 to improve the value of the objective functional. For a control to be optimal, the necessary condition given by (4.20) must be satisfied. If the control 𝐯^𝑛 is not optimal, then a correction 𝛿𝐯^𝑛 is determined so that the functional is made lager, that is, 𝛿𝐽_𝐴>0. If 𝛿𝐯^𝑛 is selected as 𝛿𝐯^𝑛=𝑤𝜕𝐽^𝑛𝜕𝐯^𝑛+𝜕𝐠^𝑛𝜕𝐯^𝑛^𝑇𝝀^𝑛+1,(5.1) where𝑤 is an arbitrary positive weighting factor, the functional variation becomes 𝛿𝐽_𝐴=^𝑡𝑓₀_Ω𝑤𝜕𝐽^𝑛𝜕𝐯^𝑛+𝜕𝐠^𝑛𝜕𝐯^𝑛^𝑇𝝀^𝑛+1^𝑇𝜕𝜕𝐽𝐯^𝑛𝑛+𝜕𝐠^𝑛𝜕𝐯^𝑛^𝑇𝝀^𝑛+1𝑑𝜎𝑑𝑡≥0.(5.2) Thus, choosing 𝛿𝐯^𝑛 as the gradient direction ensures a local improvement in the objective functional, 𝐽_𝐴. At the higher and lower bounds on 𝐯^𝑛, we must make the appropriate weighting terms equal to zero to avoid leaving the allowable region.

By substituting the discrete governing equations (3.3) and the performance index (4.1) into (5.1), the required gradient for the OCP of polymer flooding is expressed as𝛿𝐯^{𝑛𝑖,𝑗}=𝑤𝜆^{𝑛+1𝑐𝑖,𝑗}−Δ𝑡^𝑛𝜉_𝑝𝑄^{𝑛in𝑖,𝑗},(𝑖,𝑗)∈𝜅_𝑤.(5.3)

The computational algorithm of control iteration based on gradient direction is as follows.

Step 1. Make an initial guess for the control vector, 𝐯^𝑛, 𝑛=0,…,𝑁−1.

Step 2. Using stored current value of 𝐯^𝑛, 𝑛=0,…,𝑁−1, integrate the governing equations forward in time with known initial state conditions. Store the discrete states 𝐮̃^𝑛+1, 𝑛=0,…,𝑁−1.

Step 3. Calculate the profit functional with the results of the forward integration.

Step 4. Solve the adjoint equations using the stored discrete states to calculate the adjoint variables 𝝀^𝑛+1,𝑛=0,…,𝑁−1, with (4.14) and (4.17). Compute and store 𝛿𝐯^𝑛 as defined by (5.3).

Step 5. Using the evaluated 𝛿𝐯^𝑛, an improved function is computed as 𝐯^𝑛_new=𝐯^𝑛_old+𝛿𝐯^𝑛,𝑛=0,1,…𝑁−1,(5.4) where 𝟎≤𝐯^𝑛_new≤𝐯_max. A single variable search strategy can be used to find the value of the positive weighting factor 𝑤 which maximizes the improvement in the performance functional using (5.4).

Step 6. The optimization algorithm is stopped when the variation 𝛿𝐯 is too small to effectively change the performance measure, that is, when ||𝐽_new−𝐽_old||<𝜀,(5.5) where 𝜀 is a small positive number.
A penalty function method is adopted to deal with the polymer injection amount constraint (4.4). For more information about that, we are able to use this method within gradient-based algorithm, see [19].
It is clear that one forward governing equations evaluation and one backward adjoint equations evaluation are required to calculate the gradients of the performance index, irrespective of the number of controls. The time required to solve the adjoint equations is of the same order of magnitude as the governing equations. Thus, with this process, a time equivalent to approximately two integrations is all that is required to calculate any number of gradients. This is why gradient-based algorithms can be very efficient and can potentially lead to huge time savings if the number of controls is large.

6. New Algorithm for Adjoint Construction

Despite the great efficiency of the gradient-based algorithm, a major drawback of the approach is that an adjoint code is required in order to apply the algorithm. The complexity of the adjoint equations is similar to that of the governing equations. A modified approach is proposed to construct the adjoint that makes it relatively easy to create the adjoint code. The approach is due to the properties of the full-implicit finite-difference code and the forms of the performance index used in the OCP model of polymer flooding.

The main coefficient terms of the adjoint equations given by (4.14) are the two Jacobians of the governing equations: 𝜕𝐠^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛𝜕𝐮̃^𝑛,𝜕𝐠^𝑛−1𝐮̃^𝑛,𝐮̃^𝑛−1,𝐯^𝑛−1𝜕𝐮̃^𝑛.(6.1) The two terms are difficult to calculate because they are functions of the governing equations. During the forward integration of the governing equations, we solve (3.3) to determine 𝐮̃^𝑛+1 at each time step. Since these equations are nonlinear with respect to 𝐮̃^𝑛+1, the usual method to solve them is through the Newton-Raphson algorithm [20]: 𝜕𝐠^𝑛𝜕𝐮̃^𝑛+1||||𝐮̃_𝑛+1,𝑘𝐮̃𝛿𝐮̃^𝑛+1,𝑘=−𝐠^𝑛𝐮̃^𝑛+1,𝐮̃^𝑛,𝐯^𝑛,(6.2)_{𝑛+1,𝑘+1}=𝐮̃^𝑛+1,𝑘+𝛿𝐮̃^𝑛+1,𝑘,(6.3) where 𝑘 is the iteration index of the Newton-Raphson algorithm at a given time step, 𝐮̃^𝑛+1,𝑘 is the state vector of the 𝑘th iteration and (𝜕𝐠^𝑛/𝜕𝐮̃^𝑛+1)|𝐮̃^𝑛+1,𝑘 is the Jacobian of the 𝑘th iteration which is calculated by 𝜕𝐠^𝑛𝜕𝐮̃^𝑛+1||||𝐮̃_𝑛+1,𝑘=𝜕𝐅^𝑛+1𝜕𝐮̃^𝑛+1,𝑘+𝜕𝐖^𝑛+1𝜕𝐮̃^𝑛+1,𝑘−𝜕𝐀^𝑛+1𝜕𝐮̃^𝑛+1,𝑘.(6.4) Equation (6.2) is linear with respect to 𝛿𝐮̃^𝑛+1,𝑘. At the convergence of the algorithm, the Jacobian used in (6.2) is the same as the second Jacobian of (6.1).

Considering the full-implicit governing equations (3.3), the first Jacobian of (6.1) is given as 𝜕𝐠^𝑛𝜕𝐮̃^𝑛=𝜕𝐀^𝑛(𝐮̃^𝑛)𝜕𝐮̃^𝑛.(6.5) The governing equations of the previous time step is expressed as 𝐠^𝑛−1=𝐅^𝑛(𝐮̃^𝑛)+𝐖^𝑛𝐮̃^𝑛,𝐯^𝑛−1−𝐀^𝑛(𝐮̃^𝑛)−𝐀^𝑛−1𝐮̃^𝑛−1=𝟎.(6.6) Then the second Jacobian for this time step is given by 𝜕𝐠^𝑛−1𝜕𝐮̃^𝑛=𝜕𝐅^𝑛−1(𝐮̃^𝑛)𝜕𝐮̃^𝑛+𝜕𝐖^𝑛−1𝐮̃^𝑛,𝐯^𝑛−1𝜕𝐮̃^𝑛−𝜕𝐀^𝑛(𝐮̃^𝑛)𝜕𝐮̃^𝑛.(6.7) The last term of (6.7) is the same as the right side of (6.5). Thus, the first Jacobian of any given time step is calculated during the computation of the second Jacobian of the previous time step.

The rest terms of adjoint equations are the derivatives of the scalar performance index 𝐽^𝑛−1 with respect to the state vector 𝐮̃^𝑛. Due to the discrete form of performance index (4.1), the terms 𝜕𝐽^𝑛−1/𝜕𝐮̃^𝑛 are directly functions of the source terms derivatives with respective to states 𝜕𝐖^𝑛/𝜕𝐮̃^𝑛. This is also calculated within the forward integration as seen from (6.4) and is easy to extract.

Thus, all the terms required for solving the adjoint equations can be calculated during the forward governing equations evaluation itself. The computational algorithm of control iteration based on gradient direction is modified as follows.

Step 1. Make an initial guess for the control vector, 𝐯^𝑛, 𝑛=0,…,𝑁−1.

Step 2. Using stored current value of 𝐯^𝑛, 𝑛=0,…,𝑁−1, integrate the governing equations forward in time with known initial state conditions. Store the discrete states 𝐮̃^𝑛+1and the coefficients involved in the adjoint equations 𝜕𝐠^𝑛/𝜕𝐮̃^𝑛, 𝜕𝐠^𝑛/𝜕𝐮̃^𝑛+1, 𝜕𝐖^𝑛+1/𝜕𝐮̃^𝑛+1, 𝑛=0,…,𝑁−1.

Step 3. Calculate the profit functional with the results of the forward integration.

Step 4. Solve the adjoint equations using the stored Jacobians and source terms derivatives in Step 2 to calculate the adjoint variables 𝝀^𝑛+1, 𝑛=0,…,𝑁−1, with (4.14) and (4.17). Compute and store 𝛿𝐯^𝑛 as defined by (5.3).

Step 5. Using the evaluated 𝛿𝐯^𝑛, an improved function is computed as (5.4).

Step 6. The optimization algorithm is stopped when the variation 𝛿𝐯 is too small to effectively change the performance measure.

It should be noted that mathematically there is no change in the algorithm, but the adjoint equations become much simpler to code. Furthermore, the adjoint equations code can be fully consistent with the governing equations if the polymer flooding model is changed. For example, some terms reflecting new physical characteristics are added to (3.3). This is because the coefficient terms in adjoint equations are taken directly from the governing equations in the proposed method.

7. Case Study

In this section, we present a numerical example of optimal control for polymer flooding done with the proposed implicit discrete maximum principle method.

The two-phase flow of oil and water in a heterogeneous two-dimensional reservoir is considered. The reservoir covers an area of 421.02 × 443.8 m² and has a thickness of 5 m and is discretized into 90(9×10×1) grid blocks. The production model is a five-spot pattern, with one production well P1 located at the center of the reservoir (5,6) and four injection wells W1–W4 placed at the four corners (1, 10), (9,10), (1, 1), and (9,1) as shown in the permeability distribution map of Figure 2. Polymer is injected when the fractional flow of water for the production well comes to 97% after water flooding. The time domain of polymer injection is 0–1440 days and the polymer flooding project life is 𝑡_𝑓=5500 (days). Figures 3 and 4 show the contour maps of the initial water saturation 𝑆^0𝑤 and the initial reservoir pressure 𝑝⁰, respectively. The initial polymer concentration is 𝑐^0𝑝=0(g/L). In the performance index calculation, we use the price of oil𝜉_𝑜=0.0503(10⁴$/m³) [80($/bbl)], and the cost of polymer 𝜉_𝑝=2.5×10⁻⁴(10⁴$/kg). The fluid rate of the production well is 𝑄_out=60 m³/day, and the fluid rate of every injection well is 𝑄_in=15 m³/day. The PDEs are solved by full-implicit finite-difference method with step size 10 days. For the constraint (2.15), the maximum injection polymer concentration is 𝑐_max=2.2(g/L). The parameters of the reservoir description and the fluid data are shown in Tables 1 and 2, respectively.

Table 1 

Parameters of reservoir description used in the example.

Table 2 

Fluid data used in the example.

Figure 2 

Permeability (μm2) distribution.

Figure 3 

Initial water saturation contour map.

Figure 4 

Initial reservoir pressure (MPa) contour map.

The polymer injection strategies obtained by the conventional engineering judgment method (trial and error) are the same 1.8(g/L) for all injection wells. The performance index is 𝐽=$1.592×10⁷ with oil production 32429 m³ and polymer injection 155520 kg. For comparison, the results obtained by engineering judgment method are considered as the initial control strategies of the proposed implicit discrete maximum principle method. The maximum polymer injection amount in constraint (2.14) is 𝑚_𝑝max=155520 (kg). A backtracking search strategy [19] is used to find the appropriate weighting term 𝑤 and the stopping criterion is chosen as 𝜀=1×10⁻⁵. By using the proposed algorithm, we obtain a cumulative oil of 33045 m³ and a cumulative polymer of 155520.01 kg yielding the profit of 𝐽^∗=$1.623×10⁷ over the polymer flooding project life of the reservoir. The results show an increase in performance index of $3.1×10⁵. Figures 5, 6, 7, and 8 show the optimal control policies of the injection wells W1–W4. As a result, the optimal injection polymer concentration profiles of W1, W2 are significantly different from those of W3, W4. It is mainly due to the differences of the well positions and the distance to the production well, as well as the reservoir heterogeneity and the uniform initial water saturation distribution. Figures 9 and 10 show the fractional flow of water in production well and the cumulative oil production curves of the two methods, respectively. It is obvious that the fractional flow of water obtained by implicit discrete maximum principle method is lower than that by engineering judgment. Therefore, with the same cumulative polymer injection, the proposed method gets more oil production and higher recovery ratio.

Figure 5 

Injection polymer concentration of well W1.

Figure 6 

Injection polymer concentration of well W2.

Figure 7 

Injection polymer concentration of well W3.

Figure 8 

Injection polymer concentration of well W4.

Figure 9 

Fractional flow of water for the production well P1.

Figure 10 

Cumulative oil production.

8. Conclusion

In this work, a new optimal control model of DPS is established for the dynamic injection strategies making of polymer flooding. Necessary conditions of this OCP are obtained by using the calculus of variations and Pontryagin's discrete maximum principle. A gradient-based iterative computational algorithm with new adjoint construction is proposed for the determination of optimal injection strategies. It is shown that the modified approach makes it relatively easy to code the adjoint equations as compared to the standard approach. The optimal control model of polymer flooding and the proposed method are used for a reservoir example and the optimum injection concentration profiles for each well are offered. The results show that the profit is enhanced by the proposed method. Meanwhile, more oil production and higher recovery ratio are obtained. And the injection strategies chosen by engineering judgment are the same for all the wells, whereas the optimal control policies by the proposed method are different from each other as a result of the reservoir heterogeneity and the uniform initial conditions.

In conclusion, given the properties of an oil reservoir and the properties of a polymer solution, optimal polymer flooding injection strategies to maximize profit can be designed by using implicit discrete maximum principle. The approach used is a powerful tool that can aid significantly in the development of operational strategies for EOR processes.

Acknowledgments

This work was supported by the Natural Science Foundation of China (60974039), the Natural Science Foundation of Shandong Province of China (ZR2011FM002), the Fundamental Research Funds for the Central Universities (27R1105018A), and the Postgraduate Innovation Funds of China University of Petroleum (CX-1245). The authors would like to express their great appreciation to the editors and anonymous referees for their helpful and constructive comments and suggestions, which have helped to improve this paper.