Merge pull request #928 from thelfer/906-add-second-order-estimate-from-ponte-castaneda-1996

906 add tuto on affine formulation for viscoplastic polycrystal (Masson 2000)

Author: antoinecea

bindings/python/tfel/material/IsotropicModuli.cxx

@@ -92,4 +92,6 @@ void declareIsotropicModuli(pybind11::module_& m) {
   declareKGModuli<double>(m, "KGModuli");
   declareYoungNuModuli<double>(m, "YoungNuModuli");
   declareLambdaMuModuli<double>(m, "LambdaMuModuli");
+  m.def("computeKGModuli", &tfel::material::computeKGModuli<double>);
+  m.def("computeIsotropicStiffnessTensor", &tfel::material::computeIsotropicStiffnessTensor<double>);
 }

bindings/python/tfel/material/LinearHomogenizationSchemes.cxx

@@ -14,6 +14,7 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 #include "TFEL/Material/LinearHomogenizationSchemes.hxx"
+#include "TFEL/Material/HomogenizationSecondMoments.hxx"
 
 template <tfel::math::ScalarConcept StressType>
 requires(tfel::math::checkUnitCompatibility<
@@ -233,4 +234,10 @@ void declareLinearHomogenizationSchemes(pybind11::module_& m) {
           return homogenization::elasticity::computeOrientedPCWScheme(
               IM, f, IMi, n_a, a, n_b, b, c, D);
         });
+  m.def("computeMeanSquaredEquivalentStrain", &homogenization::elasticity::computeMeanSquaredEquivalentStrain<double>,
+        pybind11::arg("IsotropicModuli_of_the_matrix"),
+        pybind11::arg("volume_fraction"),
+        pybind11::arg("IsotropicModuli_of_the_inclusion"),
+        pybind11::arg("squared_of_hydrostatic_macro_strain"),
+        pybind11::arg("squared_of_equivalent_macro_strain"));
 }

docs/web/BetaRule.md

@@ -192,7 +192,7 @@ The implementation requires the integration of the local behaviours.
 These are carried out via the `@BehaviourVariable` keyword.
 The implementation of the local behaviour is explained [here](./MericCailletaudSingleCrystalPlasticity.html).
 
-All files `MericCailletaudSingleCrystalViscoPlasticity.mfront`, `BetaRule.mfront` and `BetaRule.mtest` can be downloaded [here](./downloads/BetaRule.zip).
+All files `MericCailletaudSingleCrystalViscoPlasticity.mfront`, `BetaRule.mfront` and `BetaRule.mtest` are available in the `MFrontGallery` project, [here](https://github.com/thelfer/MFrontGallery/tree/master/generic-behaviours/homogenization/).
 
 
 For the example, we assume that the composites is made of only 2 phases.

docs/web/BiphasicLinearHomogenization.md

@@ -127,13 +127,13 @@ std::array<stress,2> tab_k={k0,ki};
 auto mu0=E0/2/(1+nu0);
 auto mui=Ei/2/(1+nui);
 std::array<stress,2> tab_mu={mu0,mui};
-auto pa2=computeIsotropicHashinShtrikmanBounds<3u,2u,stress>(tab_f,tab_k,tab_mu);
+auto pa2=computeIsotropicHashinShtrikmanBounds<3u,stress>(tab_f,tab_k,tab_mu);
 auto UB=std::get<1>(pa2);
 auto kh=std::get<0>(UB);
-auto muh=std::get<1>(UB)
-kg=KGModuli<stress>(kh,muh);
-auto Enu=kg.ToYoungNu();
+auto muh=std::get<1>(UB);
+kg=tfel::material::KGModuli<stress>(kh,muh);
 };
+auto Enu=kg.ToYoungNu();
 E=Enu.young;
 }
 ~~~~

docs/web/CMakeLists.txt

@@ -95,6 +95,7 @@ pandoc_html(Taylor "--toc" "--toc-depth=2")
 pandoc_html(Sachs "--toc" "--toc-depth=1")
 pandoc_html(BetaRule "--toc" "--toc-depth=2")
 pandoc_html(PonteCastaneda1992 "--toc" "--toc-depth=2")
+pandoc_html(MassonAffineFormulation "--toc" "--toc-depth=2")
 pandoc_html(Norton-web)
 pandoc_html(Norton-full)
 pandoc_html(coverity-scan)

docs/web/MassonAffineFormulation.md

