Stress-State outside yield surface in a implict solution

[Auubinno]

Hi all,

I would like to know if it is possible to have a stress-state outside the yield surface after the convergence of the solution through a implicit scheme. Note that the yield surface is considered as a StateVariable and I'm verifying the location of the yield-stress in comparison to the yield surface using the output from the mtest simulation. The stress is approximately 1 Pa outside the yield cap.

Since the solution is computed to reduce the state-variables to a certain value defined by the @epsilon, I was not expecting to have stress-states outside the yield surface.

[thelfer]

Dear @Auubinno,

Without an example, difficult to be affirmative.

In general, the equation associated with the equivalent plastic strain as the form:

fp = (seq - Rp) / young;

where Rp is the radius of the yield surface.

So, you shall have |fp|< epsilon after convergence.

Here are some reasons to see something different are:

• you did not specify that the theta parameter is 1, to the previous equation is satisfied at an intermediate time step (by default theta equals 0.5).
• you don't compare the same stress field (ex. PK2 in the behaviour, Cauchy in the post-processing)

[Auubinno]

Hi @thelfer,

Here the example I'm currently working on. The model is a modification of the modified cam clay to reproduce the unsaturated behavior of soils. The model use both the pre-consolidation stress and saturation as hardening parameters. The yield surface is defined by the following equation:

where p'x = pc ^ beta
beta is a function of saturation

I'm using bishop effective stresses to evaluate the yield surface. Below is the equation for bishop effective stress:

The X variable is identical to saturation (Sr) in my case.

In the @ComputeStress and @ComputeFinalStress block I'm calculating the total stress using the incremental form of bishop effective stress equation.

My mfront file looks like this:

@ComputeStress {
   const auto Se_ = (Sr - SrRes) / difSat;
   const auto dSe = theta*dSr / difSat;
   sig = sig0 + theta * dsig_deel * deel + (Se_ * dPCAP + dSe * PCAP) * Stensor::Id();

}

@Integrator{
    ...
    feel = deel + depl - deto;
    fy = (q * q - M2 * peff * (pceffBeta - peff)) / (pc0 * young);
    frpc = drpc + deplV * var1 * (rpcnew - rpc_min);
    fSr = dSr - (Rep | (deel + depl)) - Qep * dPCAP;
}

@ComputeFinalStress{
    const auto Se_ = (Sr - dSr - SrRes) / difSat;
    const auto dSe = dSr / difSat;
    sig = sig0 + dsig_deel * deel + (Se_ * dPCAP + dSe * PCAP) * Stensor::Id();
}

I'm using @Theta 1.0.
I'm don't really know what PK2 stress is. Not sure if that is my case though.

[thelfer]

Ok, since I only have part of the mfront file, I may miss something.
One way to debug this is to have an auxiliary state variable yield and to save the value of the yield surface, as follows:

@Integrator{
    ...
    feel = deel + depl - deto;
    yield = (q * q - M2 * peff * (pceffBeta - peff));
    fy = yield / (pc0 * young);
    frpc = drpc + deplV * var1 * (rpcnew - rpc_min);
    fSr = dSr - (Rep | (deel + depl)) - Qep * dPCAP;
}

yield will contain the value of the yield function at convergence of the Newton.

[Auubinno]

That's interesting. I'll give it a try and get back to you.

But in any case, the yield computed with the .res file should give the same result, right?

[thelfer]

But in any case, the yield computed with the .res file should give the same result, right?

That's precisely the question

[Auubinno]

Hi @thelfer,

Indeed the yield variable shows that the yield function is zero or negative. However, the value is not the same when computed using the result file from .mtest. The ouput of the simulation is attached below

histerese_Zhou_et_al_2012_fig17.zip

The simulation is conducted at constant isotropic stress (s1 = s2 = s3 = 20 kPa) and suction (PCAP) changes over time. Change in suction causes a change in effective stress leading to a change in volume, saturation and pre-consolidation stress over time (you will notice that in the output).

You can use the following python code to read the output and compute the value of yield function using the result file and compare to the YieldStress_Value auxiliary state variable.

path = <.res file location>
df = pd.read_csv(path, header=None)
columns = df.iloc[:29,0]
df = pd.read_csv(path, skiprows=29, header=None, sep=" ")
df.columns = columns

M = 1.12 # constant
Srsat = 1.0 # constant
Srres = 0.0 # constant
a1 = 1.0 # constant
Sr = dfhist.iloc[:,20]
p = - dfhist.iloc[:,9]/1000 # kPa
q = 0.0/1000.
pc = dfhist.iloc[:,24]/1000 # kPa
pwp = dfhist.iloc[:,25]/1000 # kPa
Se = (Sr - Srres)/(Srsat - Srres)
beta = 1.0 / (1.0 - pow(1.0 - Se, a1))
M2 = M*M
pl = p + pwp*Se #pl is the mean effective stress
pcbeta = pow(pc, beta)
x = (pl - pcbeta)/1 # if x>0 the stress-state is outside the yield surface
yielding = q*q - M2*pl*(pcbeta - pl)
df['yield'] = yielding

