Skip to content

Commit b805e7c

Browse files
thelferthelfer
authored
Merge pull request #72 from thelfer/71-add-massonaffineformulation-behaviours
71 add massonaffineformulation behaviours
2 parents cd9ff30 + 5d7e2c3 commit b805e7c

9 files changed

Lines changed: 614 additions & 1 deletion

File tree

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ It is mainly decomposed in two parts:
1919
- `plasticity`
2020
- `viscoelasticity`
2121
- `viscoplasticity`
22+
- `homogenization`
2223
- `materials`: this part gathers to specific materials.
2324

2425
The last purpose of this project is to show how to build a compilation

‎generic-behaviours/homogenization/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,3 +25,5 @@ genericmtest(Homogenization BetaRule
2525
genericmtest(Homogenization Idiart_viscoelastic
2626
BEHAVIOUR Idiart_viscoelastic
2727
REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Idiart_viscoelastic.ref)
28+
29+
add_subdirectory(MassonAffineFormulation)

‎generic-behaviours/homogenization/MassonAffineFormulation/Affine_tensors.mfront‎

Lines changed: 168 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,168 @@
1+
@DSL ImplicitII;
2+
@Behaviour Affine_tensors;
3+
@Author Martin Antoine;
4+
@Date 06 / 03 / 26;
5+
@Description{"Affine formulation for homogenization of a viscoplastic polycrystal, (Masson et al. 2001.), based on Green interaction tensors"};
6+
@UseQt false;
7+
@Algorithm NewtonRaphson;
8+
@Epsilon 1e-14;
9+
10+
11+
@Includes{
12+
#include "TFEL/Material/IsotropicModuli.hxx"
13+
#include "extra-headers/TFEL/Material/tensors.hxx"
14+
#include "TFEL/Material/PolyCrystalsSlidingSystems.hxx"}
15+
16+
//! number of phases
17+
@IntegerConstant Np =10;
18+
19+
//! normalization strainrate
20+
@Parameter real gamma0 = 1.0;
21+
//! bulk modulus of the reference medium used for the interaction tensors
22+
@Parameter stress mu0 = 1;
23+
24+
@ModellingHypothesis Tridimensional;
25+
@OrthotropicBehaviour;
26+
@CrystalStructure HCP;
27+
@SlidingSystems{<1, 1, -2, 0>{1, -1, 0, 0},
28+
<-2, 1, 1, 3>{1, -1, 0, 1},
29+
<-2, 1, 1, 0>{0, 0, 0, 1},
30+
<1, 1, -2, 0>{1, -1, 0, 1}};
31+
32+
@MaterialProperty real nexp;
33+
@MaterialProperty real tau1;
34+
@MaterialProperty real tau2;
35+
@MaterialProperty real kap;
36+
//@MaterialProperty real r0;
37+
38+
@StateVariable Stensor sigma[Np];
39+
sigma.setEntryName("PhaseReferenceStress");
40+
41+
@AuxiliaryStateVariable real ne;
42+
43+
@LocalVariable Stensor4 I;
44+
@LocalVariable Stensor4 J;
45+
@LocalVariable Stensor4 K;
46+
@LocalVariable Stensor4 L0;
47+
@LocalVariable Stensor4 M0;
48+
@LocalVariable real frac[Np];
49+
@LocalVariable real tau0[Np];
50+
@LocalVariable Stensor4 M[Np];
51+
@LocalVariable Stensor4 L[Np];
52+
@LocalVariable Stensor4 Ar[Np];
53+
@LocalVariable tmatrix<6*Np,6,real> A;
54+
@LocalVariable std::array<std::array<tfel::math::st2tost2<3u,real>,Np>,Np> GAMMA;
55+
@LocalVariable tmatrix<6*Np,6,real> dsigma_deto;
56+
@LocalVariable real r0;
57+
@LocalVariable real muhom;
58+
59+
@InitLocalVariables{
60+
I=tfel::math::st2tost2<3u,real>::Id();
61+
J=tfel::math::st2tost2<3u,real>::J();
62+
K=tfel::math::st2tost2<3u,real>::K();
63+
GAMMA=Gamma<real>::get_tensor();
64+
r0=1;
65+
L0=2*mu0*K;
66+
M0=1/(2*mu0)*K;
67+
r0=max(real(1),r0);
68+
}
69+
70+
71+
72+
@Integrator {
73+
ne=nexp;
74+
using ExtendedPolyCrystalsSlidingSystems =
75+
ExtendedPolyCrystalsSlidingSystems<Np, Affine_tensorsSlipSystems<real>, real>;
76+
const auto& gs =
77+
ExtendedPolyCrystalsSlidingSystems::getPolyCrystalsSlidingSystems("polycrystal.csv");
78+
79+
//tangent modulus and polarization/////////////////////////////////
80+
using namespace tfel::math;
81+
Stensor e[Np];
82+
Stensor4 Chom;
83+
Stensor dpsi_dsigma[Np];
84+
85+
for (int r=0;r<Np;r++){
86+
M[r]=Stensor4::zero();
87+
dpsi_dsigma[r]=Stensor::zero();
88+
frac[r]=gs.volume_fractions[r];
89+
for (int k=0;k<Nss;k++){
90+
auto tau0rk = tau1;
91+
if (k>int(Nss/12)){tau0rk=tau1;}
92+
else{tau0rk=tau2;}
93+
const auto Nkr = gs.mus[r][k]^gs.mus[r][k];
94+
const auto taukr = gs.mus[r][k] | (sigma[r]+dsigma[r]);
95+
const auto puisn_1 = pow(abs(taukr)/tau0rk, nexp-1);
96+
const auto fac= nexp*gamma0/tau0rk*puisn_1;
97+
M[r]+=fac*Nkr;
98+
const auto puisn = puisn_1*abs(taukr)/tau0rk;
99+
const auto sgn= taukr< 0 ? -real(1) : real(1);
100+
dpsi_dsigma[r]+=sgn*gamma0*puisn*gs.mus[r][k];
101+
}
102+
e[r]=dpsi_dsigma[r]-M[r]*(sigma[r]+dsigma[r]);
103+
M[r]=M[r]+kap*I;
104+
L[r]=invert(M[r]);
105+
}
106+
107+
//Operators A and B/////////////////////////////////
108+
tmatrix<6*Np,6*Np,real> MAT;
109+
tmatrix<6*Np,6*Np,real> G;
110+
tmatrix<6*Np,6,real> E;
111+
for (int r=0;r<Np;r++){
112+
map_derivative<Stensor,Stensor>(E,6*r,0)=I;
113+
for (int s=0;s<Np;s++){
114+
const auto dL = L[s]-r0*L0;
115+
map_derivative<Stensor,Stensor>(MAT,6*r,6*s) =1/r0*GAMMA[r][s]*dL;
116+
map_derivative<Stensor,Stensor>(G,6*r,6*s) =1/r0*GAMMA[r][s]*L[s];
117+
}
118+
map_derivative<Stensor,Stensor>(MAT,6*r,6*r) +=I;
119+
}
120+
121+
TinyMatrixInvert<6*Np,real>::exe(MAT);
122+
A = MAT*E;
123+
tmatrix<6*Np,6*Np,real> B = MAT*G;
124+
125+
Chom=Stensor4::zero();
126+
for (int r=0;r<Np;r++){
127+
Ar[r]=map_derivative<Stensor,Stensor>(A,6*r,0);
128+
Chom+=frac[r]*L[r]*Ar[r];
129+
}
130+
131+
132+
auto KGhom=tfel::material::computeKGModuli(Chom);
133+
muhom=KGhom.mu;
134+
r0=muhom/mu0;
135+
136+
auto tau_eff=Stensor::zero();
137+
//residues/////////////////////////////////////
138+
for (int r=0;r<Np;r++){
139+
tau_eff-=frac[r]*L[r]*e[r];
140+
fsigma[r] = sigma[r]+dsigma[r]-L[r]*Ar[r]*(eto+deto)+L[r]*e[r];
141+
for (int s=0;s<Np;s++){
142+
auto Brs = map_derivative<Stensor,Stensor>(B,6*r,6*s);
143+
tau_eff+=frac[r]*L[r]*Brs*e[s];
144+
fsigma[r]-=L[r]*Brs*e[s];
145+
}
146+
dfsigma_ddsigma(r,r)=I;
147+
}
148+
149+
//macroscopic stress//////////////////////////////////////
150+
sig=Chom*(eto+deto)+tau_eff;
151+
152+
}
153+
154+
155+
@TangentOperator{
156+
for (int r=0;r<Np;r++){
157+
map_derivative<Stensor,Stensor>(A,6*r,0)=L[r]*Ar[r];
158+
}
159+
tmatrix<6*Np,6*Np,real> iJ = jacobian;
160+
TinyMatrixInvert<6*Np,real>::exe(iJ);
161+
dsigma_deto=iJ*A;
162+
Dt=Stensor4::zero();
163+
for (int r=0;r<Np;r++){
164+
const auto dsigmar_deto = map_derivative<Stensor,Stensor>(dsigma_deto,6*r,0);
165+
Dt+=frac[r]*dsigmar_deto;
166+
}
167+
168+
}

‎generic-behaviours/homogenization/MassonAffineFormulation/Affine_tensors.mtest‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
@ModellingHypothesis 'Tridimensional';
2+
@Behaviour<@interface@> @library@ @behaviour@;
3+
@ExternalStateVariable 'Temperature' {0 : 1000,1:1000};
4+
5+
@MaterialProperty<function> 'nexp' "1+0.9*(exp(2.4*t)-1)";
6+
@MaterialProperty<constant> 'tau1' 5.;
7+
@MaterialProperty<constant> 'tau2' 1;
8+
9+
@MaterialProperty<function> 'kap' "0.05";
10+
11+
@ImposedStrain 'EXX' {0 : 1, 1 : 1};
12+
@ImposedStrain 'EYY' {0 : -1, 1 : -1};
13+
@ImposedStrain 'EZZ' 0;
14+
@ImposedStrain 'EXY' 0;
15+
@ImposedStrain 'EXZ' 0;
16+
@ImposedStrain 'EYZ' 0;
17+
@Times {0.,1 in 100};
18+
19+
@Test<file> @reference_file@ 'SXX' 8 1.e-3;
20+
@Test<file> @reference_file@ 'SYY' 9 1.e-3;
21+
@Test<file> @reference_file@ 'ne' 74 1e-5;

‎generic-behaviours/homogenization/MassonAffineFormulation/CMakeLists.txt‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
mfront_behaviours_library(MassonAffineFormulation
2+
Affine_tensors
3+
)
4+
5+
STRING(COMPARE EQUAL "${CMAKE_SOURCE_DIR}" "${CMAKE_BINARY_DIR}" _insource)
6+
if(NOT _insource)
7+
configure_file(polycrystal.csv polycrystal.csv)
8+
endif()
9+
10+
genericmtest(MassonAffineFormulation Affine_tensors
11+
BEHAVIOUR Affine_tensors
12+
REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Affine_tensors.ref)

0 commit comments

Comments
 (0)