Triaxial Dense no Peak

Asked by Joao Augusto Chueire on 2021-06-09

Hello everyone,

I am a PhD Student at INSA de lyon and I am making a rectangular triaxial simulation of a real experiment that exists here.

I have managed to create a simulation following the many examples in this forum and the triax-tutorial/session1 file. But I am having problems simulating a dense soil behavior (hardening - peak - softening). I tried the following methods to create a dense specimen:

1 - Insertion through cloud with 0 friction angle and then apply radius expansion,
2 - Insertion through cloud with a high friction angle (45) and then apply radius expansion. Next gradual reduction of the friction angle (from 45 to 0.1) with piston activation to reduce the porosity
3 - Insertion through cloud with a high friction angle (45) and then apply radius expansion. Next gradual reduction of the friction angle (from 45 to 0.1) with radius expansion to reduce the porosity

In all three cases I cannot get a peak stress, it seems like the soil is not loose but also not dense enough. In the 2 and 3 cases I manage to reduce the porosity from 0.45 to 0.4 .

In this Google drive folder you will be able to find the script and also a Image of the behavior of my simulations.
https://drive.google.com/drive/folders/1Fdw4PusTu-Li1z45kz1JrEfQeYGyga69?usp=sharing

PS: in the image the porosity is not updated after the start of the real triaxial test thus leading to a constant porosity.
PS 2 : I am a little too obssessed with automating things, so I hope my script does not confuse you.

Thank you very much for your time,

GONCALVES DE O C Joao
<email address hidden>

SCRIPT ::

# encoding: utf-8
#
# Made by GONCALVES DE OLIVEIRA CHUEIRE Joao Augusto as part of PhD at INSA de Lyon
#
# The initial packing composed of a cloud of specified granulometry (radii and percentages)
# The triaxial consists of 3 stages:
#
# 1. Calibration : after grains insertion their size is increased untill sigStar value
# is reached to create the sample asked.
# 2. Isotropic consolidation : piston advancing until sigCons is reached in all directions;
# 3. Triaxial test : Constant-strain compression along the z-axis, while maintaining constant
# stress (sigCons) laterally. Run untill a certain vertical strain level is reached;
#
# Controlling of strain and stresses is performed via TriaxialStressController,
# of which parameters determine type of control and also stability condition .
# The changes between each phase is made with the checkState function and the
# variable stt_dict.
#

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Variables #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Geometry
w=0.10 # specimen targeted widith in 0 direction
l=0.35 # specimen targeted length in 1 direction
h=0.45 # specimen targeted height in 2 direction

# Material Parameters
 #grain
gDsty = 2600 # grain density
gYoung = 1e9 # grain young modulus
gPRat = 0.3 # grain poisson ratio
gFrAng = 30 # grain friciton angle
 #wall
wDsty = 0 # wall density
wYoung = 50e9 # wall young modulus
wPRat = 0.3 # wall poisson ratio
wFrAng = 0.0 # wall friciton angle

# Simulation Parameters - [convention : positive is traction]
 #insertion parameters
insT ='NEW' # insertion type : 'NEW' or 'LOAD'
tPo = 0.40 # targeted porosity
rfz = 0.33 # radius delta factor
nPart = 3000 # nb of particles
 #triaxial parameters
sigStar =-90e3 # calibration stress
sigCons =-100e3 # consolidation stress
sigComp =-100e3 # compressison stress
thk = 0.001 # piston thickness
lRate =-0.05 # compression strain rate
dplSpd = 0.1 # piston maximal displacement speed
eMax = 0.15 # strain value to stop execution
 #general
Damp = 0.2 # damp coef

# Save path - Matlab & VTK
svPath='generalvariables/'
vtkPath='vtk/'

#acceptable errors to pass to next phase
uFR = 0.03 # unbalanced force on Radius Increase phase
uFS = 0.01 # unbalanced force on Calibration phase
uFC = 0.001 # unbalanced force on Consolidation phase
stS = 5000 # stress on Calibration phase
stC = 500 # stress on Consolidation phase
mnI = 20e3 # minimal number of timesteps to nextphase

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Start files #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Auto generated variables
startT = O.time # time calculation
stt_dict={'val': 0} # state variable used on stateCheck function