I can send the code by email. Can you share you email?
I would not like to make it public already since is part of my research. The plan is to make it available in the future.

Meanwhile, I have another question. I've noticed that during a timestep the code block @InitiateLocalVariable is computed multiple times. The terminal output is shown below:

Beginning of test suite 'MTest/histerese_Zhou_et_al_2012_fig17'
** no prediction
resolution from 0 to 0.001
####################  INITIATE LOCAL VARIABLES

MCCUNSAT_SIG_SAT::integrate() : beginning of resolution
MCCUNSAT_SIG_SAT::integrate() : iteration 0 : 8.13609e-06
MCCUNSAT_SIG_SAT::integrate() : iteration 1 : 3.30487e-09
MCCUNSAT_SIG_SAT::integrate() : iteration 2 : 4.93893e-16
MCCUNSAT_SIG_SAT::integrate() : convergence after 2 iterations

Dt :
1.56678e+07 -319490 -319490 0 0 0 
-319490 1.56678e+07 -319490 0 0 0 
-319490 -319490 1.56678e+07 0 0 0 
0 0 0 7.99366e+06 0 0 
0 0 0 0 7.99366e+06 0 
0 0 0 0 0 7.99366e+06 

iteration 1 : 5.86523e-06 88.1476 (5.86523e-06 5.86523e-06 5.86523e-06 0 0 0)
####################  INITIATE LOCAL VARIABLES

MCCUNSAT_SIG_SAT::integrate() : beginning of resolution
MCCUNSAT_SIG_SAT::integrate() : iteration 0 : 7.0601e-06
MCCUNSAT_SIG_SAT::integrate() : iteration 1 : 3.60951e-18
MCCUNSAT_SIG_SAT::integrate() : convergence after 1 iterations

Dt :
2.8e+07 1.2e+07 1.2e+07 0 0 0 
1.2e+07 2.8e+07 1.2e+07 0 0 0 
1.2e+07 1.2e+07 2.8e+07 0 0 0 
0 0 0 8e+06 0 0 
0 0 0 0 8e+06 0 
0 0 0 0 0 8e+06 

iteration 2 : 2.30884e-06 120.059 (3.55639e-06 3.55639e-06 3.55639e-06 0 0 0)
####################  INITIATE LOCAL VARIABLES

MCCUNSAT_SIG_SAT::integrate() : beginning of resolution
MCCUNSAT_SIG_SAT::integrate() : iteration 0 : 7.37054e-06
MCCUNSAT_SIG_SAT::integrate() : iteration 1 : 3.80074e-18
MCCUNSAT_SIG_SAT::integrate() : convergence after 1 iterations

Dt :
2.8e+07 1.2e+07 1.2e+07 0 0 0 
1.2e+07 2.8e+07 1.2e+07 0 0 0 
1.2e+07 1.2e+07 2.8e+07 0 0 0 
0 0 0 8e+06 0 0 
0 0 0 0 8e+06 0 
0 0 0 0 0 8e+06 

iteration 3 : 2.55524e-08 1.32873 (3.53084e-06 3.53084e-06 3.53084e-06 0 0 0)
####################  INITIATE LOCAL VARIABLES

MCCUNSAT_SIG_SAT::integrate() : beginning of resolution
MCCUNSAT_SIG_SAT::integrate() : iteration 0 : 7.37416e-06
MCCUNSAT_SIG_SAT::integrate() : iteration 1 : 3.80224e-18
MCCUNSAT_SIG_SAT::integrate() : convergence after 1 iterations

Dt :
2.8e+07 1.2e+07 1.2e+07 0 0 0 
1.2e+07 2.8e+07 1.2e+07 0 0 0 
1.2e+07 1.2e+07 2.8e+07 0 0 0 
0 0 0 8e+06 0 0 
0 0 0 0 8e+06 0 
0 0 0 0 0 8e+06 

iteration 4 : 2.82795e-10 0.0147053 (3.53056e-06 3.53056e-06 3.53056e-06 0 0 0)
####################  INITIATE LOCAL VARIABLES

MCCUNSAT_SIG_SAT::integrate() : beginning of resolution
MCCUNSAT_SIG_SAT::integrate() : iteration 0 : 7.3742e-06
MCCUNSAT_SIG_SAT::integrate() : iteration 1 : 3.80224e-18
MCCUNSAT_SIG_SAT::integrate() : convergence after 1 iterations

Dt :
2.8e+07 1.2e+07 1.2e+07 0 0 0 
1.2e+07 2.8e+07 1.2e+07 0 0 0 
1.2e+07 1.2e+07 2.8e+07 0 0 0 
0 0 0 8e+06 0 0 
0 0 0 0 8e+06 0 
0 0 0 0 0 8e+06 

iteration 5 : 3.12975e-12 0.000162747 (3.53056e-06 3.53056e-06 3.53056e-06 0 0 0)
####################  INITIATE LOCAL VARIABLES

