Earthquake loading

[muhammedfurkanyilmaz]

Hi @Kratos/geomechanics,
I'm trying to make a site response analysis. I exerted the acceleration trace as acceleration at the bottom nodes of the domain, but the acceleration outputs of the corresponding nodes are quite different from the input motion. What could be the reason?

[rfaasse]

Hi @muhammedfurkanyilmaz,
Thank you for reaching out! Could you please add a (simplified) reproducer (i.e. input files), such that we could have a look and try to assist you?
Kind regards,
Richard

[muhammedfurkanyilmaz]

Hi @rfaasse,
Thanks for the quick reply. Here is the file you requested.
MaterialParameters.json
ProjectParameters_stage1.json
1D-EQResponse.mdpa.txt

Kind regards
Furkan

[rfaasse]

Hi @muhammedfurkanyilmaz,

Thank you for this clear reproducer, we'll take a look and let you know!

Kind regards,
Richard

[rfaasse]

Hi @muhammedfurkanyilmaz

I had a look into your input files and was able to reproduce your issue of acceleration values were not as expected at the nodes which had a prescribed acceleration.

The option in the apply_vector_process that should ensure the acceleration is kept at the prescribed value is "is_fixed". We've tried adding that to the process, but still the acceleration was not as expected (the process would look like this in your case):

{
    "python_module": "apply_vector_constraint_table_process",
    "kratos_module": "KratosMultiphysics.GeoMechanicsApplication",
    "process_name":  "ApplyVectorConstraintTableProcess",
    "Parameters":    {
        "model_part_name": "PorousDomain.Line_acceleration",
        "variable_name":   "ACCELERATION",
        "active":          [true,false,false],
        "is_fixed":        [true,false,false],
        "value":           [1.0,0.0,0.0],
        "table":           [1,0,0]
    }
}

