KratosMultiphysics / Kratos

Kratos Multiphysics (A.K.A Kratos) is a framework for building parallel multi-disciplinary simulation software. Modularity, extensibility and HPC are the main objectives. Kratos has BSD license and is written in C++ with extensive Python interface.
https://kratosmultiphysics.github.io/Kratos/
Other
1.01k stars 244 forks source link

[GeoMechanicsApplication] Prescribed acceleration does not result in the expected acceleration values in the respective nodes. #12184

Closed rfaasse closed 5 months ago

rfaasse commented 6 months ago

As a Kratos Geomechanics User, I would like prescribed accelerations to be reflected in the respective nodes, such that I can use this functionality in my analyses.

Acceptance Criteria

Given the user uses a vector constraint process to prescribe the acceleration on certain a certain modelpart. When the user runs a simulation. Then the acceleration at the nodes of these model parts are the same as the prescribed acceleration.

Discussed in https://github.com/KratosMultiphysics/Kratos/discussions/12178

Originally posted by **muhammedfurkanyilmaz** March 13, 2024 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 commented 6 months ago

Had another short look into this bug before we discuss it in planning: It seems to be related to the updating of the second order derivatives in the time integration scheme, where the is_fixed flag is not taken into account for the acceleration (same for the velocity). When hard-coding the following check into the update function (if (rNode.IsFixed(ACCELERATION_X)) continue;, since the reported issue was about prescribed acceleration in the x-direction), the nodal acceleration seems identical to the prescribed acceleration.

for (const auto& r_second_order_vector_variable : this->GetSecondOrderVectorVariables()) {
            if (!rNode.SolutionStepsDataHas(r_second_order_vector_variable.instance)) continue;
            if (rNode.IsFixed(ACCELERATION_X)) continue;

            noalias(rNode.FastGetSolutionStepValue(r_second_order_vector_variable.second_time_derivative, 0)) =
                ((rNode.FastGetSolutionStepValue(r_second_order_vector_variable.instance, 0) -
                  rNode.FastGetSolutionStepValue(r_second_order_vector_variable.instance, 1)) -
                 this->GetDeltaTime() * rNode.FastGetSolutionStepValue(
                                            r_second_order_vector_variable.first_time_derivative, 1) -
                 (0.5 - GetBeta()) * this->GetDeltaTime() * this->GetDeltaTime() *
                     rNode.FastGetSolutionStepValue(r_second_order_vector_variable.second_time_derivative, 1)) /
                (GetBeta() * this->GetDeltaTime() * this->GetDeltaTime());
        }

This is no solution yet (since we need to be able to only fix certain components), but it seems that correctly taking into account the component-wise 'is_fixed' in the time integration scheme will fix the reported issue.

This means however, we need to keep in mind:

FYI @avdg81 @WPK4FEM