how to get material property dependant on strain (or stress)

[pl96andy]

Hi all (and thank you for project)

sorry for noob question but i could not find answers in doc

I want to have orthotropic material model, where some of material properties are dependant on actual strain (and are updated in each iteration)
As far as i understand i could use StandardElasticity brick and try to call @ElasticMaterialProperties with 9 variables, in which some would be function. But i do not understand how to get strain (i guess strain matrix) in each iteration.
Also is @ElasticMaterialProperties function updated on each iteration?

Thank you for any help
Andrzej

[thelfer]

Hi @pl96andy

Thanks for you interest to MFront.

The trouble is that if your elastic properties do depends on the total strain, then the hypothesis of the StandardElasticity brick are no more satisfied (see https://thelfer.github.io/tfel/web/HookeStressPotential.html and https://thelfer.github.io/tfel/web/implicit-dsl.html) and the computation of the consistent tangent operator would then be wrong.

Don't hesitate to provide more information about your constitutive equations, I could help (maybe).

Best,
Thomas

[pl96andy]

Thank you for your answer
soo i cannot do it easy way :D
you made me reread docs :D
ok let me describe whole story:
I came across quite elastic orthotropic material, with variable poisson ratios - poisson ratios depends linearily on strain in particular direction. Due to not so high elasticity of material but supposed to be simulated as bended thin sheet, i assume that using hecky law would be best approach
so in terms of equation i "just" need to modify inputs of orthotropic material:
E1=E1; E2=E2; E3=E3;
nu1=C1strain1+nu1; nu2=C2strain2 +nu2; nu3=C3*strain3+nu3;
shear_modulus12=shear_modulus12; shear_modulus23=shear_modulus23; shear_modulus13=shear_modulus13
where C1, C2, C3 are constants to be defined by user (nu is poisson ratio) and model is intended to be used in calculix
well that's the idea

soo my material model should be simmilar to this?:

@dsl Implicit;
@behaviour PhysicsBasedOrthotropic;
// @ModellingHypothesis

@MaterialProperty stress E1, E2, E3;
@MaterialProperty real nu12_0, nu13_0, nu23_0;
@MaterialProperty stress G12, G13, G23;
@MaterialProperty real C12, C13, C23;

nu21 := nu12E2/E1;
nu32 := nu23E3/E2;
nu13 := nu31*E1/E3;

@TangentOperator{

StiffnessTensor De;
Stensor4 iJe;
Stensor iJd;
computeUnalteredElasticStiffness<N,Type>::exe(De, lambda, mu);
getPartialJacobianInvert(iJe, iJd);
Dt = (1 - d) * De * iJe - ((De * eel) ^ iJd);

}

@ComputeStress{
// Build stiffness matrix in Voigt form
real C[6][6] = 0.0;

real D = 1 - nu12nu21 - nu23nu32 - nu31nu13 - 2nu12nu23nu31;

C[0][0] := E1*(1 - nu23nu32)/D;
C[1][1] := E2(1 - nu31nu13)/D;
C[2][2] := E3(1 - nu12*nu21)/D;

C[0][1] := E1*(nu21 + nu31nu23)/D;
C[1][0] := E2(nu12 + nu32*nu31)/D;

C[0][2] := E1*(nu31 + nu21nu32)/D;
C[2][0] := E3(nu13 + nu23*nu12)/D;

C[1][2] := E2*(nu32 + nu12nu31)/D;
C[2][1] := E3(nu23 + nu31*nu12)/D;

// Shear terms
C[3][3] := G23;
C[4][4] := G31;
C[5][5] := G12;

// Stress = C : strain (Voigt)
sig = Stensor::Zero();
for(int i=0; i<6; ++i)
for(int j=0; j<6; ++j)
sig[voigtIndex(i)] += C[i][j]*eps[voigtIndex(j)];
}

@Integrator{

nu12=C12strain1+nu12_0
nu23=C23strain2+nu23_0
nu13=C13*strain3+nu13_0
}

and yes this code is written by AI reinforced with natural stupidity ;)
best regards
Andrzej

[thelfer]

Hi @pl96andy,
Since everything is driven by the total strain, the Default DSL seems more appropriate. The constitutive laws seems somehow unsymmetric, so I am not sure if I understood everything correctly.
You'll find attach an possible implementation. For the moment, the consistent tangent operator is approximated by the initial elastic stiffness.
Andrzej.zip

[pl96andy]

Thank you very much
I think i might understand what is happening in this config:

@dsl Default; => works for both implicit and explicit?, lets FEA solver deal with it?
@ModellingHypothesis Tridimensional; => i looked into code and what is "undefined" hypothesis? will it detect which to use?
eto is strain matrix at current iteration
deto is matrix of prediction of change to be made for next iteration (end of current one)
e_ets[0], e_ets[1], e_ets[2], are predicted strain for end of this iteration in 1,2,3 directions
auto De = StiffnessTensor{}; is declaration of stiffnes tensor

  computeOrthotropicStiffnessTensor<hypothesis, unaltered, aoc, stress, real>(
      De, E1, E2, E3, nu12, nu23, nu13, G12, G13, G23);    