@@ -75,27 +75,27 @@ In the variational approach, the homogenization problem is equivalent to a minim
 
 Ponte-Castaneda's idea is to use the dual function $f_r^*$ of $f_r$, under an hypothesis of concavity of $f_r$. Note that $w_r$ is convex relatively to $\tepsilon$, but the hypothesis is: $f_r$ concave relatively to $\epsiloneq^2$, which adds a supplementary restriction to the behaviour.
 \begin{aligned}
-    w_r(\tepsilon)=\underset{\mu_0^r(\tenseur x)}{\min} \left\{\Frac92 k_r\, \varepsilon_m^2+\Frac32 \mu_0^r(\tenseur x)\, \epsiloneq^2-f_r^*(\mu_0^r(\tenseur x))\right\}
+    w_r(\tepsilon)=\underset{\mu_0^r(\tenseur x)}{\min} \left\{\Frac92 k_r\, \varepsilon_m^2+\Frac32 \mu_0^r(\tenseur x)\, \epsiloneq^2-f_r^*(\dfrac32\mu_0^r(\tenseur x))\right\}
   \end{aligned}
   where the dual function $f_r^*$ is defined as
   \begin{aligned}
-    f_r^*(\mu_0)=\underset{e}{\min} \left\{\mu_0 e-f_r(e)\right\}
+    f_r^*(\lambda_0)=\underset{e}{\min} \left\{\lambda_0 e-f_r(e)\right\}
   \end{aligned}
 
 Considering uniform per phase shear moduli $\mu_0^r$ (also called the 'secant moduli'), and after a few manipulations, Ponte-Castaneda arrives at the following bound on the effective energy $W^{\mathrm{eff}}(\overline{\tepsilon})$:
 \begin{aligned}
-    W^{\mathrm{eff}}(\overline{\tepsilon})=\underset{\tepsilon\in\mathcal K(\overline{\tepsilon})}{\textrm{min}} \langle w(\tepsilon)\rangle \leq \overline{w}(\overline{\tepsilon})=\underset{\mu_0^r}{\min} \left\{W_0^{\mathrm{eff}}(\overline{\tepsilon})-\sum_r c_r\, f_r^*(\mu_0^r)\right\}
+    W^{\mathrm{eff}}(\overline{\tepsilon})=\underset{\tepsilon\in\mathcal K(\overline{\tepsilon})}{\textrm{min}} \langle w(\tepsilon)\rangle \leq \overline{w}(\overline{\tepsilon})=\underset{\mu_0^r}{\min} \left\{W_0^{\mathrm{eff}}(\overline{\tepsilon})-\sum_r c_r\, f_r^*(\dfrac32\mu_0^r)\right\}
   \end{aligned}
  where $c_r$ is the volume fraction of phase $r$ and $W_0^{\mathrm{eff}}(\overline{\tepsilon})$ is the effective energy relative to a 'linear comparison composite' which is heterogeneous, and whose elastic moduli are the $k_r$ and the $\mu_0^r$.  
 
   The stationarity conditions associated to the minimization on $\mu_0^r$ are:
   \begin{aligned}
-    \deriv{W_0^{\mathrm{eff}}}{\mu_0^r} = c_r\,\deriv{f_r^*}{\mu_0^r}
+    \deriv{W_0^{\mathrm{eff}}}{\mu_0^r} = \dfrac{3c_r}{2}\,\left(f_r^*\right)'(\dfrac32 \mu_0^r)
   \end{aligned}
 
 which can be shown to be equivalent to
 \begin{aligned}
-    \mu_0^r= \Frac23 \deriv{f_r}{e}\left(\langle \epsiloneq^2\rangle_r\right)
+    \mu_0^r= \Frac23 f_r'\left(\langle \epsiloneq^2\rangle_r\right)
   \end{aligned}
   where $\epsiloneq$ is relative to $\tepsilon$ solution of the homogenization problem:
 \begin{aligned}
@@ -111,34 +111,55 @@ which can be shown to be equivalent to
     \langle \epsiloneq^2\rangle_r= \Frac{2}{3c_r}\deriv{W_0^{\mathrm{eff}}}{\mu_0^r}
   \end{aligned}
 
