Full 3D Monopole & Dipole Open Hole Sonic Logging in Fast, Slow and Anisotropic Formation

Introduction: 

Acoustic wave propagation in multilayered structures is used in different industries for numerous measurements, evaluation and imaging applications. Technology based on this method is used in the Oil and Gas industry to characterise underground reservoirs during their entire life cycle; starting from exploration to their abandonment. This technology is referred to as Borehole Acoustic Measurement (Logging) and is further classified in to Sonic and Ultrasonic Logging

Sonic Logging Principle:

Borehole Sonic Logging is based on the principle of measuring the travel time of Compressional (P), Shear (S) and Stoneley waves propagating along a borehole with different radial depths into the formation. The operating frequency range for these measurements is typically from 0.5 to 20 kHz.  This helps to probe several meters into the formation. In this tutorial we will show you how to set up a sonic logging model and run it on the cloud.

The Model:

The model we will show you how to build today will be a fluid-filled borehole surrounded by a fast or slow formation. A formation is characterised as fast or slow depending on whether its shear wave velocity is slower or faster than the compressional wave velocity in the borehole fluid. 

The sonic measurement tool or in this case the pressure wave excitation will be located at the centre of borehole, the receivers will also be located along the centre of the borehole the first receiver will be placed 1.219m from the transmitter this will be referred to as the Tx-Rx spacing and the receivers are spaced periodically ever 0.152m the Rx-Rx spacing.

The model will be excited with both a monopole and a dipole source. Due to the size of this model as well we will need to utilise OnScales MPI capabilities on the cloud to run these models in a reasonable time. We will show the current workflow needed to setup an MPI model!

Download: 3D Sonic Logging Model

Step 1: Define variables 

We need to define the following variables:

• Borehole, tool & formation radii 
• Height of borehole and sonic tool 
• Tx-Rx and Rx-Rx spacing 
• Mesh settings

Using the height of the sonic tool and the spacing variables you can calculate the number of axial receivers. The following inputs that we use will yield a total of 24 receivers.

Below is the code - it can be copy and pasted directly into a new Analyst script.

c ****************************
c MODEL INPUTS
c ****************************

