Removing solid elements

[[email protected]]

Dear STKO team,

I am trying to create a model for conducting lateral dynamic pile tests. For this I need to do the 1) consolidation analysis of the soil first then 2) remove the soil from the volume of the pile, place in the pile then 3)conduct settlement analysis for the pile under self-load and applied vertical top-mass and finally 4) apply lateral dynamic loading and run a transient analysis.

How can part 2) be done in such a way that the removal of soil and placement/connection of the pile into the soil will not cause spurious stresses during my dynamic analysis in part 4)?
Do you perhaps have some e-learning material for a similar geotechnical problem?

Thank you for your time in advance!

Regards,
Bence

[STKO Team]

Here is a modified version of your model that does all the steps you required.
For the sake of simplicity, I converted all materials to linear elastic ones, so that now the only nonlinearity is contact (and it is quite a hard nonlinearity to master...).
Once you get confident with contact, go ahead with material nonlinearity.

Since it will take time to explain in detail all I did in this model, first have a look at it. Try to understand what I did, and then, if you have doubts, post them here one by one, so that it will explain them.

PSI.zip

(314.27 KiB) Downloaded 336 times

To simplify the procedure of using contact with UP element, I made a small change in the BeamContact and BeamEndContact 3d elements.
Unzip the content of this folder and copy the python files in:
C:\Program Files\STKO\external_solvers\opensees\element_properties\contact_elements
(Replace the old ones)

contact_elements.zip

(6.87 KiB) Downloaded 295 times

And the same for this:

fix.zip

(4.09 KiB) Downloaded 300 times

Unzip it, and copy it here:
C:\Program Files\STKO\external_solvers\opensees\conditions\Constraints\sp


All these changes will be available in the next (2.0.1) version of STKO

[[email protected]]

Dear STKO team,

Thank you a lot for the awesome help and the element removal script!

First, a general request. Can you please describe me the BeamContact3D penalty element? How to choose the penalty? Why we need this penalty and what does it do? How can we modify the gap and force tolerance parameters for this version of the BeamContact3D?

I have several practical questions about your process and some new issues I've encountered:


1) Why do you use the transient analysis for applying self-weight analyses? (I was unable to successfully run the model when I switched the conditions to static analysis)

2)Will applying the soil self-weight within the element instantaneously will create problems when using non-linear soil material? Maybe it's practical to apply soil weight as a body force in geotechnical OpenSees models.

3)In the attached model you have fixed the Pressure for the entire model throughout all the analyses. This makes using the UP element pointless. However, when I change this fixity to a pressure-free surface (fix P at 0 on top surface) no porewater pressure is generated in the model during any of the analyses. When I fix it to be a certain pressure (simulating sea-level) using sp the fixed porewater pressure still changes during the analysis. What might cause this?
Also, the BeamContact3D element only connects 3DOF, is this the part that you have solved with the updated python files, or do I need to implement phantom nodes with eDOF to connect the UP elements with the contact elements?

4) It seems that the frictional material assigned to the beamcontact3D element (lateral) does not activate skin friction during the vertical pile loading stage. Apparently the element slides freely along the vertical soil face and only affects the soil around the pile toe. Can you confirm this?

5) During the analysis Opensees prompts the following warnings for the contact3D elements:
BeamEndContact3Dp::setResponse(const char **argv, int argc, OPS_Stream &eleInfo): material unknown request
BeamEndContact3Dp::setResponse(const char **argv, int argc, OPS_Stream &eleInfo): localForce unknown request
BeamEndContact3Dp::setResponse(const char **argv, int argc, OPS_Stream &eleInfo): deformation unknown request
BeamEndContact3Dp::setResponse(const char **argv, int argc, OPS_Stream &eleInfo): localAxes unknown request
BeamContact3Dp::setResponse(const char **argv, int argc, OPS_Stream &eleInfo): section unknown request

Are these related to the recorder (they may not be appropriate recorder queries) or is this an issue that can cause problems with the analysis? E.g. the local Z axes assigned to the contact element all point in one direction and do not align with the interaction lines.

