RE; Use of Drucker prager material in soil

[ryancwl95]

Hi,

I am modeling a monopile embedded in sandy soil, subjected to displacement-controlled monotonic horizontal loading under drained conditions. At the end of the analysis, I aim to extract the force-displacement curve resulting from this loading and compare it with results from other FEM software for validation.

Currently, I am using the Drucker-Prager model (nDmaterial DruckerPrager) to model the soil, employing SSPbrick elements. For the monopile, I am also using SSPbrick elements to accurately represent its shape and thickness, assuming a modulus of elasticity (E) of 2E12 kPa to treat the pile as fully rigid.

The steps in my analysis include: 1) a gravity stage (purely soil), 2) pile installation, and 3) [next steps, please specify].) loading.

The purpose of this analysis is to load the pile until failure and record the corresponding force required for this displacement. Below is an image after the pile installation and during the loading process, which illustrates the model schematic. I am using a very coarse mesh for the preliminary analysis.

Screenshot 2025-03-24 235821.png (62.04 KiB) Viewed 2341 times

I am entering the weight of the pile and the soil by using the mass density input in the ndmaterials attribute. For element formation, I am specifying b3 (the body force in the z direction) as -9.81 * the density of the soil (ρ). Additionally, I am applying the specified displacement (displacement: 4 times of pile diameter) on a selected master node and measuring the reactionforce. Since the pile has a high modulus of elasticity (E = 2E12 kPa), it behaves rigidly. The density of pile is 7.85Mg/m3 and for the soil is 2Mg/m3. Below is the deformed shape.

Screenshot 2025-03-25 000444.png (69.88 KiB) Viewed 2341 times

I am using an analysis setting of:
puts " Start Loading "
# Analysis settings
#---------------------------------------------
constraints Penalty 1.0e+12 1.0e+12
test NormDispIncr 1e-4 100 2
algorithm KrylovNewton
numberer RCM
system Mumps
integrator Newmark $gamma3 $beta3
analysis VariableTransient
analyze 10000 2.000000e-03 2.000000e-04 2.000000e-03 100
puts "LOADING ANALYSIS COMPLETED.."

After the analysis, I'm receiving an fx (kN) vs ux (m) graph of:

Screenshot 2025-03-25 000930.png (9.77 KiB) Viewed 2341 times

The curve shows some yielding at the beginning (due to the curvature), but afterwards, Fx appears to increase linearly with displacement, indicating elastic behaviour, which seems odd. I'm hoping you could shed some light on this problem.
Thank you in advance! I appreciate your help and look forward to your insights.

I have run several analyses on other FEM software and limit analysis software; failure is expected at a displacement of 10m.

[STKO Team]

The output is not necessarily wrong.
At least from what I can see from the image, it shows a plastic behavior with a quite strong strain hardening.
But without any further information it's impossible to tell if there is a problem.

[ryancwl95]

Hi,
Thank you for your response; this information is helpful. I am modelling the soil based on one that exhibits no softening or strengthening/hardening behaviour, utilizing an associative flow rule for plastic volumetric change.

I would like to double-check my inputs for such behaviour, specifically for $rhobar to $theta.

"nDMaterial DruckerPrager $matTag $k $G $sigmaY $rho $rhoBar $Kinf $Ko $delta1 $delta2 $H $theta $density <$atmPressure>"

My input is as follows:
$rhoBar = $rho
$Kinf, $Ko, $delta1, $delta2, $H = 0
$theta = 1

Could I clarify the appropriate values for my intended behaviour?

[STKO Team]

$rhoBar = $rho

Yes, this means "associated" plasticity

$Kinf, $Ko, $delta1, $delta2, $H = 0

Yes, this means "perfect plasticity" not strain hardening or softening.

However, if rhoBar >= 0, your model is dilatant and will try to expand upon shear deformation. If your model is also constrained not to expand (due to the layout of your boundary conditions), the material will be subject to increasing confining stresses, thus showing an apparent strain hardening

[ryancwl95]

Hi,
Thank you for this reply.
I have tried applying this Drucker Prager ndmaterial to a simpler problem. This problem is a loading of circular footing in soil (Drucker Prager ndmaterial).
The values of $rhoBar, $Kinf, $Ko, $delta1, $delta2, $H and $tetha used in the DP soil model is the same as mentioned above. For this example, I'm still receiving similar force-displacement curves.