c All inputs should be in this section
symb r1 = 0.065 /* tool inner radius
symb r2 = 0.100 /* tool outer radius
symb r3 = 0.13 /* open hole radius
symb rform = 6 * $r1 /* formation radius estimate
symb h1 = 6 /* height of borehole
symb h2 = 5 /* height of sonic logging tool
symb hs = 0.5 * ( $h1 - $h2 ) /* start of the tool (equidistant from top & bottom)
symb hr = 0.285 - 0.152 /* location of last receiver on the tool
symb tr = 1.219 /* transmitter to 1st receiver spacing
symb rr = 0.152 /* receiver to receiver spacing

symb nr = ( $h2 - $tr ) / $rr /* no. of axial receivers
symb ncopy = $nr - 1

c read material properties
symb freq = 5e3
symb freqdamp = $freq

c Mesh Settings (element sizes)
symb boxx = 3.25e-3  /* element size in X & Y
symb boxz = 0.01     /* element size in Z

 Step 2: Keypoints and Grid 

It is beneficial to set up the grid using keypoints to allow for rapid parameterisation and specific control of the model. We won't go into great detail in this tutorial about keypoints but check out this article for more information on keypoints and why you should be using them.

The symbol code has a number of #key options that allow for a quick method of generating XYZ and IJK model keypoints these commands also help to parameterise the model grid generation to allow for quick changes to the model geometry to be made. We will use the #key commands #keyindx and #keycord you can find more information on these commands in this article. We will also make use of the #key coommand #keycopy which allows for quick replication of already defined keypoints we will use this to create keypoints for each of the receivers along the borehole axis.

This is exactly the same as the 2D model except we have added a 3rd dimension.

This code can be copy and pasted directly into the analyst script 

c Keypoints
/* #keycordt allows XYZ keypoints to be generated from absolute dimensions instead of the distance between keypoints like #keycord
symb #keycordt x 1 -$rform -$r3 -$r2 -$r1 0. $r1 $r2 $r3 $rform 
symb #get { idx } rootmax x     /* store interger suffix of last keypoint
symb #keyindx i 1 $idx 1 $boxx 1

symb #keycordt y 1 0. $r1 $r2 $r3 $rform 
symb #get { jdx } rootmax y     /* store interger suffix of last keypoint
symb #keyindx j 1 $jdx 1 $boxx 1

symb #keycord z 1 0 $hs $hr $rr /* bottom of borehole to first receiver 
symb #keycopy z 3 4 $ncopy      /* copy keypoints 3 & 4 out ncopy times 
symb #get { krend } rootmax z   /* store interger suffix of last keypoint
symb #keycord z $krend $z$krend $tr $hs /* finish keypoints in y starting from end of keycopy
symb #get { kdx } rootmax z     /* store interger suffix of last keypoint
symb #keyindx k 1 $kdx 1 $boxy 1


/* integer suffix to another variable 
symb indgrd = $idx 
symb jndgrd = $jdx 
symb kndgrd = $kdx 

/* define the number of nodes in the computational grid, the number of dimensions for the model, and the appropriate constraint relations.
grid $i$idx $j$jdx $k$kdx 

/* To define nodal coordinates for the model.
geom keypnt $idx $jdx $kdx 

Step 3: MPI Settings 

When we hear people talking about MPI, in most cases, this is a implementation associated with an application using MPI. It is widely used by most legacy simulation packages and aimed at those lucky enough to have High Performance Computing (HPC) hardware within their company to run their MPI application. Luckily, we can utilise the thousands of HPCs on the cloud. We are going to show you the necessary code to run an MPI model.

This code can be copy and pasted directly into the analyst script. For this code to work it needs to be inserted after the GRID command and before the GEOM command 

symb nptot = 500    /* total number of parts to partition the model into. parts = cores
symb npi = 10       /* total number of partitions in X 
symb npj = $nptot / $npi    /* calculate total number of parts in Y
symb xpi = $npi
symb xpj = $npj
symb xistep = $indgrd / $xpi 
symb xjstep = $jndgrd / $xpj

part
    max $nptot 
    symb ifr = 1
    /* Loop through parts and partition model
    do loopI I 1 $npi 1
       symb jfr = 1
       do loopJ J 1 $npj 1
       /* Calculate part number
       symb ip = $J + ( $I - 1 ) * $npj
       /* Calculate and j indices 
       symb xito = $I * $xistep
       symb ito = $xito
       symb xjto = $J * $xjstep
       symb jto = $xjto
       /* Check for last part
       if ( $I eq $npi ) then
       symb ito = $indgrd
       endif
       if ( $J eq $npj ) then
       symb jto = $jndgrd
       endif
       /* Define part
       sdiv $ip $ifr $ito $jfr $jto $k1 $kndgrd 
       symb jfr = $jto 
       end$ loopJ
    symb ifr = $ito 
    end$ loopI

Step 4: Defining Materials and assigning them to the grid 

OnScale has a material database with a number of common materials readily available for users to use however for this tutorial we will need to manually define the fast and slow formations using the material properties shown in the image above. First though we will use the database to save water to a file and we will edit that file. Save the file to the same working directory as the flxinp file 

We have created a copy of watr to act as an interface to define a pressure load when it gets to that point

Once you have saved the file open it in the OnScale application add the following lines of code starting on line 58 of that file. We have also included the anisotropic material properties for this model 

 /* fast formation
 wvsp on
 type elas 
 prop form_fast 2350 3658 2032 0.01 
 
 /* slow formation 
 wvsp on 
 type elas
 prop form_slow 2054 1890 508 0.01 
 
 /* anisotropic material properties 
 wvsp off 
 type lean 
 prop form_aniso 2350
 23.87 15.33 9.79 0 0 0 23.87
 9.79 0 0 0 15.33 0 0 
 0 2.77 0 0 2.77 0 4.27

Now that materials have been defined, we can now assign those materials to elements of our grid to do this we use the SITE command. This is where keypoints come in handy because now we can define regions between start and end keypoints in a given direction to assign material to all the elements encapsulated by that window.

For this model instead of assigning a cylindrical region of nodes to assign the fast formation material to we will assign the formation material to the full grid and assign water and steel to the appropriate cylindrical regions. Why? Because it is much easier to apply the appropriate boundary conditions to the model if the materials lie on the extents.

The model we are setting up here will take advantage of symmetry conditions on the YMIN boundaries

The code can be copy and pasted directly into the analyst script you are working with 

c **************************************************
c MATERIAL PROPERTIES
c **************************************************
symb #read sonic_logging_mat.prjmat

site
     regn form_slow                  /* All materials must be assigned a material initalise them all with void
     cyln watr stnd z z1 $z29 $x5 $x5 $x8 $x8 * * 0 360  
     cyln steel stnd z $z1 $z29 $x5 $x5 $x7 $x7 $x6 $x6 0 360 
     regn watr2 $i5 $i5+1 $j1 $j1+1 $k28-1 $k28+1
 
grph
     plot matr                       /* plot geometry  
term

Step 5: Boundary Conditions and Loads 

It is best to simplify the model as much as possible. If you are modelling a symmetrical object and the same conditions are applied to both sides of the model, there is no reason to model the entire structure. It is more efficient to model half (or less) of the structure and achieve the same results for a fraction of the computing time. We will apply a symmetry boundary on the minimum X extent of the grid. For more information on the topic of boundary conditions check out the following articles:

Boundary Conditions & Loading

How to Apply your Own Boundary Conditions

IMPD (Impedance) vs. ABSR (Absorbing) Boundary Conditions

We will also apply a pressure load on the interface of the materials watr and watr2 this will act as a monopole pulse we will also provide you with the text file that contained the time and amplitude values of the pulse we used. 

For this model we will excite it with a monopole pulse and a dipole pulse. The below code snippet contains both options, but you will need to comment out the ones you don't want to use.

c External boundary conditions
boun
    side xmin impd 2054 1890 508 
    side xmax impd 2054 1890 508
    side ymin symm
    side ymax impd 2054 1890 508
    side zmin absr
    side zmax absr
    end

c Define time varying funciton
data hist pulse * pulse_5kHz.txt    /* open the text file and save it to an array called pulse
func hist pulse                     /* assign array to a time history func 

/* Pressure Loads (monopole pulse)
plod 
 pdef pld1 func 
 matr mat1 watr2 out 
 sdef pld1 mat1 $i5 $i5+1 $j1+1 $j1+1 $k28-1 $k28+1

/* Pressure Loads (dipole pulse)
plod 
 pdef pld1 func 
 vctr vct1 1 0 0 
 sdef pld1 vct1 $i5 $i5 $j1 $j1+1 $k28-2 $k28+2     

Step 6: Calculations and Outputs 

The receivers situated on the borehole axis measure the time it takes for the acoustic pulse to travel in the formation in various forms while undergoing dispersion and attenuation. When the sound energy arrives at the receiver, having passed through the formation, it does so at different times in the form of different types of waves. This is because these waves travel with different velocities we want to capture these waveforms so we need to set up outputs to capture the acoustic pressure measured on each receiver 

c Request calculation of additional arrays
calc pres aprs                           /* open an additional APRS array - acoustic pressure 

c Request time domain outputs saved to a flxhst file
pout
     hist func                           /* write function to a output file  
     symb jrx = 3                        /* keypoint suffix of first receiver 
     do loop_rx J 1 $nr 1                /* open do loop to loop through receivers 
        hist aprs $i5 $i5 $j$jrx $j$jrx  /* loop through receiver locations to output acoustic pressure at these points
        symb jrx = $jrx + 1              /* add one to the keypoint suffix
     end$ loop_rx                        /* close loop

Step 7: Process and Model Execution 

Before running the model, you should always issue the prcs command. In preparation for processing the time-step information, it computes the model time step for stability and opens the necessary arrays for field-data storage.

The basic information required prior to executing the model:

• Model timestep (step)
• Simulation runtime (simtime)
• Simulation runtime in terms of number of timesteps (nexec)

c ***********************
c PROCESS MODEL
c ***********************
c Time settings
time * * 0.8

c Process model
prcs 

grph 
     set imag mpeg                                /* set to generate a video 
     nvew 1 
     mirr x on

c ***********************
c MODEL EXECUTION
c ***********************
symbx ncycles = 250 
symb simtime = $ncycles / $freq /* 0.01
symb nloops = 100
symb #get { step } timestep                       /* extract calculated timestep from prcs command  
symb nexec_loop = ( $simtime / $step ) / $nloops  /* partition full simulation into smaller segments for plotting

c Runtime calculations
proc exec_loop save                               /* name procedure exec_loop and save the code
exec $nexec_loop                                  /* run model partially 
grph 
     mirr x on
     plot aprs rang 0                             /* plot acoustic pressure - shows wave propagation 
     imag

end$ proc                                          /* end of procedure

c run model by calling 'plot' procedure 'nloops' times
proc exec_loop $nloops

c Save symbols to file for later use
symb #get { labl } jobname                         /* find name of run
symb #save '$labl.symb'                            /* save in symb file

c End of input file
stop

Step 8: Running the model on the cloud and downloading the results 

Before we run the simulations on the cloud, we need to enable the Advanced Cloud Job Submission option in the settings to open the settings window click  in the top right hand corner of the application window. Once the settings are open click Cloud Account Configuration and enable the setting 

Now that the models have been set up they can now  be ran on the cloud - To do this simply open the cloud scheduler (Run on Cloud )  those familiar with OnScale will notice that there is no estimate option when you open the scheduler, this is because we are running an MPI simulation. In the estimators place is the Estimate Override options all you need to do as a user is request an appropriate amount of RAM for the simulation, for these simulations a RAM request of 50GB will be sufficient, so request 50GB and click Run to submit the files to the cloud and the simulation will start.

To access the storage, click  locate the job you have just run and download the results. 

Step 8: Post-Processing the Results 

Similar to the 2D Borehole Sonic Logging Simulations we have post-processed all of our results in MATLAB if you would like to see what the results can look like if done properly see the linked 2D tutorial there you can also download the MATLAB script that we used to process these results