it will compute stiffness tensor but De is treated as an input? and what does stress and real men in this case?, is it similar to @ComputeStiffnessTensor?
sig = De * e_ets; => this is stress tensor from hookes law

hmm i guess i was trying to overshot with hecky law :D

Sorry if i am asking too much questions
Best regards
Andrzej

[thelfer]

@pl96andy

@dsl Default; works for both implicit and explicit?, lets FEA solver deal with it?

Yes. Every DSLs works in implicit and explicit solvers.

Note that the Implicit DSL is meant to integrate the consistutive equations using an implicit scheme, it has nothing to do with the solver

I looked into code and what is "undefined" hypothesis

We are talking about two different things: the argument of the @ModellingHypothesis keyword and the C++ enum defined in the TFEL/Material library.

The undefined hypothesis is defined in this enum. It is used internaly to denote all the modelling hypotheses by MFront. More precisely, variables can be defined for specific modelling hypothesess, but if we declare them with the "undefined" hypothesis, those variables will be defined for all modelling hypotheses.

eto is strain matrix at current iteration
deto is matrix of prediction of change to be made for next iteration (end of current one)
e_ets[0], e_ets[1], e_ets[2], are predicted strain for end of this iteration in 1,2,3 directions

Yes. but matrix is misleading. eto is a symmetric tensor and the storage convention is described here: https://thelfer.github.io/tfel/web/tensors.html#vector-notations-for-symmetric-tensors

Note that by default, the definition of eto depends on the solver. For small strain analysis, it is the linearised strain.
For finite strain analysis, it depends on the choice of the solver. For instance, for Abaqus, it will be defined by the choice of the Zaremba-Jaumann framework. For CalculiX, it is the Green-Lagrange strain.

But, you may change this by specifiying a strain measure in MFront. See mfront --help-keyword=Default:@StrainMeasure for details.

auto De = StiffnessTensor{}; is declaration of stiffnes tensor

Yes, but a this stage, all components are null

it will compute stiffness tensor but De is treated as an input?

This function takes De as a mutable reference and updates it

and what does stress and real men in this case?,

I don't know if it is meaningful here. stress, real, temperature are type aliases. But default, they are equivalent to the floating point type required by the interface used (double for Calculix). But if you add @UseQt true; then those type will be different and the C++ compiler will check the consistency of the operations used. See:

• mfront --help-keyword=Default:@UseQt
• https://thelfer.github.io/tfel/web/mfront-types.html

for details.

hmm i guess i was trying to overshot with hecky law :D

Don't know. For the moment, I need to be sure of what you want to do.

[pl96andy]

well in my case material is (thin) reinforced "rubber", which during tensile tests seems to behave like linearily elastic material, but is meant to be used for live hinge, where bending are rather large, thats why i am constantly looking if i by accident do not need to use hencky material model which would be a bit harder for me to grasp. But i guess native calculix approach of NLGEOM should be fine.
(by the way maybe i am just describing incorrectly what i am trying to achieve)

so s far as i understand (i will try to make my post better looking):

1. is is quite important to add failsafe of @StrainMeasure because without warning it could get unreasonable output - in my case it would be: @StrainMeasure GreenLagrange;
2. in this line, hypothesis is taken from @ModellingHypothesis Tridimensional; ?

computeOrthotropicStiffnessTensor<hypothesis, unaltered, aoc, stress, real>
because in ModellingHypothesis.hxx there is :
enum Hypothesis { AXISYMMETRICALGENERALISEDPLANESTRAIN, AXISYMMETRICALGENERALISEDPLANESTRESS, AXISYMMETRICAL, PLANESTRESS, PLANESTRAIN, GENERALISEDPLANESTRAIN, TRIDIMENSIONAL, UNDEFINEDHYPOTHESIS };

3. honestly i am confused where exactly during (calculix) calculation, mfront library is executed. because for example in calculix doc it is :

A linear elastic material, as implemented in CalculiX simulates a linear elastic relationship between the Lagrange strain and the Piola-Kirchhoff stress ofthe second kind (PK2)
so does this means that in mfront i need to present stress in Piola-Kirchhoff stress of the second kind, which in turn will be converted to cauchy stress by calculix? or is it dealt by solver interface
But honestly i am not sure if actually i would see any difference in mfront config
When i looked at umat_abaqusnl.f (delivered with calculix) it seemed that material model stress, strain and other states needs to be exported back and forth to calculix convention

4. so with using implicit DSL, somebody could even create his own solver? (offtopic)

[thelfer]

well in my case material is (thin) reinforced "rubber", which during tensile tests seems to behave like linearily elastic material, but is meant to be used for live hinge, where bending are rather large, thats why i am constantly looking if i by accident do not need to use hencky material model which would be a bit harder for me to grasp. But i guess native calculix approach of NLGEOM should be fine. (by the way maybe i am just describing incorrectly what i am trying to achieve)

