Buckley-Leverett two-phase problem

The Buckley-Leverett test problem is a classical reservoir simulation benchmark that demonstrates the nonlinear displacement process of a viscous fluid being displaced by a less viscous fluid, typically taken to be water displacing oil.

Problem definition

This is a simple model without wells, where the flow is driven by a simple source term and a simple constant pressure boundary condition at the outlet. We define a function that sets up a two-phase system, a simple 1D domain and replaces the default relative permeability functions with quadratic functions:

$k_{r α} (S) = min {(\frac{S - S_{r}}{1 - S_{r}}, 1)}^{n}, S_{r} = 0.2, n = 2$

In addition, the phase viscosities are treated as constant parameters of 1 and 5 centipoise for the displacing and resident fluids, respectively.

The function is parametrized on the number of cells and the number of time-steps used to solve the model. This function, since it uses a relatively simple setup without wells, uses the Jutul functions directly.

julia

using JutulDarcy, Jutul
function solve_bl(;nc = 100, time = 1.0, nstep = nc)
    T = time
    tstep = repeat([T/nstep], nstep)
    domain = get_1d_reservoir(nc)
    nc = number_of_cells(domain)
    timesteps = tstep*3600*24
    bar = 1e5
    p0 = 100*bar
    sys = ImmiscibleSystem((LiquidPhase(), VaporPhase()))
    model = SimulationModel(domain, sys)
    kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0], [0.2, 0.2])
    replace_variables!(model, RelativePermeabilities = kr)
    tot_time = sum(timesteps)
    pv = pore_volume(domain)
    irate = 500*sum(pv)/tot_time
    src  = SourceTerm(1, irate, fractional_flow = [1.0, 0.0])
    bc = FlowBoundaryCondition(nc, p0/2)
    forces = setup_forces(model, sources = src, bc = bc)
    parameters = setup_parameters(model, PhaseViscosities = [1e-3, 5e-3]) # 1 and 5 cP
    state0 = setup_state(model, Pressure = p0, Saturations = [0.0, 1.0])
    states, report = simulate(state0, model, timesteps,
        forces = forces, parameters = parameters)
    return states, model, report
end

solve_bl (generic function with 1 method)

Run the base case

We solve a small model with 100 cells and 100 steps to serve as the baseline.

julia

n, n_f = 100, 1000
states, model, report = solve_bl(nc = n)
print_stats(report)

Jutul: Simulating 1 day as 100 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │  Avg/step │  Avg/ministep │    Total │
│                │ 100 steps │ 100 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton         │      3.31 │          3.31 │  331 (0) │
│ Linearization  │      4.31 │          4.31 │  431 (0) │
│ Linear solver  │      3.31 │          3.31 │  331 (0) │
│ Precond apply  │       0.0 │           0.0 │    0 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type   │   Each │   Relative │  Total │
│               │     ms │ Percentage │      s │
├───────────────┼────────┼────────────┼────────┤
│ Properties    │ 0.0145 │     0.21 % │ 0.0048 │
│ Equations     │ 0.3319 │     6.33 % │ 0.1431 │
│ Assembly      │ 0.0068 │     0.13 % │ 0.0030 │
│ Linear solve  │ 2.7361 │    40.09 % │ 0.9057 │
│ Linear setup  │ 0.0000 │     0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │     0.00 % │ 0.0000 │
│ Update        │ 0.2619 │     3.84 % │ 0.0867 │
│ Convergence   │ 0.0115 │     0.22 % │ 0.0049 │
│ Input/Output  │ 0.0166 │     0.07 % │ 0.0017 │
│ Other         │ 3.3521 │    49.11 % │ 1.1095 │
├───────────────┼────────┼────────────┼────────┤
│ Total         │ 6.8257 │   100.00 % │ 2.2593 │
╰───────────────┴────────┴────────────┴────────╯
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │  Avg/step │  Avg/ministep │    Total │
│                │ 100 steps │ 100 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton         │      3.31 │          3.31 │  331 (0) │
│ Linearization  │      4.31 │          4.31 │  431 (0) │
│ Linear solver  │      3.31 │          3.31 │  331 (0) │
│ Precond apply  │       0.0 │           0.0 │    0 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type   │   Each │   Relative │  Total │
│               │     ms │ Percentage │      s │
├───────────────┼────────┼────────────┼────────┤
│ Properties    │ 0.0145 │     0.21 % │ 0.0048 │
│ Equations     │ 0.3319 │     6.33 % │ 0.1431 │
│ Assembly      │ 0.0068 │     0.13 % │ 0.0030 │
│ Linear solve  │ 2.7361 │    40.09 % │ 0.9057 │
│ Linear setup  │ 0.0000 │     0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │     0.00 % │ 0.0000 │
│ Update        │ 0.2619 │     3.84 % │ 0.0867 │
│ Convergence   │ 0.0115 │     0.22 % │ 0.0049 │
│ Input/Output  │ 0.0166 │     0.07 % │ 0.0017 │
│ Other         │ 3.3521 │    49.11 % │ 1.1095 │
├───────────────┼────────┼────────────┼────────┤
│ Total         │ 6.8257 │   100.00 % │ 2.2593 │
╰───────────────┴────────┴────────────┴────────╯