-## Summary
-
-The resolution hence consists in 
-\[
- \text{Find}\,e^r=\langle \epsiloneq^2\rangle_r\,\text{such that}\quad\langle \epsiloneq^2\rangle_r= \Frac{2}{3c_r}\deriv{W_0^{\mathrm{eff}}}{\mu_0^r}
-\]
-where \(W_0^{\mathrm{eff}}\) is the effective energy of a linear comparison composite whose elastic moduli are \(k_r\) annd \(\mu_0^r\) such that:
-\[
-\mu_0^r= \Frac23\deriv{f_r}{e}\left(\langle \epsiloneq^2\rangle_r\right).
-\]
- 
 ## Macroscopic stress and tangent operator
 
 The macroscopic stress $\overline{\tsigma}$ is also shown to be (see [@ponte_castaneda_nonlinear_1998], eq. (4.41)):
   \begin{aligned}
-    \tsigma=\deriv{\overline{w}}{\overline{\tepsilon}}=\deriv{W_0^{\mathrm{eff}}}{\overline{\tepsilon}}=\tenseurq C_0^{\mathrm{eff}}\dbldot\overline{\tepsilon}
+    \overline{\tsigma}=\deriv{\overline{w}}{\overline{\tepsilon}}=\deriv{W_0^{\mathrm{eff}}}{\overline{\tepsilon}}=\tenseurq C_0^{\mathrm{eff}}\dbldot\overline{\tepsilon}
   \end{aligned}
 where $\tenseurq C_0^{\mathrm{eff}}$ is the effective elasticity of the linear comparison composite. In our case, we use Hashin-Shtrikman bounds.
 
 The tangent operator is given by
 \begin{aligned}
-    \dfrac{\mathrm{d}\tsigma}{\mathrm{d}\overline{\tepsilon}}=\tenseurq C_0^{\mathrm{eff}}+\dfrac{\mathrm{d}\tenseurq C_0^{\mathrm{eff}}}{\mathrm{d}\overline{\tepsilon}}
+    \dfrac{\mathrm{d}\overline{\tsigma}}{\mathrm{d}\overline{\tepsilon}}=\tenseurq C_0^{\mathrm{eff}}+\dfrac{\mathrm{d}\tenseurq C_0^{\mathrm{eff}}}{\mathrm{d}\overline{\tepsilon}}
   \end{aligned}
   Here, the derivative of $\tenseurq C_0^{\mathrm{eff}}$ w.r.t. $\overline{\tepsilon}$ can be computed by derivating the Hashin-Shtrikman
   moduli w.r.t. the secant moduli $\mu_0^r$ and the derivatives of these moduli w.r.t. $\overline{\tepsilon}$. However, it is tedious and in the implementation, we see that the convergence of the Newton Raphson algorithm is good if we only retain the first term \(\tenseurq C_0^{\mathrm{eff}}\) in the tangent operator.
+  
+  
+## Summary and possible implementations
+
+The resolution hence consists in 
+\[
+ \text{Find}\,e^r=\langle \epsiloneq^2\rangle_r\,\text{such that}\quad\langle \epsiloneq^2\rangle_r= \Frac{2}{3c_r}\deriv{W_0^{\mathrm{eff}}}{\mu_0^r}
+\]
+where \(W_0^{\mathrm{eff}}\) is the effective energy of a linear comparison composite whose elastic moduli are \(k_r\) annd \(\mu_0^r\) such that:
+\[
+\mu_0^r= \Frac23\deriv{f_r}{e}\left(\langle \epsiloneq^2\rangle_r\right).
+\]
+
+### Possible implementations
+
+The iterative resolution of the non-linear equation can be summarized as
+\[
+\underset{\substack{\\ \\ \\ \Downarrow\\ \\ \\\text{analytic / finite difference}\\ \\ \\\text{\texttt{@BehaviourVariable} / directly provided}}}{\mu_0^r= \Frac23\deriv{f_r}{e}\left(\langle \epsiloneq^2\rangle_r\right)}\quad\rightarrow\quad \underset{\substack{\\ \\ \\ \Downarrow\\ \\ \\\text{mean-field scheme / morphological tensors}\\ \\ \\\texttt{tfel::material::homogenization}}}{W_0^{\mathrm{eff}}\left(\mu_0^r\right)= \dfrac12\overline{\tepsilon}\dbldot\tenseurq C_0^{\mathrm{eff}}\dbldot\overline{\tepsilon}}\quad\rightarrow\quad\underset{\substack{\\ \\ \\ \Downarrow\\ \\ \\\text{analytic / finite difference}\\ \\ \\\texttt{tfel::material::homogenization}}}{\langle \epsiloneq^2\rangle_r= \Frac{2}{3c_r}\deriv{W_0^{\mathrm{eff}}}{\mu_0^r}}
+\]
+
+And the macroscopic stress must also be computed:
+\[
+\overline{\tsigma}=\tenseurq C_0^{\mathrm{eff}}\dbldot\overline{\tepsilon}
+\]
+
+The first step of the resolution consists in computing the secant modulus $\mu_0^r$ and it can be done analytically (as below) or via finite difference or automatic differentiation. Moreover, for both strategies, we could use a `BehaviourVariable`, defining the function $f_r$ in an external file (in our example, the function is to simple to do that).
+
+The second step consists in computing the effective energy $W_0^{\mathrm{eff}}$. This will be done via `tfel::material::homogenization`. This `namespace` provides classical mean-field schemes (we use a Hashin-Shtrikman bound in our example below). This could also be done with interaction tensors, providing externally. In another [tutorial](./MassonAffineFormulation.html) we show that we can take into account more specifically the morphology of a polycrystal by providing interaction tensors.
+
+The third step consists in the computation of the second-moments. This is also done via `tfel::material::homogenization`. Here again, there are two strategies: analytical computation when possible, and finite difference or automatic differentiation when the analytical derivation is too tedious or even impossible. In our example below, the computation is analytic.
 
 # Implementation in MFront
 