For large rotations, small strain, the St Venant Kirchhoff law is sufficient. If large strain has to be considered, then we have several hyperelastic behaviours implemented.

in this line, hypothesis is taken from @ModellingHypothesis Tridimensional; ?

The behaviour is "instantiated" for every supported modelling hypotheses (defined by @ModellingHypothesis or @ModellingHypotheses), i.e. we generate the code for each hypothesis. In each instantiation, hypothesis as the corresponding value.

honestly i am confused where exactly during (calculix) calculation

MFront is called by calculix at each integration point during the assembly of the inner forces and stiffness matrix

Most behaviours in MFront are of the "strain-based" kink. Now we need to know what is strain measure used.

If NLGEOM is not set, then eto in MFront is the linearised strain (and @StrainMeasure shall not be set)

If NLGEOM is set, then eto in MFront depends if @StrainMeasure is set or not:

• If not set, Green-Lagrange strain is used
• If set to GreenLagrange, Green-Lagrange strain is used
• If set to Hencky, Hencky strain is used following the paper of Miehe et al. (Anisotropic additive plasticity in the logarithmic strain space: modular kinematic formulation and implementation based on incremental minimization principles for standard materials)

so with using implicit DSL, somebody could even create his own solver?

No, MFront only deal with what is happening at an integration point

[pl96andy]

I guess Code you have given me is more or less what i am looking for (i am not sure about details)

maybye i am just overthinking but if calculix returns green-lagrange strain (not real/engineering one), line: const auto nu12 = C12 * e_ets[0] + nu12_0; will lead to incorrect solution. Because i will try to calculate C12 constant from linear interpolation of testing machine results, while e_ets is different type of stiffness/unit

If i read google correctly, it might seem that line :
const auto nu12 = C12 * e_ets[0] + nu12_0;
for nl geom in calculix (green lagrange strain)should look like somewhat:
const auto nu12 = C12 * sqrt(2*e_ets[0] + 1) - 1+ nu12_0;
(and simmilar for rest of directions) am i right?
also performance-wise should not i define nu12 before @Integrator and only change its value within integrator so that this constant wont be recreated upon each iteration (but i am not sure about mfront syntax)?

Best Regards
Andrzej

[thelfer]

I am really sorry, but I still have troubles to understand the consistutive equations that you want to implement

[pl96andy]

Let me try to speak my mind
The constitutive model I want to implement is the following:

– Orthotropic (almost linear) elastic material
– Finite strain formulation (green lagrange or hencky)

– Strain measure: Green–Lagrange strain E
– Stress measure: second Piola–Kirchhoff stress S

calculix with nlgeom will give me Green-Lagrange strain
for calculix i should return PK2 stress (as i understand)

Standard equation of PK2 stress for green - lagrange stress: (yes i am undecided)
as far as i understand i should aim at hyperelastic definition:
S=∂W​/∂E
(W - energy density function)
but my material looks linear so i should just stick to simple:
S = C(E) : E
(C - stiffness tensors)

i hope i did not messed up those equations

stiffness tensor C depends on strain through the Poisson ratio:
ν12(E11) = ν12_0 + a1 · e11
ν12(E23) = ν23_0 + a2 · e23
ν12(E13) = ν13_0 + a3 · e13
a is linear constant measured from strain-stress tensile test data
Where e is normal strain, not green lagrange strain (of course alternatively i can convert a constant to be green lagrange strain compatible)
Of course i am not sure if for my case there is difference between E and e but it seems so (i think)

All other elastic constants are constant.

The goal is to implement this law in MFront and update the stiffness tensor
during the Newton iterations based on the strain.

I hope this less black box description would make more sense

my guess:
Andrzej-V2.mfront.txt

Best regards
Andrzej

[thelfer]

Except the last line, this seems like a valid MFront file.

However, I am still doubtful about the definition of the Poisson's ratios:

This seems more consistent:

ν12(E12) = ν12_0 + a1 · e12
ν12(E23) = ν23_0 + a2 · e23
ν12(E13) = ν13_0 + a3 · e13

However, in MFront:

• e12 is given by e_ets[3] (up to a factor sqrt(2), see https://thelfer.github.io/tfel/web/tensors.html#vector-notations-for-symmetric-tensors),
• e13 is given by e_ets[4]
• e23 is given by e_ets[5]

Please find attached a corrected version.

Andrzej-V2_TH.zip

[pl96andy]

well thank you
Honestly i cannot hind what is wrong with config - i copied this in raw.
than run test of it on ai created, one element simulation in calculix
than i run test on sample and calculix just crashes (no error code)
for a while it seemed to be almost ok but i cannot find any error in input file
My test file is simple tensile test
SAMPLE-v.inp.txt

edit: i asked linux kernel about errors and it seems that process has been killed due to exceeding RAM capacity... 20GB RAM