MCCUNSAT_SIG_SAT::integrate() : beginning of resolution
MCCUNSAT_SIG_SAT::integrate() : iteration 0 : 7.3742e-06
MCCUNSAT_SIG_SAT::integrate() : iteration 1 : 3.80224e-18
MCCUNSAT_SIG_SAT::integrate() : convergence after 1 iterations

Dt :
2.8e+07 1.2e+07 1.2e+07 0 0 0 
1.2e+07 2.8e+07 1.2e+07 0 0 0 
1.2e+07 1.2e+07 2.8e+07 0 0 0 
0 0 0 8e+06 0 0 
0 0 0 0 8e+06 0 
0 0 0 0 0 8e+06 

iteration 6 : 3.46377e-14 1.80116e-06 (3.53056e-06 3.53056e-06 3.53056e-06 0 0 0)
convergence, after 6 iterations, order 1

resolution from 0.001 to 0.002
####################  INITIATE LOCAL VARIABLES

This behavior is not what I understood reading the Mfront Book as shown in the figure below. Figure 10.1 makes me understand that the code block @InitiateLocalVariable will only be computed again after @ComputeFinalStress, @UpdatedAuxiliaryStateVariable and @TangentOperator. Since @UpdatedAuxiliaryStateVariable is only computed at the end of the timestep, the expected behavior would be to compute @InitiateLocalVariable only after the all iterations from resolution 0 to 0.001 have been finished. Can you explain what am I missing here?

[thelfer]

The figure shows in which order the code blocks are called at each call of the behaviour. In particular, @InitLocalVariables is the first user defined code block called at each call to the behaviour. The arrow explicits that order.

Figure 10.1 makes me understand that the code block @InitiateLocalVariable will only be computed again after @ComputeFinalStress, @UpdatedAuxiliaryStateVariable and @TangentOperator.

No, that's not how this figure shall be read. I'll try to add explicit label like "call the behaviour" and "exit of the behaviour" to make it clearer.

[thelfer]

Indeed the yield variable shows that the yield function is zero or negative.

That means that the stress used at the final iteration for the evaluation of the yield surface is different from the final stress.

You can check this by introducing another auxiliary state variables sig_iter to save the stress at the last iteration as follows:

@AuxiliaryStateVariable StressStensor sig_iter;

@Integrator {
   sig_iter = sig;
   .....
}

If this confirms a difference you must try to understand why there is such a différence. Is @ComputeFinalStress inconsistent with @ComputeStress ?

[thelfer]

I've just notice that the external state variable PCAP have the suction value at the end of the time step while running Newton iteration

I don't think that this statement is true. In @Integrator, PCAP shall have the value at the beginning of the time step.

Everything is documented, but the simpliest way to be sure of what is done by MFront is to display the code blocks after the treatement by MFront. Consider the following file: ImplicitNorton.zip

$ mfront-query --code-blocks ImplicitNorton.mfront
 - ComputeFinalThermodynamicForces
- ComputeFinalThermodynamicForcesCandidate
- ComputeTangentOperator
- ComputeThermodynamicForces
- InitializeLocalVariables
- Integrator
$ mfront-query --code-block=ComputeFinalThermodynamicForces ImplicitNorton.mfront 
- ComputeFinalThermodynamicForces: (no description available)
  - used variables: T eel lambda mu sig 
  - code:
  #line 70 "ImplicitNorton.mfront"
  std::cout << (this->T+this->dT) << std::endl;
  #line 71 "ImplicitNorton.mfront"
  this->sig = this->lambda*trace(this->eel)*Stensor::Id()+2*this->mu*this->eel;

Here, @ComputeFinalStress was implemented as:

@ComputeFinalStress {
  std::cout << T << std::endl;
  sig = lambda*trace(eel)*Stensor::Id()+2*mu*eel;
}

You can see T was replaced by this->T+this->dT, which is the value of the temperature at the end of the time step.

If you have a look at the Integrator code block, you will see that T is replaced by this->T which is the value at the beginning of the time step:

#line 54 "ImplicitNorton.mfront"
  std::cout << (this->T) << std::endl;

[Auubinno]

You're right. I thought that both @ComputeStress and @Integrator used T at the beginning of the time step. But now I see that the @Integrator uses the value of T at the beginning while @ComputeStress and @ComputeFinalStress use the value of T at the end.
My understand was not correct because of miss interpretation of the flowchart from Figure 10.1.

[thelfer]

Hope this helps you to correctly set up your simulations.
Conventions in @ComputeFinalStress are indeed involved, but there was no good solution, so I tried to take the better. Hopefully, almost nobody needs to use @ComputeFinalStress (except when doing incremental elasticity, of course)

[Auubinno]

Thank you for your help! I'll review the implementation and correct accordingly.

To be honest, the context of implementing a numerical model is fairly new to me and the only way to implement without the use of @ComputeFinalSress that I'm aware of is using @brick StandardElasticity. Is there another way to do it?

[thelfer]

In most cases, only implementing @ComputeStress is enough as MFront will interpolate the values (at t+theta*dt or t+dt) but that does not work for incremental laws...