Incompressibility in VUMAT

[MarcoSteva]

Hello,

I am a graduating student working on a VUMAT for modeling the behaviour of the mitral valve. I need to enforce incompressibility to the model, but I'm not sure about the proper way to do it. I'm working with ABAQUS/Explicit on shell elements using an exponential strain-energy function. Do you have any suggestion about it?

Thanks,
Marco

C0 = props(1)
C1 = props(2)
C2 = props(3)

DO i=1,nblock

C11 = stretchNew(i,1)**2 + stretchNew(i,4)**2
C22 = stretchNew(i,2)**2 + stretchNew(i,4)**2
C12 = stretchNew(i,4)* ( stretchNew(i,1)+stretchNew(i,2) )
C33 = 1 / ( C11*C22 - C12**2 )

I1 = C11 + C22 + C33
I4 = C11

W1 = 2*C0*C1*(I1-3)* exp( C1*(I1-3)**2 + C2*(I4-1)**2 )
W4 = 2*C0*C2*(I4-1)* exp( C1*(I1-3)**2 + C2*(I4-1)**2 )
p = -2*W1*C33

S11 = 2*W1 + 2*W4 + p*(C22*C33)
S22 = 2*W1 + p*(C11*C33)
S12 = p*(C12*C33)

stressNew(i,1) = S11*defgradNew(i,1)**2 + S12*2*defgradNew(i,1)*defgradNew(i,4) + S22*defgradNew(i,4)**2

stressNew(i,2) = S11*defgradNew(i,5)**2 + S12*2*defgradNew(i,2)*defgradNew(i,5) + S22*defgradNew(i,2)**2

stressNew(i,3) = S11*defgradNew(i,1)*defgradNew(i,5) + S12*(defgradNew(i,1)*defgradNew(i,2)+defgradNew(i, 4)*defgradNew(i,5)) + S22*defgradNew(i,2)*defgradNew(i,4)

END DO

[Jorgen]

First, you cannot write a general purpose VUMAT for an incompressible material. However if you are only interested in shell elements that it is OK to use an incompressible material model.

There is nothing particularly different about using an incompressible material when working with a VUMAT for shell elements. The perhaps most difficult part when working with shell elements is to determine the through thickness strain (which is actually easier to determine when the material is incompressible!).

What aspects of the implementation do you have a problem with?

- Jorgen

[MarcoSteva]

Dear Jorgen,

first of all, thank you for your reply. I realize that in my previous letter my explanation of the problem was not exhaustive.

The VUMAT subroutine I am implementing is intended to be used only with shell elements and should describe an incompressible material characterized by preferentially oriented fibers within the its matrix, so that its overall behaviour can be described as transversely isotropic. Thus, from this point of view incompressibility should not be a problem. Besides, incrompessibility is hypothesized because the material we want to describe is actually nearly incompressible (it is a biological soft tissue) and because this hypothesis allows us to easily calculate the through thickness strain. Moreover, the model I took from the literature was reported to be successfully implemented in a UMAT subroutine, so that the model per se should not be an issue.

However, even if I am pretty sure that I have translated correctly the constitutive equations into my fortran code, the calculation of strains and stresses seems to be unstable. I say so on the basis of two problems I see on a simple test-case I am currently using to test the subroutine. It consists in an uniaxial traction of a patch of material in the preferential direction of the fibers. The traction is diplacement-controlled and the maximum applied deformation is 0.4.

The first problem is that the simulation I run goes on as I expect during its initial part, but when a critical point is reached a part of the mesh gets abruptly and hugely distorted and convergence is lost. Since I know that in Abaqus/Explicit the adopted time increment is a key factor, I repeated the same simulation with four different values of the maximum time increment set in my input file, these values ranging from 1E-07 to 5E-02. Surprisingly a gradual reduction of the maximum time-increment does not imply I gradually better convergence. On one hand, this evidence confuses me a little. On the other one, it suggests me that the key point could be at the beginning of the simulated process; in the last part of the simulation the adopted time increment goes down to 10E-07 in any case, so that using much higher values for the maximum time increment should only affect the initial iterations.

The second problem is that in three of the four simulations I run shear stresses and strains are not negligible, while they should be so given the boundary conditions of the problem. Could this outcome be due to round off errors or numerical instability? Moreover, given the nature of the constitutive model, the response of the material in the directions 2 and 3 should be the same (the material is transversely isotropic), but significant differences occurr between the two directions.

I do not know if you have ever had a similar problem in your experience, but I thought that, given your expertise, you might have an idea of what may be the critical aspects of the code. Besides, here is the part of the code that calculated the through thickness strain, which I did not include in the previous email.

Thank you again for your time and help. Sincerely,

Marco


C11 = stretchNew(i,1)**2 + stretchNew(i,4)**2
C22 = stretchNew(i,2)**2 + stretchNew(i,4)**2
C12 = stretchNew(i,4)* ( stretchNew(i,1)+stretchNew(i,2) )
C33 = 1 / ( C11*C22 - C12**2 ) ! incompressibility

stretchNew(i,3)=C33**(0.5)
strainInc(i,3) = stretchNew(i,3)-stretchOld(i,3)

[Jorgen]

It sounds to me like you are doing everything right (although I have to admit that I have not looked carefully at your code).

The problems that you describe should not occur, and my guess would be that you have some small but important error in your Fortran code :?

What I often do is to add print statements in my code at strategic locations in order to figure out what is going on when the subroutine is called. That way you can hopefully figure out where the problem comes from.

- Jorgen

[MarcoSteva]

Dear Jorgen,

I tried to run the simulation printing all the variables of interest. The results are quite weird. Everything goes straight until the 60%, then the deformation in direction 2 stops decreasing and starts to diverge exponentially, while deformations in direction 1 (the direction of the fibers, and the one i'm pulling) and 3 stay constant. In few iterations, the time increment dt falls to zero and the simulation exits because of an illegal floating point operation error. The Fortran code seems to be ok, and the numerical stability condition i found in literature is satisfied. I begin to suspect that the problem is a sort of instability within Abaqus/Explicit. I'm using version 6.5-4. Do you have any suggestion?

Thank you again. Sincerely,

Marco

[Jorgen]

I am not sure what is causing your problem, but here are a few suggestions (perhaps you have already tried most of these):

- Can you run your simulation with a newer version of ABAQUS (the current version if 6.7)?
- Do you use mass scaling? If so, can you use less mass scaling?
- How many increments does the simulation take before reaching the instability?
- Do you know that the instability is not a real geometric or material instability?

Jorgen