We've had a quick look at the code, but couldn't find the problem straightaway, I created an issue to look into this case further and fix it (see #12184). We'll discuss in the team when we can pick this up and will keep you posted.

One last question: With the original 'solver_type' (amgcl), the reproducer model did not complete successfully for me and stopped after one time step. Does it complete all time steps for you? And if so, which operating system are you running the calculation on?

Kind regards,
Richard

[muhammedfurkanyilmaz]

Hi @rfaasse,
Thanks for the feedback. amgcl solver is completing the whole analysis on Windows 10. I've tried to use bigcstab and this also worked.

Kind regards
Furkan

[muhammedfurkanyilmaz]

By the way, I've tried to validate this model under vertical harmonic (Coswt) loading applied at the top by the analytical solution. After reaching the steady state the calculated pore pressures and vertical displacements showed a good agreement with the exact solution. The only difference was the drainage was permitted at the top. I also solved the same problem(vertical loading) on Plaxis2D results was similar.

Kind regards
Furkan

[rfaasse]

@muhammedfurkanyilmaz An update from our side: We are working to fix the issue regarding prescribed accelerations (to follow the progress, see #12184) and will let you know when it is on our master branch.

We are also very interested in the other vertical harmonic loading example. If you could provide us with the input files for that same model, we can put it in our automated test-suite and make sure the results will stay as expected when making changes.

Kind regards,
Richard

[rfaasse]

Hi @muhammedfurkanyilmaz,

We have a PR ready to fix this issue on the master branch. For that we would like to add an adaptation of the test case you added earlier for inspection to the automated test suite in the repository (such that we ensure this keeps working in the future). Do you agree with this? We will add your name as author of the test-model if you agree.

Kind regards,
Richard

[muhammedfurkanyilmaz]

Hi @rfaasse,

Of course, you can use it. I'd be happy with that.

Kind Regards
Furkan.

[rfaasse]

Hi @muhammedfurkanyilmaz,

The fix has been merged, so for the latest version of master, the prescribed accelerations (and velocities) should work correctly. We also added automated test cases (for the documentation, see GitHub).

To get your case running on my machine, I had to make a few adjustments, I attached the .zip here containing the project files and mdpa.

We faced an issue when running this model with the diff order elements (the system matrix seems to be singular, with the following error: RuntimeError: Error: THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT 2). Since you didn't seem to face this issue, you might just be able to run using diff order elements, but just FYI if you happen to see this issue. We'll have a look at it somewhere in the coming time.

The change that is necessary in your model to get the prescribed accelerations to work after our fix, is to add the 'is_fixed' entry to the vector constraint process of the acceleration, such that the process looks like this:

{
    "python_module": "apply_vector_constraint_table_process",
    "kratos_module": "KratosMultiphysics.GeoMechanicsApplication",
    "process_name":  "ApplyVectorConstraintTableProcess",
    "Parameters":    {
        "model_part_name": "PorousDomain.Line_acceleration",
        "variable_name":   "ACCELERATION",
        "active":          [true,false,false],
        "is_fixed":        [true,false,false],
        "value":           [0.0,0.0,0.0],
        "table":           [1,0,0]
    }
}

This ensures the prescribed acceleration will be fixed to that value and no updates will be done for this specific component.

Furthermore, when using the UPw solver, it's always a good idea to add a process to set boundaries/values for the water pressures. Even if the water pressures are not relevant for a certain case, we have seen leaving it undefined leading to issues. Therefore, we added the following process for your model (to just set the pressure to 0 everywhere):

{
        "python_module": "apply_scalar_constraint_table_process",
        "kratos_module": "KratosMultiphysics.GeoMechanicsApplication",
        "process_name":  "ApplyScalarConstraintTableProcess",
        "Parameters":    {
          "model_part_name": "PorousDomain.Soil_two_phase",
          "variable_name":   "WATER_PRESSURE",
          "active":          true,
          "is_fixed":        true,
          "value":           0.0,
          "table":           0
        }
}

Note that in this case the Soil_two_phase model part, also needs to be added to the process model part list:
"processes_sub_model_part_list": ["BottomCorners","Bottom","Side","Line_acceleration","Soil_two_phase"]

When running the attached model with the fixed accelerations, the values at the bottom nodes reflected the prescribed values as found in the table at their respective time-steps.

Feel free to test the newest version of master and see if the fix also works for you!

Kind regards,
Richard

[muhammedfurkanyilmaz]

Hi @rfaasse,
Firstly I'd like to thank you for your great effort. When I fixed the increase factor to 1, the solution did not converge. And my older models also are not fixing the time step after the update.

Kind Regards
Furkan

[rfaasse]

Hi @muhammedfurkanyilmaz,

We had another look and indeed the solver seems to have trouble converging when we calculate every time-step. It took a few tries, but when tweaking the convergence criteria a bit, I got it to converge in the end with the following input files: 1D-EQResponse.zip. FYI, for me it took about 35 minutes to run with a release build of geomechanics.

When plotting the acceleration in the x-direction at one of the bottom nodes, we see the following pattern through time:

Which is visually equal (I didn't do a quantitative check on the entire domain, just a few time steps) to the input:

Our idea of why the solver has difficulties with converging, is that the prescribed acceleration adds a lot of energy to the system, without it being able to dissipate, which might lead to these difficulties in convergence. This is e.g. seen when looking at the x acceleration at a node in the middel of the model (green curve):

It might be an idea to experiment with absorbing boundaries?

In any case, from our side the prescribed acceleration seems to work as expected. If this is also your observation, feel free to close this discussion. Thank you for reporting this problem to us and of course, if there are any more questions, issues, suggestions or insights you encounter, always feel free to reach out!

Kind regards,
Richard

FYI @WPK4FEM

[muhammedfurkanyilmaz]

Hi @rfaasse,

There is one more question related to this issue. The prescribed acceleration function returns as expected in terms of acceleration but not for displacement. The displacements do not match with the second integration of input acceleration as shown below. What could be the reason for that?

[rfaasse]

Hi @muhammedfurkanyilmaz,

Thank you for your thorough testing. We can reproduce your findings and are looking into the issue.

Kind regards,
Richard

[rfaasse]

Hi @muhammedfurkanyilmaz,

We found that the displacements were overwritten by the solver for the bottom nodes with the prescribed accelerations. A way to quickly fix this, is to not only fix the y-displacements, but also the x-displacements at the bottom corner and bottom model parts (see json snippet at the bottom). This is a bit counter-intuitive, because what happens is that the x-displacements are then fixed as far as the solver is concerned, but still updated using the prescribed acceleration in our newmark time integration scheme.

After making this change, the displacement follows the expected curve:

As a bonus it converges a lot quicker, since it only needs ~1 iteration for each time step.

Snippet for the correct processes.

    {
        "python_module": "apply_vector_constraint_table_process",
        "kratos_module": "KratosMultiphysics.GeoMechanicsApplication",
        "process_name":  "ApplyVectorConstraintTableProcess",
        "Parameters":    {
            "model_part_name": "PorousDomain.BottomCorners",
            "variable_name":   "DISPLACEMENT",
            "active":          [true,true,false],
            "is_fixed":        [true,true,false],
            "value":           [0.0,0.0,0.0],
            "table":           [0,0,0]
        }
    },{
        "python_module": "apply_vector_constraint_table_process",
        "kratos_module": "KratosMultiphysics.GeoMechanicsApplication",
        "process_name":  "ApplyVectorConstraintTableProcess",
        "Parameters":    {
            "model_part_name": "PorousDomain.Bottom",
            "variable_name":   "DISPLACEMENT",
            "active":          [true,true,false],
            "is_fixed":        [true,true,false],
            "value":           [0.0,0.0,0.0],
            "table":           [0,0,0]
        }
    }

Hope this helps!

Kind regards,
Richard

[muhammedfurkanyilmaz]

Thanks @rfaasse ,

I appreciate your help.

Best Regards
Furkan