Question about approach for DDSDDE calculation

[ashu28]

Hello all!
I am having some trouble formulating Jacobian for one of the plasticity models (Mulliken-Boyce model). Please have a look at it. I have simplified the equations to point to what I wanna ask:

F = Fe*Fp
Fe = Ve * Re
ε = ln (Ve)

T = f (ε)

T_back = g (F)

σ= T + T_back

Now, d(ΔT)/d(Δε) = Elastic Jacobian (Jac_1, say)

and that d(ΔT_back)/d(Δε) = Jac_2

So, DDSDDE = d(Δσ)/d(Δε) = Jac_1 + Jac_2

Am I correct in my understanding? I would really appreciate is someone can tell whether it is the right wau to calculate Jacobian...

Regards,
Ashu

[Jorgen]

Ashu,

Your approach for calculating the Jacobian is perfectly valid. When working with an advanced material model, such as the MB model, the tricky part of the Jacobian calculation is in the details of the implementation.

What technique do you intend to use to actually determine the terms of the Jacobian? A numerical approach? or closed-form derivations?

- Jorgen

[ashu28]

Hello Jorgen,
A very Happy New Year to you and everyone else.
At the current time step (t+dt), I update stress using explicit discretization (I didnt go for implicit because I coudnt figure out the discretization)
For Jac_1 (elastic Jacobian), I know the exact form: [λ+2µ 0 0 0 0 0 0;0 λ+2µ 0 0 0 0;0 0 0 λ+2µ 0 0 0;0 0 0 µ 0 0] (so on)

The issue remains to calculate Jac_2 (part coming from back stress) :
The model says:
T_back = c1*(1/lambda_p_chain)*Lang_inv(c2*lambda_p_chain)*B_bar_ dev_transpose
where, B_bar=(det(F)^(-2/3))FF'
lambda_p_chain = sqrt(trace(B_bar)/3)

At present I cant figure out how to get d(T_back)/dε. Although closed form would be exact, I would be satisfied if the numerical scheme converges fast enough as well. Right now I use only elastic jacobian and so it takes long to converge even with a single element.

Can you provide a few comments/starting steps on calculating d(T_back)/dε. I would really appreciate your help.

Looking forward to your reply.

[Jorgen]

There are basically two different ways to determine the Jacobian. In this case the complexity of the equations makes an attempt to derive an exact closed-form Jacobian quite challenging. Instead of persuing such a challenging task, you can derive closed-form approximate Jacobians with different levels of sophistication. You mentioned a linear elastic approach, which would be the easiest solution. As an alternative you can, for example, use a Neo-Hookean approach which can also be derived in closed-form. In order to improve the robustness of this approach you might want to take the current rate of viscoelastic flow into consideration, since that flow can significantly reduce the effective stiffness.

The other technique is to use a purely numerical approach in which you apply individual increments in strain in the various directions, and calculate the resulting stress increments. From that information you can then formulate an approximate Jacobian using a finite difference approach.

I have used both of these techniques with success when developing various user-material models.

- Jorgen

[ashu28]

Hello Jorgen,
Thanks a lot for the reply. Let me reiterate what you said about the numerical technique, and also write down what I understood.

Using finite difference approach:
Jacobian --> d(Δσ)/d(Δε) = [σ(t+dt) - σ(t)]/ [ε(t+dt) - ε(t)]

Now as I have already coded my constitutive model using explicit scheme, I know both σ and ε at both t and t+dt (storing the previous time step values in state variables). So I can evaluate each component of
d(Δσ_ij)/d(Δε_ij) using above scheme.

Is the expression written above correct ? and also does this approach significantly reduce the simulation time ?

I think I would go with this first as it seems easier.

Waiting for your advice.

Regards,
Ashu

[Jorgen]

The approach that you outlined is valid, although not the most numerically efficient. Note that you need to perform strain increments in various directions in order to vary only one term at a time for the numerical approximation.

- Jorgen

[ashu28]

Note that you need to perform strain increments in various directions in order to vary only one term at a time for the numerical approximation

Is this what you mean or something else which perhaps I did not get:

Cijkl = [σ_ij (t+dt) – σ_ij(t)] / [ε_kl(t+dt) – ε_kl(t)]
Now if stresses are denoted as [11 22 33 12] and so are strain (12 being shear)
and I take DDSDDE to be a 4x4 matrix this gives,
DDSDDE(1,1) = [σ_11 (t+dt) – σ_11(t)] / [ε_11(t+dt) – ε_11(t)]
DDSDDE(1,2) = [σ_11 (t+dt) – σ_11(t)] / [ε_22(t+dt) – ε_22(t)]
DDSDDE(1,4) = [σ_11 (t+dt) – σ_11(t)] / [ε_12(t+dt) – ε_12(t)]
DDSDDE(2,1) = [σ_22 (t+dt) – σ_22(t)] / [ε_11(t+dt) – ε_11(t)]
DDSDDE(4,4) = [σ_12 (t+dt) – σ_12(t)] / [ε_12(t+dt) – ε_12(t)]
and so on

