[StructuralMechanicsApplication] Some comments about the calculation of the DAMPING MATRIX

[AlejandroCornejo]

Dear @KratosMultiphysics/structural-mechanics team,

After some calculations and time to study some dynamic cases I've come up with the following conclusions:

1. The computation of the Damping matrix D using the Rayleigh formula depends on the M matrix and the stiffness matrix K. This is OK but there are some concerns regarding WHICH stiffness matrix should be used...

2. If we use the TANGENT stiffness matrix of the current iteration (this is what we do now) the solution becomes a bit unstable and it iterates a lot in constitutive non-linear problems (since K is different at each iteration).

3. After discussing this with some professors of the school I've concluded that the best option is to use the original stiffness matrix Ko or the secant/tangent stiffness from the previous time step (this means that the K used for computing the D is constant along the time step)

Knowing this, maybe could be interesting to implement a new method like ComputeLHSforDamping that flexibilizes this operation. Thank you, Alejandro

[philbucher]

Hi @AlejandroCornejo Thanks for the explanation! Do you have some link to literature as well? I only have one comment how to test/facilitate this: I would not add another (and specialized) method to the Elements/Conditions, since they are already quite heavy/have a large interface. Instead such an option could be nicely integrated into the schemes I think.

[RiccardoRossi]

Or through the Calculate function of the element.

[philbucher]

Going trough Calculate requires Element specific modifications (fine I think same is done for the damping matrix) But also then the scheme would need to be modified

[AlejandroCornejo]

Thank you @philbucher and @RiccardoRossi , If you want I can search some literature but i'm afraid that this is a "programatic" issue and it is going to be difficult to find something. I'll try to find something anyway. I'm open to hear options to enhance this!

[loumalouomega]

OK from my side, I have questions how to do it in a clean way, because you can add a calculate which uses the initial configuration, but in case the CL considers something from the geometry (like volume) it can be a little bit more difficult to control.

El vie., 10 may. 2019 9:09, Philipp Bucher notifications@github.com escribió:

Hi @AlejandroCornejo https://github.com/AlejandroCornejo Thanks for the explanation! Do you have some link to literature as well? I only have one comment how to test/facilitate this: I would not add another (and specialized) method to the Elements/Conditions, since they are already quite heavy/have a large interface. Instead such an option could be nicely integrated into the schemes I think.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/KratosMultiphysics/Kratos/issues/4833#issuecomment-491182913, or mute the thread https://github.com/notifications/unsubscribe-auth/AEYQZAFVFYJ5OYMJBHJN4NLPUUNSJANCNFSM4HLYL2QA .

[AlejandroCornejo]

Yes! I know that there is no direct solution but this is a "serious problem" because it can lead to perform 3 or 4 more iter or even diverge.. :)

[AndreasWinterstein]

What is the status here? @AlejandroCornejo with which type of elements are you encountering those problems?

[AlejandroCornejo]

I'm working with solid elements! I talked with @RiccardoRossi about this problem but we haven't reach a consensus about how to solve it..

[AlejandroCornejo]

Can we discuss this pls? I thinnk that this is importart :) @loumalouomega @RiccardoRossi

[philbucher]

@AlejandroCornejo why don't you give it a try by creating a scheme? I am quite sure that this should be not difficult to do, you basically have to compute the damping-matrix only once (I would say in InitializeSolutionStep) and then use it in the iterations instead of computing it (e.g CalculateSystemContributions would not calculate the DampingMatrix but get it from the stored value)

[loumalouomega]

@AlejandroCornejo why don't you give it a try by creating a scheme? I am quite sure that this should be not difficult to do, you basically have to compute the damping-matrix only once (I would say in InitializeSolutionStep) and then use it in the iterations instead of computing it (e.g CalculateSystemContributions would not calculate the DampingMatrix but get it from the stored value)

I think it is a good idea, except for the fact that the damping matrices are dense...(and you need to store it for each element). The schemes add D directly to the LHS, so you cannot differentiate what is D or K or M contributions.