+## Example used for the implementation
+
 In the application of the implementation, we take the example of the following behaviour:
   \begin{aligned}
     w_r(\tepsilon)=\dfrac{9}{2}k_r\,\varepsilon_{m}^2+\dfrac{\sigma_r^0}{n+1}\,\varepsilon_{eq}^{n+1}
@@ -161,7 +182,7 @@ which means that we only need the derivative of the Hashin-Shtrikman bound. This
 
 ## Details of implementation
 
-All the files are available [here](./downloads/PonteCastaneda1992.zip)
+The file `PonteCastaneda1992.mfront` is available in the `MFrontGallery` project, [here](https://github.com/thelfer/MFrontGallery/tree/master/generic-behaviours/homogenization/).
 
 For the jacobian, we adopt the `Numerical Jacobian`, which
 means that the beginning of the `mfront` file reads:
@@ -171,7 +192,7 @@ means that the beginning of the `mfront` file reads:
 @Behaviour PC_VB_92;
 @Author Martin Antoine;
 @Date 12 / 12 / 25;
-@Description{"Ponte Castaneda second-order estimates for homogenization of non-linear elasticity (one potential), based on second-moments computation"};
+@Description{"Ponte Castaneda variational bounds for homogenization of non-linear elasticity (one potential), based on second-moments computation"};
 @UseQt true;
 @Algorithm NewtonRaphson_NumericalJacobian;
 @PerturbationValueForNumericalJacobianComputation 1e-10;
@@ -214,8 +235,9 @@ The second moment is given by the function `computeMeanSquaredEquivalentStrain`
 
 ~~~~ {#Integrator .cpp .numberLines}
 //second moments/////////////////////////////////
-const auto em2 = tfel::math::trace(eto+deto)/3.;
+const auto em = tfel::math::trace(eto+deto)/3.;
 const auto ed = tfel::math::deviator(eto+deto);
+const auto em2 = em*em;
 const auto eeq2 = 2./3.*(ed|ed);
 using namespace tfel::material;
 const auto kg0 = KGModuli<stress>(k_m,mu0);

docs/web/PonteCastaneda1992.md

@@ -94,6 +94,8 @@ The jacobian matrix is:
 
 # Implementation in MFront
 
+The files `Plasticity.mfront` and `Sachs.mfront` are available in the `MFrontGallery` project, [here](https://github.com/thelfer/MFrontGallery/tree/master/generic-behaviours/homogenization/).
+
 As for Taylor's scheme, we will use `BehaviourVariable`
 for the integration of local behaviours.
 The integration of global behavior, which requires solving a non-linear equation,

docs/web/Sachs.md

@@ -51,6 +51,8 @@ The macroscopic tangent operator is given by
 
 # Implementation in MFront
 
+The files `Plasticity.mfront` and `Taylor.mfront` are available in the `MFrontGallery` project, [here](https://github.com/thelfer/MFrontGallery/tree/master/generic-behaviours/homogenization/).
+
 It is possible to implement this homogenization scheme through a `Behaviour`.
 This `Behaviour` must call the behaviour laws of each phase.
 So, an interesting solution is to use a `BehaviourVariable`.