Run refined version (1000 cells, 1000 steps)

Using a grid with 100 cells will not yield a fully converged solution. We can increase the number of cells at the cost of increasing the runtime a bit. Note that most of the time is spent in the linear solver, which uses a direct sparse LU factorization by default. For larger problems it is recommended to use an iterative solver. The high-level interface used in later examples automatically sets up an iterative solver with the appropriate preconditioner.

julia

states_refined, _, report_refined = solve_bl(nc = n_f);
print_stats(report_refined)

Jutul: Simulating 23 hours, 60 minutes as 1000 report steps
╭────────────────┬────────────┬────────────────┬──────────╮
│ Iteration type │   Avg/step │   Avg/ministep │    Total │
│                │ 1000 steps │ 1000 ministeps │ (wasted) │
├────────────────┼────────────┼────────────────┼──────────┤
│ Newton         │      3.265 │          3.265 │ 3265 (0) │
│ Linearization  │      4.265 │          4.265 │ 4265 (0) │
│ Linear solver  │      3.265 │          3.265 │ 3265 (0) │
│ Precond apply  │        0.0 │            0.0 │    0 (0) │
╰────────────────┴────────────┴────────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type   │   Each │   Relative │  Total │
│               │     ms │ Percentage │      s │
├───────────────┼────────┼────────────┼────────┤
│ Properties    │ 0.0829 │     3.53 % │ 0.2708 │
│ Equations     │ 0.0620 │     3.44 % │ 0.2643 │
│ Assembly      │ 0.0475 │     2.64 % │ 0.2027 │
│ Linear solve  │ 2.0238 │    86.02 % │ 6.6078 │
│ Linear setup  │ 0.0000 │     0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │     0.00 % │ 0.0000 │
│ Update        │ 0.0254 │     1.08 % │ 0.0828 │
│ Convergence   │ 0.0176 │     0.98 % │ 0.0751 │
│ Input/Output  │ 0.0414 │     0.54 % │ 0.0414 │
│ Other         │ 0.0418 │     1.78 % │ 0.1364 │
├───────────────┼────────┼────────────┼────────┤
│ Total         │ 2.3527 │   100.00 % │ 7.6814 │
╰───────────────┴────────┴────────────┴────────╯
╭────────────────┬────────────┬────────────────┬──────────╮
│ Iteration type │   Avg/step │   Avg/ministep │    Total │
│                │ 1000 steps │ 1000 ministeps │ (wasted) │
├────────────────┼────────────┼────────────────┼──────────┤
│ Newton         │      3.265 │          3.265 │ 3265 (0) │
│ Linearization  │      4.265 │          4.265 │ 4265 (0) │
│ Linear solver  │      3.265 │          3.265 │ 3265 (0) │
│ Precond apply  │        0.0 │            0.0 │    0 (0) │
╰────────────────┴────────────┴────────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type   │   Each │   Relative │  Total │
│               │     ms │ Percentage │      s │
├───────────────┼────────┼────────────┼────────┤
│ Properties    │ 0.0829 │     3.53 % │ 0.2708 │
│ Equations     │ 0.0620 │     3.44 % │ 0.2643 │
│ Assembly      │ 0.0475 │     2.64 % │ 0.2027 │
│ Linear solve  │ 2.0238 │    86.02 % │ 6.6078 │
│ Linear setup  │ 0.0000 │     0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │     0.00 % │ 0.0000 │
│ Update        │ 0.0254 │     1.08 % │ 0.0828 │
│ Convergence   │ 0.0176 │     0.98 % │ 0.0751 │
│ Input/Output  │ 0.0414 │     0.54 % │ 0.0414 │
│ Other         │ 0.0418 │     1.78 % │ 0.1364 │
├───────────────┼────────┼────────────┼────────┤
│ Total         │ 2.3527 │   100.00 % │ 7.6814 │
╰───────────────┴────────┴────────────┴────────╯

Plot results

We plot the saturation front for the base case at different times together with the final solution for the refined model. In this case, refining the grid by a factor 10 gave us significantly less smearing of the trailing front.

julia

using GLMakie
x = range(0, stop = 1, length = n)
x_f = range(0, stop = 1, length = n_f)
f = Figure()
ax = Axis(f[1, 1], ylabel = "Saturation", title = "Buckley-Leverett displacement")
for i in 1:6:length(states)
    lines!(ax, x, states[i][:Saturations][1, :], color = :darkgray)
end
lines!(ax, x_f, states_refined[end][:Saturations][1, :], color = :red)
f

Example on GitHub

If you would like to run this example yourself, it can be downloaded from the JutulDarcy.jl GitHub repository as a script, or as a Jupyter Notebook

This example took 14.400415149 seconds to complete.

This page was generated using Literate.jl.

Buckley-Leverett two-phase problem ​

Problem definition ​

Run the base case ​

Run refined version (1000 cells, 1000 steps) ​

Plot results ​

Example on GitHub ​

Buckley-Leverett two-phase problem

Problem definition

Run the base case

Run refined version (1000 cells, 1000 steps)

Plot results

Example on GitHub