Or to put in words, for the any row of Cijkl, stress indices are constant but the strains vary

Please correct me if I am wrong in my interpretation and understanding

Regards,
Ashu

[Jorgen]

Your experission might be right, although I am not quite sure I understand what you mean by [TeX:7d559bc00a]\epsilon_{kl}(t+dt)[/TeX:7d559bc00a].

Remember, one definition of C is:
[TeX:7d559bc00a]C_{ijkl} = \frac{\partial[\Delta\sigma_{ij}]}{\partial[\Delta\epsilon_{kl}]}[/TeX:7d559bc00a]
which means that you need to take the partial derivative, and hence only vary one of the strain components at a time. These strain perturbations do not have to be directly related to the applied strain increment [TeX:7d559bc00a]\epsilon_{ij}[/TeX:7d559bc00a].

- Jorgen

[ashu28]

Hello Jorgen,
Thanks for the suggestion. I am writing down what I understand from your reply. I should admit that I really think I might be completely wrong as this technique seems to require solving constitutive law for every C_ijkl component in itself. I am writing down the steps I'd need to follow:

Let the current time step values be denoted as t+dt and the previous time step values (obtained from state variables be t)

1. Let s1 = σ_11 at previous time step (obtained from state variables)
e1 = ε_11 at previos time step (obtained from state variables)
e2 = ε_11 at current time step

2. Get the total strain vector at the previous time step and replace the ε_11 component by ε_11 value at the current time step

3. Find σ corresponding to this new strain tensor and let s2 = σ_11 component of this new stress

4. C_{1111} = {s2-s1}/{e2-e1}


This has really confused me about how to discretize C_1111. I would be highly obliged if you can outline the steps required to obtain any one component of C matrix once.


Waiting for your reply,
Ashu

[Jorgen]

Ashu,

You are sort of right...

Here are a few more pointers:
- Yes, you should know the stress and strain at time t.
- If you then add a increment to e_ij then you can find one whole column in C_ijkl if you calculate the stress corresponding to that strain increment.
- The increments do not have to be related to the applied strain increments. This is important since otherwise you could run into problems with your Eq. (4) if e1==e1.

/Jorgen

[ashu28]

Okay, so inorder to avoid e1 (strain at time t+dt) = e1 (strain in previous step, ie time t), I stopped using the values of previous iteration altogether. Instead I approached it this way:

1. Take all the values of e and sigma at the current time step.
2. Construct a new e tensor by varying strain values one at a time, e_ij_new = e_ij_old (current time step tensor)*(1+x) where x is any small number, say 0.01
3. Find the stresses corresponding to this new strain tensor. This will give one column of Jacobian

The problem with this approach, or for that matter any approach where I do not use the previous time step values (for which I have strains and F) is that I need deformation gradient F to find out backstress, and inturn Sigma but I cant reconstruct F from e. So, even though, I perturb e to get the new strain tensor, I dont have F to get the backstress.

Can you advise on how to attack this problem ?

Regards,
Ashu

[Jorgen]

If you know the strain tensor, then you can calculate the left stretch tensor (v). Then, if you know the rotation tensor (R) from the increment that you just calculated in the UMAT, then you can calculate F that corresponds to that strain tensor.

Another approach for calculating the Jacobian, that you should be aware of, is that you can numerically determine the Jacobian not only from perturbations in the strain, but also directly from perturbations in F. If you want to use this approach then you start with the equation:
[TeX:e27b6846a9]\frac{d\tau}{dt} = J \mathbf{c} : \mathbf{d}[/TeX:e27b6846a9].

- Jorgen

[ashu28]

Thanks for the reply. I understood how to perturb e to get new strain vector and to use R of the current time step to calculate perturbed F.

One more question:

With new e and F, I can calculate backstress (which is a direct function of F), but to get stress(elastic), I need e_elastic or equivalently F_elastic. So the subrouting Jacobian itself becomes dependant on UMAT for splitting new F into Fe and Fp. How should I do that ? I hope you understand what I mean.


Regards,
Ashu

[Jorgen]

For a research-quality UMAT (as opposed to a commercial quality UMAT), you can assume that the strain perturpation takes place in zero time. That is, Fp does not change during the strain perturbations for the Jacobian matrix calculations.

- Jorgen