Brief description of this example:

Screenshot 2025-04-14 151805.png (227.99 KiB) Viewed 2146 times

- A quarter of the circular footing is modelled.
- All continuum elements are modelled using SSPbrick elements
- Fixities are applied at the base of the model, x = 0 and y =0 (plane) using:
fixX 0.000000 1 0 0 -tol 1.000000e-10
fixY 0.000000 0 1 0 -tol 1.000000e-10
fixZ -2.000000 1 1 1 -tol 1.000000e-10
- Fixities applied in the x and y direction are applied at the circumferential nodes bounding the model
- A weightless footing is modelled using nDMaterial ElasticIsotropic, with a high value of E (kPa) to simulate fully rigid behaviour.

Screenshot 2025-04-14 151942.png (77.71 KiB) Viewed 2146 times

- Nodes belonging to the footing elements are fixed in the x and y direction so that they can only move vertically.
- Nodes at the base of the footing are tied to the soil nodes below using equalDof in 1 2 3 dof (x, y, z direction).
- a point displacement (uz = -0.02m)is applied at the top centre of the circular footing (x=0, y=0). The applied displacement is applied linearly in this format:

pattern Plain $patternTag {Series -time {0 50 51} -values {0 1 1.020000e+00} -factor 1} {
sp 4000001 3 -2.000000e-02
}

For the analysis stages:
1 - Gravity stage (soil only), at this stage I applied InitialStateAnalysis On
2 - Footing installation (InitialStateAnalysis is OFF before the start of this stage)
3 - Loading of the footing

Thereafter, I used the STKO post-processing to retrieve the reaction force (without inertia) and displacement in the z-direction at the 3rd stage.

And this is the graphplot that I've gotten.

Screenshot 2025-04-14 151423.png (12.6 KiB) Viewed 2146 times

According to what you have mentioned, this apparent strain hardening is reflected by confining stresses, which could be caused by boundary and fixity conditions. I would like to check if the bounding conditions I have applied are appropriate.

Thank you for your help!

[ryancwl95]

Hi,

I was wondering if there's any advice you can give on this situation, any will be helpful. Thank you in advance!

Sincerely,
Ryan

[STKO Team]

I think your model is too small in the vertical direction to allow the slip plane to happen

[ryancwl95]

Hi

Thank you for the advice. I have attempted to expand the model's depth to about 20 meters, increasing it from the previous depth of 2 meters. However, I am still unable to observe any yielding behaviour.

Below is the analysis geometry showing the model at a vertical depth of 20 meters, along with the results of the Fz vs. Uz plot at different depths: 2 meters (initially shown), 4 meters, 10 meters, and 20 meters.

Screenshot 2025-05-11 222956.png (140.21 KiB) Viewed 1974 times

Screenshot 2025-05-11 223115.png (33.11 KiB) Viewed 1974 times

Could you advise on the next steps to demonstrate the yielding behaviour of soil using the Drucker-Prager soil model for a rigid footing subjected to monotonically vertical loading?

[STKO Team]

I would suggest you to take a look at the literature, and find a reference benchmark to reproduce.
The failure type will depend on other material properties (like internal friction angle and dilatancy angle).

[ryancwl95]

Hi,

Thank you for the suggestions; I will review the existing literature.

The current literature does not fully address the plastic behavior of a monotonically loaded monopile problem. The Fx-Ux curve presented in the literature typically only shows displacement up to 0.1D, which is generally sufficient for design purposes.

However, I am attempting to simulate the complete plastic failure of this monopile problem under displacement-controlled loading conditions. This preliminary analysis aims to assess the program's capability to model complete plastic failure before examining the transition from elastic to fully plastic behavior.

I have modified the material properties and set both the friction angle and the dilatancy angle to 20 degrees (on the low side) and 20 degrees (for associated flow), respectively.

Am I overlooking something necessary to model complete plastic failure, or is the program incapable of representing the plateau curve in the Fx vs. Ux plot? If so, could you explain the reason for this limitation?

Sincerely,
Ryan Chia