Hello,

I've posted, at the end of message, the script I use for the stress initialization of a deep tunnel.
I've noticed that the initial y stress value (isyy1 parameter) is reset after solving to equilibrium.
Is ther something wrong with the script?
Thank You




// create new/fresh model
// the model is suited for deep tunnel wiyh only 1 layer
newmodel()
// set unit to m, kN, kPa, ton/m3
set("steplimit",100000,"unit","stress-kpa")

//model dimension
var lmax = 80 //model width
var dmax = 50 //model height
var water_depth = 600
var tunnel_crown_depth_from_top_of_model = 25
var initial_y_stress = 7680 //vertical stress at the topo of model in kPa

// MODEL MATERIALS
//layer 1 characterization
var rho1 = 2.7 //mass density, ton/m3
var nu1 = 0.3 //poisson's ratio
var shear_modulus1 = 15385000 // elastic shear modulus, kPa
var bulk_modulus1  = 33333000 //elastic bulk modulus, kPa
var ucs = 45000 //Hoek-Brown parameter, σci, kPa
var mb = 0.337 //Hoek-Brown parameter, mb
var s = 0.0002 //Hoek-Brown parameter, s
var a = 0.506 //Hoek-Brown parameter, a
var t = 3000  //current value of tensile strength, σt, kPa
var s3cv = 500 //Hoek-Brown parameter, σcv
//masonry characterization
var rho_mas = 1.835 //ton/m3
var shear_modulus_mas = 2000000 //kPa
var bulk_modulus_mas  = 2666700 //kPa
//concrete characterization
var rho_c = 2.45 //ton/m3
var shear_modulus_c = 15000000 //kPa
var bulk_modulus_c  = 20000000 //kPa
//ballast characterization
var rho_b = 2.035 //ton/m3
var shear_modulus_b = 65400 //kPa
var bulk_modulus_b = 70800 //kPa
var c_b = 0 //cohesion kPa
var phi_b = 38 //friction angle ?
var t_b = 0.0 //tensile strenght kPa

// geotechnical model
rect("startPoint",-lmax/2,-dmax,"endPoint",lmax/2,0)

// generate geometry
arc("startPoint",-4.26381e-17,-25,"midPoint",-0.983496,-25.1956,"endPoint",-1.81726,-25.7527,"numSeg",20)
arc("startPoint",1.81726,-25.7527,"midPoint",0.983496,-25.1956,"endPoint",-4.26381e-17,-25,"numSeg",20)
arc("startPoint",-5.89989e-18,-24.4,"midPoint",-0.820456,-24.508,"endPoint",-1.585,-24.8247,"numSeg",20)
arc("startPoint",1.585,-24.8247,"midPoint",0.820456,-24.508,"endPoint",-5.89989e-18,-24.4,"numSeg",20)
arc("startPoint",-1.81726,-25.7527,"midPoint",-2.46234,-26.6655,"endPoint",-2.76573,-27.7412,"numSeg",20)
arc("startPoint",2.76573,-27.7412,"midPoint",2.46234,-26.6655,"endPoint",1.81726,-25.7527,"numSeg",20)
line("startPoint",7.1e-15,-31.5,"endPoint",-2.11573,-31.5)
line("startPoint",-6.9e-15,-31.5,"endPoint",2.11573,-31.5)
line("startPoint",-2.11573,-31.5,"endPoint",-3.40904,-31.5)
line("startPoint",2.11573,-31.5,"endPoint",3.40904,-31.5)
arc("startPoint",-2.76573,-27.7412,"midPoint",-2.7017,-29.6658,"endPoint",-2.11573,-31.5,"numSeg",20)
arc("startPoint",2.11573,-31.5,"midPoint",2.7017,-29.6658,"endPoint",2.76573,-27.7412,"numSeg",20)
arc("startPoint",-1.585,-24.8247,"midPoint",-1.96567,-25.0791,"endPoint",-2.30987,-25.3809,"numSeg",20)
arc("startPoint",2.30987,-25.3809,"midPoint",1.96567,-25.0791,"endPoint",1.585,-24.8247,"numSeg",20)
arc("startPoint",-2.30987,-25.3809,"midPoint",-3.4728,-27.2063,"endPoint",-3.75531,-29.3522,"numSeg",20)
arc("startPoint",3.75531,-29.3522,"midPoint",3.4728,-27.2063,"endPoint",2.30987,-25.3809,"numSeg",20)
arc("startPoint",-3.75531,-29.3522,"midPoint",-3.62117,-30.4324,"endPoint",-3.40904,-31.5,"numSeg",20)
arc("startPoint",3.40904,-31.5,"midPoint",3.62117,-30.4324,"endPoint",3.75531,-29.3522,"numSeg",20)