# Imports
from datetime import date
import errno
from yade import pack, qt, plot,os,timing

#Create folder
try:
 os.mkdir(svPath)
except OSError as exc:
 if exc.errno != errno.EEXIST:
  raise
 pass
try:
 os.mkdir(vtkPath)
except OSError as exc:
 if exc.errno != errno.EEXIST:
  raise
 pass

#Create StressForce files
fnm=svPath+'servoforce.txt'
f=open(fnm,'w')
f.write('Yade triaxial simulation - Piston files \n')
f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
f.write('Step\tSx\tSy\tSz\tUz\tUy1\tUy2\tUx1\tUx2 \n' )
f.close()

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Materials and geometries #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# first choose a high gFrAng0
if gFrAng>45:
 gFrAng0=float(gFrAng)
else:
 gFrAng0=45.

#Create materials for spheres and plates
SphereMat = O.materials.append(FrictMat(young = gYoung, poisson = gPRat,
 frictionAngle = radians(gFrAng0), density = gDsty))
WallMat = O.materials.append(FrictMat(young = wYoung, poisson = wPRat,
 frictionAngle = radians(wFrAng), density = wDsty))

if insT=='NEW':
 # Create box - box are created before spheres, thus the triax_ID is not redefined
 mn,mx = Vector3(-w/2,0.,0.),Vector3(w/2,l,h)
 wallIds=O.bodies.append(aabbWalls([mn,mx],thickness=thk,oversizeFactor=2,material=WallMat))
 #Changing following creating type
 sp=pack.SpherePack()
 sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=-1,rRelFuzz=rfz,
  num=nPart,periodic=False)
 sp.toSimulation(material = SphereMat)
elif insT=='LOAD':
 #Try to load file
 loadNm=insT.upper()+'ins.yade.bz2'
 try:
  sp.load(loadNm)
 except:
  print('#**#** File not found in the root of the execution script **#**#')
  exit()
else:
 print('#**#** New/Load Preference not set **#**#')
 exit()

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Engines #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Triaxial Controller
triax1= TriaxialStressController(
 # specify target values and whether they are strains or stresses
 goal1=sigStar,
 goal2=sigStar,
 goal3=sigStar,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is zyx, i.e.
 ## "001" (1) means "x", "011" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask=7,
 thickness = thk,
 max_vel=dplSpd, # velocity
 wall_back_activated = False, #turn off movment on the wall direction 2(z) movment +
 internalCompaction = True, #If true the confining pressure is generated by growing particles
 maxMultiplier=1.+1.5e5/gYoung, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+4e3/gYoung,
 computeStressStrainInterval = 10,
 )

# Contact friction
newton=NewtonIntegrator(damping=Damp)
O.dt = 0.25*PWaveTimeStep()

#Stuff?
O.usesTimeStepper=True
O.timingEnabled=True
#O.trackEnergy=True # enable energy tracking in the code

# Start engine
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],
 ),
 #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.4),
 triax1,
 newton,
 PyRunner(iterPeriod=1000,command='addPlotData()'),
 PyRunner(iterPeriod=1000,command='cmdUpdate()'),
 PyRunner(iterPeriod=2000,command='checkState()'),
 PyRunner(iterPeriod=2000,command='saveFiles()'),
 #VTKRecorder(iterPeriod=5000,fileName=vtkPath+'vt',recorders=['spheres','bstresses','coordNumber','mass','velocity'])
 #TriaxialStateRecorder(iterPeriod=2000,file=svPath+'WallStresses')
]

# Add plots
plot.plots={'i':('ub',),'i ':('s1','s2','s3'),' i':('e1','e2','e3'),' i ':('por'),
 # energy plot' i ':(O.energy.keys,None,'totE'),
}
# Plot
plot.plot()

