Problem of convergence with PDMY02 with high friction angle

[francisconazar]

Hi all,

I'm currently experimenting with PDMY's (PressureDependMultiYield) materials in a single stdBrick element, subjecting it to a compressive displacement controlled load with a confining pressure, just to explore and to maybe calibrate the model with some triaxial tests.

I came with a curious problem; if I use a friction angle of 46 degrees or less, everything runs smoothly and converges, and the results make sense. But, if I change the friction angle to 47 or greater -with the same other parameters intact-, the model does not converge with any algorithm (nor with other tests).

I don't believe this is a problem of the model, but I don't know if it is a problem of the material (that maybe it was programmed with low friction angles in mind) or the analysis's parameters I'm using (I have tried the majority of them). Examining the info that the test throws, I notice that at some point in the nonlinear phase the CTestRelativeEnergyIncr ratio starts to diverge and then it becomes 'nan'.

I find it hard to believe that a very simple model with one element has convergence problems with the nonlinear phase, yet I can't find any other explanation.

Any help is welcomed. I attach the code below.

Regards.




wipe

#############################################################
# BUILD MODEL

#create the ModelBuilder
model basic -ndm 3 -ndf 3

node 1 0 0 0
node 2 1 0 0
node 3 1 1 0
node 4 0 1 0
node 5 0 0 1
node 6 1 0 1
node 7 1 1 1
node 8 0 1 1

fix 1 1 1 1
fix 2 0 1 1
fix 3 0 0 1
fix 4 1 0 1
fix 5 1 1 0
fix 6 0 1 0
fix 7 0 0 0
fix 8 1 0 0

equalDOF 5 6 3
equalDOF 5 7 3
equalDOF 5 8 3


# define material and properties

set rho 2100; # kg/m3
#set rho 0
set E 220000000; # Pa
set poisson 0.1;
set refShearModul [expr $E/(2*(1+$poisson))];
set refBulkModul [expr $E/(3*(1-2*$poisson))];
set pres -166000
# set mpres 10000;
# set refShearModul 130; # MPa
# set refBulkModul 200 # MPa
set frictionAng 46 ;# grados
set peakShearStra 0.1 ;# grados
set refPress 80000;# Pa
set pressDependCoe 1.0 ;#
set PTAng 27 ;#
#set contrac 0 ;#
#set dilat1 50 ;#
set liquefac1 0 ;#
set liquefac2 0 ;#
set liquefac3 0 ;#
#$contrac1 $contrac3 $dilat1 $dilat3
set contrac1 1
set contrac2 1
set contrac3 1
set dilat1 0
set dilat2 0 ;#
set dilat3 0
set e 0.6
set cs1 0.0
set cs2 0.0
set cs3 0.0
set pa 101000
set ce 40000; #Pa

#nDMaterial PressureDependMultiYield02 1 2 $rho $refShearModul $refBulkModul $frictionAng $peakShearStra $refPress $pressDependCoe $PTAng $contrac $dilat1 $dilat2 $liquefac1 $liquefac2 $liquefac3 <$noYieldSurf=20 $contrac2=5. $dilat2=3. $liquefac1=1. $liquefac2=0. $e=0.6 $cs1=0.9 $cs2=0.02 $cs3=0.7 $pa=101 <$c=0.1>>
nDMaterial PressureDependMultiYield02 1 3 $rho $refShearModul $refBulkModul $frictionAng $peakShearStra $refPress $pressDependCoe $PTAng $contrac1 $contrac3 $dilat1 $dilat3 20 $contrac2 $dilat2 $liquefac1 $liquefac2 $e $cs1 $cs2 $cs3 $pa $ce; # 0 1.
#nDMaterial PressureDependMultiYield02 100 2 $rho $refShearModul $refBulkModul $frictionAng $peakShearStra $refPress $pressDependCoe $PTAng $contrac1 $contrac3 $dilat1 $dilat3 20 5. 3. 1. 0. 0.6 0.9 0.02 0.7 101000 20000 0 1.

# ele# thick maTag bulk mDensity perm1 perm2 gravity
#element Quad 1 1 2 3 4 5 6 7 8 9 1.0 1 $bulk $fmass $hperm $vperm 0 0

element stdBrick 1 1 2 3 4 5 6 7 8 1 $pres $pres $pres

recorder Node -file NodosDisp.out -time -node 1 2 3 4 5 6 7 8 -dof 1 2 disp;
recorder Node -file NodosReac.out -time -node 1 2 3 4 5 6 7 8 -dof 1 2 reaction;

recorder Element -file Elementostress1.out -ele 1 -time material 1 stress
recorder Element -file Elementostrain1.out -ele 1 -time material 1 strain