6) I've changed the dynamic analyses according to my problem. It's a dynamic pile load test. (dynamic loading at pile head), however I had convergence issues after the first few seconds of the analysis. Which parameters are wrong in my attached model?

7) Finally, a suggestive question. If I use separate volumes for the soil domain and the soil to be removed then attach them with eDOF for the geostatic analysis, then remove all the elements and eDOFs can we evade the trouble of selecting and keeping the interface nodes that are shared in the merged model (model you provided)?

I am attaching my slightly modified file.

Thank you very much for your help in advance!

Bence

[STKO Team]

Can you please describe me the BeamContact3D penalty element? How to choose the penalty? Why we need this penalty and what does it do? How can we modify the gap and force tolerance parameters for this version of the BeamContact3D?

The BeamContact3D uses the Lagrange Multipliers method, while the BeamContact3Dp uses the Penalty method to enforce Contact condition.
The Lagrange Multiplier can be tricky in dynamic analysis, that's why I tried with the Penalty.

1) Why do you use the transient analysis for applying self-weight analyses? (I was unable to successfully run the model when I switched the conditions to static analysis)

Because the UP elements work only in dynamic analyses, so if you want to use those in a quasi-static problem, you have to run a dynamic analysis with very large time steps so that iniertia effect are zero.

2)Will applying the soil self-weight within the element instantaneously will create problems when using non-linear soil material? Maybe it's practical to apply soil weight as a body force in geotechnical OpenSees models.

Of course, in that case, use the VolumeForce tool in STKO. It is input a force per unit volume in STKO, but then STKO converts it to equivelent nodal loads that can be used in load patterns.

3)In the attached model you have fixed the Pressure for the entire model throughout all the analyses. This makes using the UP element pointless. However, when I change this fixity to a pressure-free surface (fix P at 0 on top surface) no porewater pressure is generated in the model during any of the analyses. When I fix it to be a certain pressure (simulating sea-level) using sp the fixed porewater pressure still changes during the analysis. What might cause this?

I fixed all P DOFs to make a comparison with the standard U element.

Also, the BeamContact3D element only connects 3DOF, is this the part that you have solved with the updated python files, or do I need to implement phantom nodes with eDOF to connect the UP elements with the contact elements?

Yes, now if the contact finds a node with num dofs != 3, it will create a new 3DOF node at the same position, linked with EDOF at the common DOFs.

4) It seems that the frictional material assigned to the beamcontact3D element (lateral) does not activate skin friction during the vertical pile loading stage. Apparently the element slides freely along the vertical soil face and only affects the soil around the pile toe. Can you confirm this?

Yes because you gave a small G (tangential stiffness in the contact material) and also because there is not large enough confining stress to create the stick condition. Try to increase G to 1.0e10, and use large (1.0e200) C and t. I know it is not real, but will enforce the stick condition even if there is no compressive stress (this is just to show how it works)

Are these related to the recorder (they may not be appropriate recorder queries) or is this an issue that can cause problems with the analysis? E.g. the local Z axes assigned to the contact element all point in one direction and do not align with the interaction lines.

They are just warnings, don't bother

6) I've changed the dynamic analyses according to my problem. It's a dynamic pile load test. (dynamic loading at pile head), however I had convergence issues after the first few seconds of the analysis. Which parameters are wrong in my attached model?

You are using very large penalty values for constraints. Also use the Newton with linesearch, it is better for contact problems.

7) Finally, a suggestive question. If I use separate volumes for the soil domain and the soil to be removed then attach them with eDOF for the geostatic analysis, then remove all the elements and eDOFs can we evade the trouble of selecting and keeping the interface nodes that are shared in the merged model (model you provided)?

It will work, it's just more complicated because in the TCL script for the removal, you have to remove MP constraints as well.



Here is a working file. During you final analysis it will have convergence problem. Probably it is due to forces that are very high. Try reducing the time step, or using smooth timeSeries (the pulse one creates shocks)