if insT=='LOAD':
 # If load jump directly to the consolidation
 prepareCalibration()
 prepareConsolidation()

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Main Functions #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def checkState():
 global stt_dict,uFR,uFS,uFC,stS,stC,mnI,emax,gFrAng0,Porosity
 # Function that controls the triaxial engine behavior in three phases
 if stt_dict['val']==0:
  # Diam change : increase the diameter of the particles to create a sample
  # where are particles are in contact. Using a high gFrAngle0
  unb=unbalancedForce()
  strXY =(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2
  strZ=triax1.stress(triax1.wall_front_id)[2]
  if unb<uFR and abs(sigStar-strXY)<stS and abs(sigStar-strZ)<stS and O.iter>mnI:
   prepareCalibration()
 if stt_dict['val']==1:
  # Calibration : gradually lower the value of gFrAng0 until the target porosity is reached
  unb=unbalancedForce()
  if unb<uFS:
   if Porosity<=tPo:
    prepareConsolidation()
   else:
    if gFrAng0==0.20:
     prepareConsolidation()
     print('Targeted porosity could not be reached, continuing with p=' + str(Porosity))
     return
    elif gFrAng0>5.:
     gFrAng0=5.
    elif gFrAng0==5.:
     gFrAng0=1.
    elif gFrAng0==1.:
     gFrAng0=0.20
    O.materials[0].frictionAngle=radians(gFrAng0)
    print('Ang='+str(gFrAng0)+' Porosity :'+ str(Porosity))
    setContactFriction(radians(gFrAng0))
    stt_dict.update({'consoSt':O.iter})
 elif stt_dict['val']==2:
  # Consolidation : insotropic control untill all pistons are applying the same
  # stress to the grains.
  unb=unbalancedForce()
  strXY =(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2
  strZ=triax1.stress(triax1.wall_front_id)[2]
  if unb<uFC and abs(sigCons-strXY)<stC and abs(sigCons-strZ)<stC and (O.iter-stt_dict['consoSt'])>mnI:
   prepareCompression()
 else:
  # Compression : piston direction 2 control changed to strian. Run untill the
  # maximal strain defined above is applyed.
  e2=-triax1.strain[2]
  if e2>eMax:
   O.save(svPath+insT.upper()+'compEnd.yade.bz2')
   stt_dict.update({'compEnd':O.iter})
   print('#**#** End Compression at '+ str(stt_dict['compEnd'])+' **#**#')
   # Prepare to end simulation
   #Save plot data
   plot.saveDataTxt(svPath+'PlotData.txt.bz2')
   #Save last files
   saveFiles()
   #create file for Matlab
   endFile()
   O.pause()

def prepareCalibration():
 global stt_dict
 # Count nb of spheres
 nbSph=len(O.bodies)-6
 # Save variables for later usage
 stt_dict.update({'val': 1,'consoSt':O.iter,'nbS':nbSph})
 # add gravity
 newton.gravity = (0.0,0.0,-9.81)
 # Prind in cmd line
 print('#**#** End Diametre change at '+ str(O.iter)+' **#**#')

def prepareConsolidation():
 global stt_dict
 # Prepare for consolidation
 # Save state
 O.save(svPath+insT.upper()+'ins.yade.bz2')
 # Save variables for later usage
 stt_dict.update({'val': 2,'consoSt':O.iter,'Pz':O.bodies[5].state.pos[2],
  'Py1':O.bodies[3].state.pos[1],'Py2': O.bodies[2].state.pos[1],
  'Px1':O.bodies[1].state.pos[0],'Px2':O.bodies[0].state.pos[0]})
 # Change piston states
 triax1.goal1 = sigCons
 triax1.goal2 = sigCons
 triax1.goal3 = sigCons
 triax1.depth0 = triax1.depth
 triax1.height0 = triax1.height
 triax1.width0 = triax1.width
 triax1.internalCompaction = False
 #update material values
 O.materials[0].frictionAngle=radians(gFrAng)
 setContactFriction(radians(gFrAng))
 # restart plot data
 #if insT!='LOAD':
  #plot.saveDataTxt(svPath+'InsertionPlotData.txt.bz2')
  #plot.resetData()
 # Prind in cmd line
 print('#**#** End Calibration at '+ str(stt_dict['consoSt'])+' **#**#')
 print('w: '+str(triax1.depth)+'\t l: '+str(triax1.height)+'\t h: '+ str(triax1.width))

def prepareCompression():
 global stt_dict
 # Save variables
 O.save(svPath+insT.upper()+'consoEnd.yade.bz2')
 stt_dict.update({'val': 3,'consoEnd':O.iter})
 # Change piston states
 triax1.goal1 = sigComp
 triax1.goal2 = sigComp
 triax1.goal3 = lRate
 triax1.stressMask=3 # integer for binary 011
 # Prind in cmd line
 print('#**#** End Consolidation at '+ str(stt_dict['consoEnd'])+' **#**#')

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Execution Data #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def cmdUpdate():
 sxy=(-triax1.stress(triax1.wall_right_id)[0]-triax1.stress(triax1.wall_top_id)[1])/2
 print('I: %d\t ubf: %.4f\tsxy: %.0f\tsz: %.0f\tez: %.3f'%(O.iter,unbalancedForce(),
  sxy,-triax1.stress(triax1.wall_front_id)[2],-triax1.strain[2]))

# Plots preparation
def addPlotData():
 global Porosity
 #mStress = -(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1]+
 # triax1.stress(triax1.wall_front_id)[2])/3.0
 #Vol = ((O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*
 # (O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*
 # (O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2]))
 # Volume spheres
 Vsph = sum(4.0/3.0*pi*b.shape.radius**3 for b in O.bodies if isinstance(b.shape,Sphere))
 Porosity = 1-Vsph/(h*w*l)
 plot.addData(e1=-triax1.strain[0], e2=-triax1.strain[1], e3=-triax1.strain[2],
  ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
  s1=-triax1.stress(triax1.wall_right_id)[0],
  s2=-triax1.stress(triax1.wall_top_id)[1],
  s3=-triax1.stress(triax1.wall_front_id)[2],
  ub=unbalancedForce(),
  i=O.iter,
  por=Porosity
  #totE=O.energy.total(),**O.energy
  #dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
  #disx=triax1.width,
  #disy=triax1.height,
  #disz=triax1.depth,
  #porosity=Porosity,
  #ke=utils.kineticEnergy(),
  )

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Save Functions #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

def saveGrains():
 # Save grains details into a txt file for matlab to read
 global stt_dict
 # per grain stress calculation following https://yade-dem.org/doc/user.html#micro-stress-and-micro-strain
 TW=TesselationWrapper()
 TW.setState()
 TW.computeVolumes()
 stress=bodyStressTensors()
 # start writing
 fnm=svPath+'grains'+str(O.iter)+'.txt'
 f=open(fnm,'w')
 f.write('Yade triaxial simulation - Grains files for step'+str(O.iter)+'\n')
 f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
 f.write('Number of grains : '+str(stt_dict['nbS'])+'\n')
 f.write('ID\tX\tY\tZ\tRadius\tUx\tUy\tUz\tSx\tSy\tSz\tSxy\tSxz\tSyz \n' )
 for k in O.bodies:
  if isinstance(k.shape, Sphere):
   stp = stress[k.id]*4.*pi/3.*k.shape.radius**3/TW.volume(k.id)
   [ID,r,pos,vel]=[k.id,k.shape.radius,k.state.pos,k.state.vel]
   f.write(str(ID)+'\t'+str(pos[0])+'\t'+str(pos[1])+'\t'+str(pos[2])+
    '\t'+str(r)+'\t'+str(vel[0])+'\t'+str(vel[1])+'\t'+str(vel[2])+
    '\t'+str(stp[0,0])+'\t'+str(stp[1,1])+'\t'+str(stp[2,2])+
    '\t'+str((stp[0,1]+stp[1,0])/2)+'\t'+str((stp[0,2]+stp[2,0])/2)+
    '\t'+str((stp[1,2]+stp[2,1])/2)+'\n'
    )
 f.close()

def saveContacts():
 # Save contacts details into a txt file for matlab to read
 global stt_dict
 fnm=svPath+'contForce'+str(O.iter)+'.txt'
 f=open(fnm,'w')
 f.write('Yade triaxial simulation - Contact files step'+str(O.iter)+'\n')
 f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
 f.write('Number of grains : '+str(stt_dict['nbS'])+'\n')
 f.write('ID1\tID2\tFx\tFy\tFz \n' )
 for k in O.interactions:
  [ID1,ID2,force]=[k.id1,k.id2,k.phys.normalForce+k.phys.shearForce]
  f.write(str(ID1)+'\t'+str(ID2)+'\t'+str(force[0])+'\t'+str(force[1])+'\t'+str(force[2])+'\n')
 f.close()

def savePiston():
 # Save pistons details into a txt file for matlab to read
 fnm=svPath+'servoforce.txt'
 f=open(fnm,'a')
 f.write(str(O.iter)+'\t'+str(triax1.stress(triax1.wall_right_id)[0])+'\t'+str(triax1.stress(triax1.wall_top_id)[1])+'\t'+
  str(triax1.stress(triax1.wall_front_id)[2])+'\t'+str(O.bodies[5].state.pos[2])+'\t'+
  str(O.bodies[3].state.pos[1])+'\t'+str(O.bodies[2].state.pos[1])+'\t'+
  str(O.bodies[1].state.pos[0])+'\t'+str(O.bodies[0].state.pos[0])+'\n'
  )
 f.close()

def saveFiles():
 # Command to lauch all three save files in one pyrunner.
 # Only launch after the calibration part
 if stt_dict['val']>1:
  saveGrains()
  saveContacts()
  savePiston()

def endFile():
 # Print file with important information for matlab
 global stt_dict
 print('Time taken:',(O.time - startT))
 #YADEtoMatlab file
 f=open('YADEtoMatlab.txt','w')
 f.write('The following values give the Matlab app information about the simulaiton\n')
 f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
 f.write('Number of grains : '+str(stt_dict['nbS'])+'\n')
 f.write(' \n')
 f.write('ExpWidth = '+str(triax1.width0)+'\n' )
 f.write('ExpLength = '+str(triax1.height0)+'\n' )
 f.write('ExpHeight = '+str(triax1.depth0)+'\n' )
 f.write('Timestep = '+str(O.dt)+'\n' )
 f.write('StartCons = '+str(stt_dict['consoSt'])+'\n' )
 f.write('EndCons = '+str(stt_dict['consoEnd'])+'\n' )
 f.write('EndComp = '+str(stt_dict['compEnd'])+'\n' )
 f.write('vtkFolder = '+svPath[:-1]+'\n' )
 f.write('genFolder = '+svPath[:-1]+'\n' )
 f.write('PositionPz = '+str(stt_dict['Pz'])+'\n' )
 f.write('PositionPy1 = '+str(stt_dict['Py1'])+'\n' )
 f.write('PositionPy2 = '+str(stt_dict['Py2'])+'\n' )
 f.write('PositionPx1 = '+str(stt_dict['Px1'])+'\n' )
 f.write('PositionPx2 = '+str(stt_dict['Px2'])+'\n' )
 f.write(' \n')
 f.close()

############################# Auto Run #############################
#O.run(); O.wait()

#For info :
# %F.Pt => F stand for flag (nb of characters), P stand for precision (nb of val after .) and t for the
#type of data (d signed integer decimal, floating point decimal, e floating point exponential format)
#f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16d %-16g\n'%(pos[0],pos[1],pos[2],rad,displ[0],displ[1],displ[2],b.id,b.material.density))

Question information

Language:
English Edit question
Status:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jérôme Duriez (jduriez) said :
#1

Hi,

If I were you, I would try to get back known DEM results showing that dense behavior [*].

By adopting the same particle size distribution (at least), and checking whether your generation procedures enable you to get similar initial porosities than those previous DEM simulations.

[*] https://onlinelibrary.wiley.com/doi/full/10.1002/nag.2774, Figure 3 in case you would appreciate an example. Dense-like behavior is obtained therein, under 10 kPa confining pressure, for a 0.38 initial porosity

Can you help with this problem?

Provide an answer of your own, or ask Joao Augusto Chueire for more information if necessary.

To post a message you must log in.