#set material to elastic for gravity loading
updateMaterialStage -material 1 -stage 0

#constraints Plain ; # Penalty 1.0e18 1.0e18 ;#
constraints Plain
test NormDispIncr 1.e-5 300 0
#algorithm Newton
#algorithm BFGS
algorithm ModifiedNewton
#algorithm KrylovNewton
numberer Plain
system FullGeneral
#system ProfileSPD
integrator LoadControl 1;
#integrator DisplacementControl 5 3 -0.0001; # displacement control algorithm seking constant increment of 0.1 at node 1 at 2'nd dof.
analysis Static

analyze 1


updateMaterialStage -material 1 -stage 1

#constraints Plain ; # Penalty 1.0e18 1.0e18 ;#
constraints Plain
test NormDispIncr 1.e-5 300 0
#algorithm Newton
#algorithm BFGS
algorithm ModifiedNewton
#algorithm KrylovNewton
numberer Plain
system FullGeneral
#system ProfileSPD
integrator LoadControl 1;
#integrator DisplacementControl 5 3 -0.0001; # displacement control algorithm seking constant increment of 0.1 at node 1 at 2'nd dof.
analysis Static

analyze 1

pattern Plain 1 Linear {
load 5 0 0 -1 ;
# load 2 -1000 0;
# load 3 -1000 0;
#
}

#constraints Plain ; # Penalty 1.0e18 1.0e18 ;#
constraints Plain
test NormDispIncr 1.e-5 300 1
test RelativeEnergyIncr 1.e-5 10000 1
#test NormUnbalance 1.e-5 300 1
#algorithm Newton
#algorithm BFGS
#algorithm ModifiedNewton
algorithm KrylovNewton
numberer Plain
#system FullGeneral
system ProfileSPD
integrator DisplacementControl 5 3 -0.0001; # displacement control algorithm seking constant increment of 0.0002 at node 5 at 3'nd dof.
analysis Static

analyze 100

wipe; #flush ouput stream

[skamalzare]

Hi Francisco,

Friction angle plays a very important role in PDMY and PDMY02 materials. It defines the size of yield surface. The higher the phi, the larger the yield surface.

Another, very important factor in these models is the angle of Phase Transformation Line - (PTL - your model: PTAng). This line defines the boundary between compressive (i.e compression in drained tests, and +ve PWP in undrained tests) and dilative (i.e dilation in drained tests, and -ve PWP in undrained tests) volumetric behaviors. Soil state below PTL shows compressive behavior, and above it shows dilative behavior. Therefore, if you shear a soil with initial isotropic condition, it compresses till reaching PTL, thereafter it dilates to hit the yield surface.

Now, if you define phi very large, and PTAng considerably smaller. You will increase the distance between PTL and yield surface. Therefore, you will have lots of dilation. This is the reason why you are not converging with phi>=47. Because your PTAng is 27, and simply 47-27=20 is too large of a distance for your model. If you increase your PTAng, you will immediately converge.

I think, what probably is happening in your model is that you never reach he displacement control because you have an extremely dilative soil!

I hope this could help you,
Bests,
Soheil

[francisconazar]

Dear Soheil,

Thank you for your response! Now a lot of things are more clear.

Regards.

[francisconazar]

Dear Soheil,

I'm trying to calibrate the model to this experimental data and model (École centrale Paris multimechanism model)
https://i.imgsafe.org/5b06f06.png

This is so far my model
http://hostcode.sourceforge.net/view/5305

and the results are this:
https://i.imgsafe.org/791a74c.png

The matlab code for generating the figure is here if it is of help:
http://hostcode.sourceforge.net/view/5306

I have 2 questions. In your experience:
1. Is it possible to obtain a better match in the volumetric behaviour and the sigma-epsilon law (particularly the first)? I have played around the contrac and dilat parameters and they behave very strange, e.g., sometimes a slight change in contrac1 makes the model diverge. The pressDependCoe also at some point over 2.0 does not converge.
2. Is it possible to obtain softening in the sigma-epsilon law? I read the model is hyperbolic but maybe the dilat and contrac parameters can be chosen to provide some softening (when I played around with them they sometimes they showed softening but I don't know if it was the model or a convergence problem)

Regards.

[sinapourabbasi]

Hi all
Can anyone explain to me what the $Ce coefficient is in the code below?


nDMaterial PressureDependMultiYield02 1 3 $rho $refShearModul $refBulkModul $frictionAng $peakShearStra $refPress $pressDependCoe $PTAng $contrac1 $contrac3 $dilat1 $dilat3 20 $contrac2 $dilat2 $liquefac1 $liquefac2 $e $cs1 $cs2 $cs3 $pa $ce; # 0 1.