PSI_try_pile_load.zip

(281.81 KiB) Downloaded 315 times

[[email protected]]

Thank you for your prompt response!

You are using very large penalty values for constraints. Also use the Newton with linesearch, it is better for contact problems.

How did you define the 1e+12 penalty? In your seminars I've heard that for geotech problems it should be soil Young's modulus * element volume and add 1e+8 order. In my meshed model, the largest element has a volume of 40m3, which considering E = 2e+6 gives me order of magnitude e+7 + e+8 = e+15.

The smaller penalty with sine wave and 1/1000 of the load still doesn't converge. I'll keep trying to play around with the parameters. Thank you for your time!

[marafini.f]

Dear bence,

I have not been involved in reviewing your model. But regarding penalty parameters, in theory they should be infinite, as it is no possible in numerical problems, in practice they cannot be too large nor too small. A common rule to defined them is:

• Evaluate (approximately) the largest entry on the global (assembled) stiffness matrix diagonal... not the exact value, just an approximation is fine. Which means the largest stiffness you find in your model.

• Get the order of magnitue: O = ceil(log10(Kmax)).

• Penalty = 10^(O+8). Where 8 is half the significant number of digits for a double-precision floating point number

Remember that the penalty method does not work well with the norm unbalance test.

If you want to try with other setups, in terms of constraints handler, algorithms or integrator, I can suggest you to try using the search box in the forum. You will find explanations from problems faced by other users , and maybe some advice on how to improve your model.

Hope this helps.

Enjoy your modeling.

[STKO Team]

How did you define the 1e+12 penalty? In your seminars I've heard that for geotech problems it should be soil Young's modulus * element volume and add 1e+8 order. In my meshed model, the largest element has a volume of 40m3, which considering E = 2e+6 gives me order of magnitude e+7 + e+8 = e+15.

This is true, but only approximately. Let's say that adding 8 to your maximum order of magnitude is the upper bound.
In fact, in your previous model, you were using 1.0e17 which was higher and it was not converging at all.
Sometimes, in models with a lot of constraints (in your model for example you have a lot of edof, and contact), you have a lot of these penalty stiffnesses, and in those cases, it's better to lower the upper bound.

Typically what I do is to first estimate the upper bound, then I try it and scale it down one order of magnitude at a time, to find the best one.
As soon as possible we will implement a tool in OpenSees that computes the best Penalty parameter for you.

[[email protected]]

Thank you very much for both of your replies.

I will continue browsing the forum to try to find solutions to my convergence issue.

Is there any reading material/online resource you could suggest that describes the pros/cons and applicability of the system solvers, integrators and iterating algorithms which can help me understand how to choose and parameterize them for non-linear dynamic problems? Perhaps a user manual of a mature FEM software?

[STKO Team]

Is there any reading material/online resource you could suggest that describes the pros/cons and applicability of the system solvers, integrators and iterating algorithms which can help me understand how to choose and parameterize them for non-linear dynamic problems? Perhaps a user manual of a mature FEM software?

You won't find online a unique user manual that will explain to you all the possible options you have in OpenSees. Especially regarding the transient integrators.

Regarding the algorithms this is a brief guideline:

I generally use the standard Newton-Raphson, which is the typical one explained in every FEM book. It guarantees quadratic convergence, but only when your nonlinear function is quite smooth so that its derivatives are quite stable.
Newton is great at predicting the direction of the solution (using the tangent), however it may fail in predicting the magnitude of the solution (when for example the tangent approaches zero, or when the tangent changes rapidly in non-smooth problems like Contact).

For contact problems, I prefer Newton+LineSearch. LineSearch is a technique that tries to perform a better prediction about the magnitude of the solution.

Instead for problems involving material instabilities (such as softening) I prefere the krylov-newton, or ModifiedNewton with initial tangent.

Regarding the integrators, I suggest you look at the OpenSees wiki and refer to the references given therein for each integrator type.

[[email protected]]

That's very helpful advice, thanks a lot for the kind help!