// discretize/mesh
triangle("elemtype","T3")
discretize("maxedge",1.5)
triangle("maxedge",1.5)
// create new materials
material("create","Hoek-Brown","matid",1,"matname","Calcschists","density",rho1 ,"shear",shear_modulus1,"bulk",bulk_modulus1,"sigci",ucs,"mb",mb,"s",s,"a",a,"s3cv",s3cv)
material("create","IsoElastic","matid",99,"matname","Masonry","density",rho_mas,"shear",shear_modulus_mas,"bulk",bulk_modulus_mas)
material("create","IsoElastic","matid",98,"matname","Concrete","density",rho_c,"shear",shear_modulus_c,"bulk",bulk_modulus_c)
material("create","Mohr-Coulomb","matid",97,"matname","Ballast","density",rho_b,"shear",shear_modulus_b,"bulk",bulk_modulus_b,"coh",c_b,"fric",phi_b,"dil",0,"tens",t_b)

// assign materials
material("assign","matid",1,"region",0.0,- 0.1)
material("assign","matid",1,"region",0.0,-tunnel_crown_depth_from_top_of_model + 0.1)
material("assign","matid",1,"region",0.0,-tunnel_crown_depth_from_top_of_model - 0.1)
watertable("add","dens",1,"elev",-water_depth)

// assign boundary condition
applybc("xfix","xlim",-lmax/2-0.1,-lmax/2+0.1,"ylim",-dmax-0.1,0.1)
applybc("xfix","xlim",lmax/2-0.1,lmax/2+0.1,"ylim",-dmax-0.1,0.1)
applybc("xyfix","xlim",-lmax/2-0.1,lmax/2+0.1,"ylim",-dmax-0.1,-dmax+0.1)

// enable nodal mixed discretization (NMD) to increase accuracy for T3 elements
set("useNMD","on")

// initialize stress condition
//layer 1
var k01= nu1/(1-nu1)
var dsyy1 = rho1*9.81
//var isyy1 = -dsyy1*top_of_model_depth
var isyy1 = -initial_y_stress
initial("syy",    isyy1,"yvar",dsyy1,"xlim",-lmax/2,lmax/2,"ylim",-dmax-0.1,0)
initial("sxx",k01*isyy1,"yvar",k01*dsyy1,"xlim",-lmax/2,lmax/2,"ylim",-dmax-0.1,0)
initial("szz",k01*isyy1,"yvar",k01*dsyy1,"xlim",-lmax/2,lmax/2,"ylim",-dmax-0.1,0)

// apply gravity force
set("gravity",0,9.81)

//initialize model
//solve()
initial("xydisp",0)
savemodel("filename","st0_init.adf")

//plot principal stresses
plot("contour","sig1")

Hi,

there are couple of issues in your script:

1- I don't see any reason for using stress gradient in x-direction, when you try to assign initial stress. I think you should remove the "xlim",-lmax/2,lmax/2" from the command and only make sure stress varies in y-direction which make more sense.

2- more importantly based on your assumption the stress at the surfaces of the model is -7.6E+3. there is no problem with specifying such stress at the surface if you want to remove some overburden from the top but you have to make sure correct boundary condition is specified to prevent upward movement. In your model, ground surface is free to move upward and with high compression stress in the surface model cycles until stress goes to zero due to expansion.

So I don't see any issue in the program itself so I'll move this topic into "ADONIS Users" category.

-Roozbeh

Hello Roozbeh,

there is no stress gradient in x-direction and the parameter "xlim" is superfluous.
I do not understand why "yvar" for sxx and szz should not vary with k0 (see also tutorial 2).
Ok for the point nr. 2 (copy-paste error without output check).

You are right, yvar needs to be factored and your script was correct. Modified my response.