I would create a class called "ConstantDampedElement" (or something like that): (so you don't need to rewrite or derivate all the elements). I would not register it, it should be something more on runing time

template<Element TElement>
class ConstantDampedElement : public TElement
....

void CalculateDampingMatrix(MatrixType& rDampingMatrix, ProcessInfo& rCurrentProcessInfo) override
{
    if (!mDampingInitialized) {
        TElement.CalculateDampingMatrix(mSavedDampingMatrix, rCurrentProcessInfo);
    }

    // Add resize and stuff
    rDampingMatrix = mSavedDampingMatrix;
}

....
MatrixType mSavedDampingMatrix; /// We save the damping matrix, it is expensive, but better than revert mesh moving....
bool mDampingInitialized = false; /// To check if the damping matrix is initialized 

[loumalouomega]

Then you need to replace your elements for this "constantly" damped elements

[loumalouomega]

But again, this is expensive, you are saving a dense matrix for each element...

[AlejandroCornejo]

Why don't we pass a FLAG like the ones used

    ConstitutiveLawOptions.Set(ConstitutiveLaw::COMPUTE_STRESS);
    ConstitutiveLawOptions.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR);

but with a new setting like:

ConstitutiveLawOptions.Set(ConstitutiveLaw::DAMPING_CALCULATION); that identificates wheter we enter the CalculateLeftHandSide to compute D or K?

[loumalouomega]

The problem is that dynamic systems those are multiply by the derivatives, I don't think it may work. Certainly ..., if we could save a sparse D matrix, it is cheaper to save. But then I don't know how to use it.

[AlejandroCornejo]

What do you mean? I don't understand you

[philbucher]

• I think we should really avoid having this in the elements, this has nothing to do with the elements (and will require implementation for each elements and hence is more effort and error prone
• It is also not related to the CLaw flags
• Since this is EXPERIMENTAL, I think for TESTING the easiest, fastest and least intrusive option is the scheme. If it gives good results then we can think abt integrating it better (e.g. a sparse damping matrix)

[loumalouomega]

* I think we should really avoid having this in the elements, this has nothing to do with the elements (and will require implementation for each elements and hence is more effort and error prone

It is an auxiliar class, you don't need to register anything, it is only called on running time. And indeed it has to do with the elements. Check the implementation of the schemes. The Damping matrices are local, and assembled all together on the LHS. You can manage them only in the element level.

• It is also not related to the CLaw flags True
• Since this is EXPERIMENTAL, I think for TESTING the easiest, fastest and least intrusive option is the scheme. If it gives good results then we can think abt integrating it better (e.g. a sparse damping matrix) OK, and how does the scheme to keep constant the contribution of the D matrix?, it is computed in a element level, the scheme cannot do anything about this.

[philbucher]

* I think we should really avoid having this in the elements, this has nothing to do with the elements (and will require implementation for each elements and hence is more effort and error prone

It is an auxiliar class, you don't need to register anything, it is only called on running time. And indeed it has to do with the elements. Check the implementation of the schemes. The Damping matrices are local, and assembled all together on the LHS. You can manage them only in the element level.

Ofc the element computes this, what I meant is that the element does not care if the problem is static, dynamic, quasi static, it just gives back what it is asked for (e.g. the LHS, the Damping-Matrix etc) => The element does not care how often the DampingMatrix is computed, whether in each timestep or in each iteration

• Since this is EXPERIMENTAL, I think for TESTING the easiest, fastest and least intrusive option is the scheme. If it gives good results then we can think abt integrating it better (e.g. a sparse damping matrix)

OK, and how does the scheme to keep constant the contribution of the D matrix?, it is computed in a element level, the scheme cannot do anything about this.

Yes, but the scheme can store it, what I was explaining in my comments above, the D is stored in InitializeSolutionStep and then used in the iterations (instead of recomputing)

[loumalouomega]

Yes, but the scheme can store it, what I was explaining in my comments above, the D is stored in InitializeSolutionStep and then used in the iterations (instead of recomputing)

You will need to store them in a map, with the ID of the element. This would be very expensive.

[philbucher]

I think for testing it would be ok

[AlejandroCornejo]

what is the state of this issue? I think we should reopen this :)

[philbucher]

have you tried my suggestion @AlejandroCornejo ?