Set a transverse isotrope/orthotropic constitutive law for SOVS

[alexisbu]

Hello,
I have a very basic understanding of mfront, and I am using it to model sintering through the castem interface. To that end, I used the SOVS webpage for mfront, and the paper of Balaguer et al. https://doi.org/10.1016/j.jeurceramsoc.2024.05.035 to adapt the proposed code to my material and adding his lambda (grain growth) parameter. I have done a few tutorial to be able to understand what was written, and Thomas H. gave me in the past some reading about mfront.
The current implementation is isotropic, and from my understanding, this is done by the line
dfeel_ddeel += 2 * mu * theta * dt * A * Stensor4::K()

I would like to set up a transverse isotropic behaviour, as the sintering strain is anisotropic in the Z direction.
I saw the discussion #469 but I am not sure that I should modify the StiffnessTensor as proposed in the shared code.

In papers related to this anisotropy, I saw that either we could modify the strain rate tensor (see Manière et al. [2020] https://doi.org/10.1016/j.actamat.2019.10.032) or the sintering stress tensor (see Abanobi et al [2025] https://doi.org/10.1016/j.jmapro.2025.11.030, but should the jacobian be re-written ?)

In brief my question is : what are the steps to follow to implement transverse isotropy in the code shared in the SOVS mfront web page ?

I zipped the working mfront and mtest files here
20251127_balaguer_lambda_val_optim.zip

Thanks for your help

Alexis B.

[thelfer]

Hi @alexisbu,
Would the material elasticity be transverse isotropic ?

[thelfer]

Hi @alexisbu,

With theta equal 1, one can see the computation of the jacobian is wrong:

mfront --obuild --interface=generic --@CompareToNumericalJacobian=true --@PerturbationValueForNumericalJacobianComputation=1e-8 --debug SOVS_sigmas_aniso_theta0_lamda.mfront 

[thelfer]

Instead of using the Hill tensor, I simply set the tensor sigS (which is supposed to represent an hypothetic stress linked to laplace pressure) with a different value for the ZZ component.

This only makes sense in the material frame, which you can't assume in general to be aligned with the global frame in general, no ?

[thelfer]

Tell me if this is easier for you by mail/skype or something.

Yes, that would help. I also have colleague interested by the topic

[alexisbu]

1. For theta=1, I did not get how you can see it using the command line proposed
2. Regarding the frame, I guess you must be right
3. ok let's do that

[thelfer]

For theta=1, I did not get how you can see it using the command line proposed

You shall get outputs like this:

SkoroAnisoTheta0::integrate() : iteration 19 : 1.66857e-12
2.35662 3.6e-13
dfeel_ddeel :
[[3.68426,-1.34213,-1.34213,0,0,0]
 [-1.34213,3.68426,-1.34213,0,0,0]
 [-1.34213,-1.34213,3.68426,0,0,0]
 [0,0,0,5.02639,0,0]
 [0,0,0,0,5.02639,0]
 [0,0,0,0,0,5.02639]]
ndfeel_ddeel :
[[3.94611,-1.08028,-1.08028,0,0,0]
 [-1.08028,3.94611,-1.08028,0,0,0]
 [-1.08028,-1.08028,3.94611,0,0,0]
 [0,0,0,5.02639,0,0]
 [0,0,0,0,5.02639,0]
 [0,0,0,0,0,5.02639]]
dfeel_ddeel - ndfeel_ddeel :
[[-0.261847,-0.261847,-0.261847,0,0,0]
 [-0.261847,-0.261847,-0.261847,0,0,0]
 [-0.261847,-0.261847,-0.261847,0,0,0]
 [0,0,0,0,0,0]
 [0,0,0,0,0,0]
 [0,0,0,0,0,0]]

• dfeel_ddeel is the jacobian that you have computed
• ndfeel_ddeel is the numerical approximation

You also seen the difference between the two values

[alexisbu]

After discussions and various trials, here is a working solution if sintering is considered anisotropic, with a larger shrinkage in the Z direction:
Still using the implementation proposed in the mfront SOVS webpage, we consider the sintering stress is a diagonal tensor with different values for the Z direction.

///// Definition of diagonal tensor sigS
stensor<N,real> sigS;
if constexpr (N==1){ 
   sigS  = {1, 1.2, 1};
   } else if constexpr (N==2){
   sigS = {1, 1.2, 1, 0};`
   } else {
   sigS = {1,1,1.2,0,0,0};
  }
// The value 1.2 is supposed to give 20% larger strain in Z than in X or Y.
// setting the value to 1 gives identical results than in the isotropic case

Then, the implicit equations and jacobians are changed as follows (A, B and C are unchanged)

// Implicit system
  feel = deel - deto + dt * A * s + trace(sig) * dt * B * id - C * dt * sigS;
  fe   = de - 3 * dt * B * trace(sig) + C*trace(sigS)*dt;

// Jacobian
  dfeel_ddeel += 2 * mu * theta * dt * A * Stensor4::K()+dt*theta*B*(3*lambda+2*mu)*(id^id);
  dfeel_dde = dt * dA_dde * s +dt*trace(sig)*dB_dde*id - dt*sigS*dC_dde;

  dfe_ddeel =  -3 * (3 * lambda + 2 * mu) * theta * B * dt * id;
  dfe_dde = 1- 3 * dt * dB_dde * trace(sig) + dt*trace(sigS)*dC_dde;

This implementation was tested for now using Castem for axisymetrical and plane strain conditions.

Additionnally, I am wondering whether the SOVS webpage contains a mistake when it is said

"The implicit system is then readily implemented:

  feel += de/3 * id + A * s;
  fe   -= 3 * dt * (B * trace(sig) - C);

.....
The previous lines are thus equivalent to

  feel = deel - deto + de/3 * id + A * s;
  fe   = de + 3 * dt * (B * trace(sig) - C);

I think the last equation (fe) should be with a "minus" before the 3*dt according to the -= above."

Work are on-going on this topic.

Thanks @thelfer for your help.

Best regards.

[thelfer]

I think the last equation (fe) should be with a "minus" before the 3*dt according to the -= above.

Yes, there is a clear typo here. Thanks !

[thelfer]

Fixed !

[thelfer]

Just a suggestion, I would replace this code

stensor<N,real> sigS;
if constexpr (N==1){ 
   sigS  = {1, 1.2, 1};
   } else if constexpr (N==2){
   sigS = {1, 1.2, 1, 0};`
   } else {
   sigS = {1,1,1.2,0,0,0};
  }

by:

auto sigS = Stensor{};
sigS[2] = 1.2;

[alexisbu]

Hello Thomas,
I observed some weird results and struggled finding the root cause of it. It seems I copy pasted your solution from github and not the one you sent me, so actually :

auto sigS = Stensor{};
sigS[2] = 1.2;

should be

auto sigS = Stensor::Id();
sigS[2] = 1.2; // for example if one wants an anisotropy around 20%

Best

[thelfer]

Sorry, my bad