How to Set Up the Fluidity Parameter γ of the Hoek-Brown Model with the Softening Model
How to Set Up the Fluidity Parameter γ of the Hoek-Brown Model with the Softening Model
How to Set Up the Fluidity Parameter γ of the Hoek-Brown Model with the Softening Model
Fluidity Parameter γ of
the Hoek-Brown Model
with the Softening Model
Richard Witasse
Principal Product Manager,
Geotechnical,
Bentley Systems
Introduction
In the framework of numerical analyses of geomaterials, one of the classical problems for modeling the
development of shear bands is the pathological mesh-dependence of the computed solution which
implies failure without energy dissipation. To avoid this unphysical behavior, an internal length must be
introduced to govern the evolution of the shear band thickness in the post-peak regime of the material
response.
More specifically, it is well known that introducing too large a difference between the dilatancy angle
𝜓𝜓𝑚𝑚𝑚𝑚𝑚𝑚 and the equivalent friction angle 𝜑𝜑𝑒𝑒𝑒𝑒 will lead to a mesh-dependent solution which then will be
almost systematically associated with:
The numerical solution then becomes ill-posed and the conditioning of the numerical solution further
worsens as the difference between the dilatancy angle and the friction angle increases or as the element
size in which shear bands develop as a result of shear failure becomes smaller.
In this context, the Hoek-Brown with softening (HBS) model uniquely proposes a feature to restore
the mesh-objectivity of the numerical solution through a visco-plastic regularization based on the
over-stress theory of Perzyna, thus enabling the introduction of an internal length through a temporal
gradient. The use of the viscous regularization technique in the Hoek-Brown model with softening
is a definite advantage when one wants to consider a high level of non-associativity and maintained
mesh-objective solution as shear plastic strains further develop.
• A loading rate which is created by associating a physical time duration in the calculation phase
this technique is introduced.
• An extra model parameter called the fluidity γ (which is the inverse of a viscosity) that needs to
be properly calibrated in relation to the previously introduced rate of loading.
The idea is to simulate a circular tunnel excavation in an infinite Hoek-Brown medium subjected to an
in-situ isotropic stress field 𝜎𝜎0 . Using a given set of HBS model parameters, this simple excavation
analysis will be run multiple times, each time with a different value for the fluidity γ (all other model
parameters remaining constant, of course). Once all runs are terminated, the radial stress evolution
over the radius can be plotted and the influence of the fluidity parameter in terms of stress
development can be properly evaluated and quantified.
Most specifically, we are interested in the computation of the stress evolution over the plastic radius
surrounding the excavated area. The plastic radius can be evaluated from the closed form solution
provided by Carranza-Torres and Fairhurst (1999).
𝜎𝜎0 𝑠𝑠 𝑝𝑝𝑖𝑖 𝑠𝑠
𝑆𝑆0 = + 2 and 𝑃𝑃𝑖𝑖 = + 2
𝑚𝑚𝑏𝑏 𝜎𝜎𝑐𝑐𝑐𝑐 𝑚𝑚𝑏𝑏 𝑚𝑚𝑏𝑏 𝜎𝜎𝑐𝑐𝑐𝑐 𝑚𝑚𝑏𝑏
Where 𝜎𝜎0 is the initial stress in the rock and 𝑝𝑝𝑖𝑖 the support pressure onto the excavation circular
contour. Then the scaled critical internal pressure, 𝑃𝑃𝑖𝑖𝑐𝑐𝑐𝑐 , at which the elastic limit of the stress state
is reached, can be calculated as:
1 2
𝑃𝑃𝑖𝑖𝑐𝑐𝑐𝑐 = �1 − �1 + 16𝑆𝑆0 �
16
𝑠𝑠
𝑝𝑝𝑖𝑖𝑐𝑐𝑐𝑐 = �𝑃𝑃𝑖𝑖𝑐𝑐𝑐𝑐 − � 𝑚𝑚𝑏𝑏 𝜎𝜎𝑐𝑐𝑐𝑐
𝑚𝑚𝑏𝑏2
A plastic region developing uniformly around the hole will develop only if 𝑝𝑝𝑖𝑖𝑐𝑐𝑐𝑐 is larger than the
supporting pressure. The extent of the failure zone is:
𝑟𝑟 = �𝐴𝐴⁄𝜋𝜋
R has to be sufficiently large to avoid boundary effects. For instance, R = 5 rpl could be accounted for.
pi
Project properties will define the model contour ranging from xmin = 0 to xmax = R. For the vertical model
extent, one can define ymin = -0.1 to ymax = 0.
One soil borehole containing a single layer ranging from 0 down to -0.1 m could then be introduced.
Place the water level at an arbitrary low level to ensure the soil formation will remain dry (for instance,
Head = -1 m).
Soil material properties are defined using the user-defined soil model HB-GSM or HB-SSM contained in
the hbs64.dll. The drainage type will be Drained. It is important to set both saturated and unsaturated
unit weight to zero as previously explained. Model parameters should then be defined by the user.
Note that each fluidity value that the user is willing to investigate should lead to the definition of a
new material set (such that there are as many material sets as fluidity values to be considered).
In the Structures mode, it is important to partition the geometry and introduce geometry lines at:
Boundary conditions are defined by prescribing fixed vertical direction over the bottom and top
horizontal lines and fixed horizontal displacement at x = R. A distributed load should also be applied
at x = r to represent the (optional) supporting pressure (see Figure 2a).
(a) Geometry
(b) Mesh
In the calculation mode, initial stresses are defined during the initial phase by means of field stress with
a uniform state of stress σ1 = σ2 = σ3 = -σ0. Then, as many plastic analyses should be independently
created (from the initial phase) as there are fluidity parameters (and, therefore, independent material
sets) to be investigated (see Figure 1). For each calculation phase, it is important to:
• Activate for each calculation phase the prescribed displacement along with the distributed load
if non-zero.
It is also very important to associate to each calculation phase a duration Δt = 10 days that will define
the rate of loading. Ultimately, the loading rate effect will solely be a function of γ×Δt. Still, we prefer
to give to Δt a reasonable physical meaning, although it would not be necessarily required, as finally
γ×Δt = (10γ) × (Δt/10)!
Finally, by setting to just 1 the maximum number of cores to be used for each phase (numerical control
parameters → max core to use), then independent calculation phases could be run simultaneously,
resulting in an even shorter overall calculation time.
Once the calculation has terminated, one can enter PLAXIS Output, only displaying the plastic area
cluster (by hiding the other soil cluster from x = rpl to x = R) and visualize the Cartesian Total Stress σxx
(radial stress) as shown in Figure 4.
Then one can draw a cross section from x = r to x = rpl and then apply the cross section tool that will
require the selection of the last steps of each plastic calculation phase (as shown in Figure 5).
Results obtained in this context are displayed in Figure 6 from which it can be evaluated the effect of
introducing fluidity onto the stress development around the excavated zone.
Figure 6: Radial stress evolution in the plastic zone as a function of the fluidity parameter γ.
The process documented in the previous section has been fully automated in a Python script called
fluidityVariationStudy.py. All functions being called are presented in the Appendix.
fluidityVariationStudy.py
from plxscripting.easy import *
from getInputAxiTunnel import param
from tunnelAxiModel import createAndRunModel
from postProcessFluidity import displayResults
password = "9>a>LNFyt4CGvcWh"
discretization = 50
#
# Check if no fluidity is defined. If not add it
if 0 not in param['Fluidity']:
param['Fluidity'].append(0)
param['Fluidity'].sort()
# Check if plasticity will develop (otherwise solution remains elastic and there is no
point doing sensitivity analysis)
if param['pcr_i'] < param['supportPressure']:
sys.exit("The solution of the excavation problem will remain elastic!")
g_o.close()
Figure 7: Automated presentation of the radial stress evolution in the plastic zone as a function of the
fluidity parameter γ.
References
Carranza-Torres, C., and Fairhurst, C. (1999). “The elasto-plastic response of underground excavations in
rock masses that satisfy the Hoek-Brown failure criterion” - International Journal of Rock Mechanics and
Mining Sciences, 36(6), 777-809.
The input parameters are summarized in a dictionary param being defined by the user in
getInputAxiTunnel.py.
getInputAxiTunnel.py
import math
param = {}
# USER_INPUT
param['eqRadius'] = 6.
param['supportPressure'] = 0.
param['excavationTime'] = 10.
# Material properties
param['E'] = 7.8e6
param['v'] = 0.3
param['k'] = -1
param['sig_ci'] = 80e3
param['m_i'] = 18
param['GSI'] = 49
param['D'] = 0.2
param['alpha'] = 1
param['m_psi0'] = 0.
param['GSI_res'] = 49
param['F_psi'] = 1
#param['m_res'] = 9
#param['s_res'] = 1
#param['m_psi_res'] = 1
#param['B_mb'] = 9
#param['B_s'] = 1
#param['B_psi'] = 1
param['m_b'] = param['m_i']*math.exp((param['GSI']-100.)/(28.-14.*param['D']))
param['s'] = math.exp((param['GSI']-100.)/(9.-3.*param['D']))
param['S_0'] = param['sigIni']/(param['m_b']*param['sig_ci']) +
param['s']/(param['m_b']*param['m_b']) # scaled far-field stress
param['P_i'] = param['supportPressure']/(param['m_b']*param['sig_ci']) +
param['s']/(param['m_b']*param['m_b']) # scaled internal pressure
param['pcr_i'] = (param['Pcr_i'] -
param['s']/(param['m_b']*param['m_b']))*param['m_b']*param['sig_ci'] # critical
internal pressure
param['b_pl'] = param['eqRadius']*math.exp(2.*(math.sqrt(param['Pcr_i'])-
math.sqrt(param['P_i']))) # extend of the failure zone
def HBSParam(param):
i = 0
matSetList = []
if param['softeningType'] == "GSM":
("Colour", 16761441),
("SoilModel", 100),
("UserDLLName", "hbs64.dll"),
("UserModel", "HB-GSM"),
("DrainageType", 0),
("gammaUnsat", 0.),
("gammaSat", 0.),
("User1", param['E']),
("User2", param['v']),
("User3", param['D']),
("User4", param['GSI']),
("User5", param['alpha']),
("User6", param['m_i']),
("User7", param['m_psi0']),
("User8", param['GSI_res']),
("User9", param['sig_ci']),
("User10", param['B_GSI']),
("User11", param['F_psi']),
("User12", param['Fluidity'][i]),
("EoedRef", param['E']),
("cref", 1),
("phi", 30),
("K0Primary", 1),
("K0Secondary", 1),
("UDPower", 0),
("UDPRef", 100)]
else:
("Colour", 16761441),
("SoilModel", 100),
("UserDLLName", "hbs64.dll"),
("UserModel", "HB-SSM"),
("DrainageType", 0),
("gammaUnsat", 0.),
("gammaSat", 0.),
("User1", param['E']),
("User2", param['v']),
("User3", param['D']),
("User4", param['GSI']),
("User5", param['alpha']),
("User6", param['m_i']),
("User7", param['m_psi0']),
("User8", param['m_res']),
("User9", param['s_res']),
("User10", param['m_psi_res']),
("User11", param['sig_ci']),
("User12", param['B_mb']),
("User13", param['B_s']),
("User14", param['B_psi']),
("User16", param['k']),
("EoedRef", param['E']),
("cref", 1),
("phi", 30),
("K0Primary", 1),
("K0Secondary", 1),
("UDPower", 0),
("UDPRef", 100)]
i = i + 1
matSetList.append(matset)
return matSetList
soilLayerThickness = 0.1
s.new()
# Project properties
g.Project.setproperties("Title", projecttitle,
"UnitForce", "kN",
"UnitLength", "m",
"UnitTime", "day",
"ModelType", "Axisymmetry")
g.Soilcontour.initializerectangular(0, 0, modelWidth, 1)
matSetList = HBSParam(param)
# nMat = len(matSetList)
soilmat_set_parameters.make_soilmat(g, mat)
# SOIL
firstBorehole = g.borehole(0)
firstBorehole.Head.set(-1)
g.soillayer(1)
g.setsoillayerlevel(firstBorehole, 0, 0)
g.setsoillayerlevel(firstBorehole, 1, -soilLayerThickness)
g.setmaterial(g.Soils[-1], g.Materials[0])
# STRUCTURES
g.gotostructures()
topLineDispl = g.linedispl(topLine)
g.set(topLineDispl.Displacement_x, "Free")
g.set(topLineDispl.Displacement_y, "Fixed")
g.arrayr(topLine, 2, 0, -soilLayerThickness)
g.set(pressureLoad.qx_start, param['supportPressure'])
g.set(pressureLoad.qy_start, 0.)
bcDispl = g.linedispl(bcLine)
g.set(bcDispl.Displacement_y, "Free")
g.set(bcDispl.Displacement_x, "Fixed")
# MESH
g.gotomesh()
g.mesh(0.06)
# STAGED CONSTRUCTION
g.gotostages()
initialPhase = g.Phases[0]
if mat.TypeName == "SoilMat":
excavationPhase = g.phase(g.InitialPhase)
g.set(excavationPhase.MaxCores, 1)
if mat.UserModel == "HB-GSM":
else:
else:
if (mat.User12.value == 0):
else:
g.set(excavationPhase.TimeInterval, param['excavationTime'])
g.set(excavationPhase.Deform.UseDefaultIterationParams, False)
g.set(excavationPhase.Deform.MaxLoadFractionPerStep, 0.1)
g.deactivate(g.Soils[0], excavationPhase)
g.deactivate(g.Deformations, excavationPhase)
g.activate(g.LineDisplacements, excavationPhase)
g.activate(g.LineLoads, excavationPhase)
g.calculate()
soilmat_set_parameters.py
"""
When using the soilmat command (or setproperties command) not all defined
- some parameters are derived parameters, and in the soilmat command these are
not set
The function make_soilmat will try to overcome this by using the soilmat
command and then checks if each of the defined parameter values are set
correctly.
Disclaimer
This Python script is made available as a service to PLAXIS users and can
You use this Python script at your own responsibility and you are solely
responsible for its results and the use thereof. Plaxis does not accept
any responsibility nor liability for the use of this Python script nor
"""
""" Creates or updates a soil material dataset with the defined parameter
values.
"""
# to be expanded
# create material
if material is None:
material = g_i.soilmat(*parameters)
else:
# update properties:
g_i.setproperties(material, *parameters)
soil_attributes = dir(material)
param_name = item[0]
param_value = item[1]
paramname_input = None
if soil_attribute_name.lower() == param_name.lower():
paramname_input = soil_attribute_name
try:
# param_obj.value,
# param_value))
param_obj.set(param_value)
except:
return material
def trial(g_i):
("Colour", 16761441),
("SoilModel", 3),
("DrainageType", 0),
("cref", 0.1),
("phi", 37.5),
("psi", 0),
("gammaUnsat", 18.5),
("gammaSat", 19.5),
("DilatancyCutoff", True),
("einit", 0.49),
("emin", 0.43),
("emax", 0.87),
("Rinter", 0.67),
("OCR", 1),
("POP", 0),
("DefaultValuesAdvanced", True),
("nu", 0.2),
("E50ref", 60000),
("powerm", 0.5),
("K0nc", 0.3912)
def trial2(g_i):
("Colour", 9323),
("SoilModel", 100),
("UserDLLName", "shansep64.dll"),
("DrainageType", 1),
("gammaUnsat", 10.1),
("gammaSat", 10.1),
("Gref", 22.7272727272727),
("User1", 250),
("User2", 0.25),
("User3", 2.5),
("User4", 15),
("User5", 0),
("User6", 0),
("User7", 0.45),
("User8", 0.85),
("User9", 20),
("User10", 3.15),
("User11", 1),
("EoedRef", 250),
("phi", 15),
("K0Primary", 0.741180954897479),
("K0Secondary", 0.741180954897479),
("UDPower", 0),
("UDPRef", 100)]
if __name__ == "__main__":
trial2(g_i)
postProcessFluidity.py
from plxscripting.easy import *
import math
A = 0.
for i in range(len(stressVal)-1):
A = A + abs(rVal[i+1]-rVal[i])*(stressVal[i]+stressVal[i+1])/2
return A
resultsX = []
resultsY = []
for i in range(sample_count):
start[1] + i * step[1])
resultsX.append(position[0])
resultsY.append(float(result_string))
SRR_Res = []
Table_Arel = []
plt.style.use('seaborn-bright')
plt.figure(figsize=(18, 9))
plot_left.set_xlabel("Radius")
plot_left.set_ylabel("Radial Stress")
print(phase.Identification)
SXX = g_o.ResultTypes.Soil.SigxxT
SRR_Res.append(SRR)
Decomp = phase.Identification.value.split()
if Decomp[0] == 'No':
else:
A = area(R, SRR)
plot_right.axis('off')
table = plot_right.table(cellText=Table_Arel,
colLabels=column_labels,
cellLoc='center',
loc='center')
table.scale(1.2, 3.)
if key[0] == 0:
cell.set_fontsize(14)
else:
cell.get_text().set_fontstyle('italic')
cell.set_fontsize(12)
plt.show()
© 2022 Bentley Systems Incorporated. Bentley, the Bentley Logo, and PLAXIS are either registered or unregistered trademarks or service marks of Bentley Systems, Incorporated, or one of
its direct or indirect wholly-owned subsidiaries. Other brands and product names are trademarks of their respective owners.