diff --git a/bindings/python/tfel/material/IsotropicModuli.cxx b/bindings/python/tfel/material/IsotropicModuli.cxx index c2f8a56e9a..fb440e34d6 100644 --- a/bindings/python/tfel/material/IsotropicModuli.cxx +++ b/bindings/python/tfel/material/IsotropicModuli.cxx @@ -92,4 +92,6 @@ void declareIsotropicModuli(pybind11::module_& m) { declareKGModuli(m, "KGModuli"); declareYoungNuModuli(m, "YoungNuModuli"); declareLambdaMuModuli(m, "LambdaMuModuli"); + m.def("computeKGModuli", &tfel::material::computeKGModuli); + m.def("computeIsotropicStiffnessTensor", &tfel::material::computeIsotropicStiffnessTensor); } diff --git a/bindings/python/tfel/material/LinearHomogenizationSchemes.cxx b/bindings/python/tfel/material/LinearHomogenizationSchemes.cxx index f4fafea920..60de9cabd7 100644 --- a/bindings/python/tfel/material/LinearHomogenizationSchemes.cxx +++ b/bindings/python/tfel/material/LinearHomogenizationSchemes.cxx @@ -14,6 +14,7 @@ #include #include #include "TFEL/Material/LinearHomogenizationSchemes.hxx" +#include "TFEL/Material/HomogenizationSecondMoments.hxx" template 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, + 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")); } diff --git a/docs/web/BetaRule.md b/docs/web/BetaRule.md index 4ccf3d6924..01d9ad1983 100644 --- a/docs/web/BetaRule.md +++ b/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. diff --git a/docs/web/BiphasicLinearHomogenization.md b/docs/web/BiphasicLinearHomogenization.md index a3b731c442..9f0e43432a 100644 --- a/docs/web/BiphasicLinearHomogenization.md +++ b/docs/web/BiphasicLinearHomogenization.md @@ -127,13 +127,13 @@ std::array tab_k={k0,ki}; auto mu0=E0/2/(1+nu0); auto mui=Ei/2/(1+nui); std::array 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(kh,muh); -auto Enu=kg.ToYoungNu(); +auto muh=std::get<1>(UB); +kg=tfel::material::KGModuli(kh,muh); }; +auto Enu=kg.ToYoungNu(); E=Enu.young; } ~~~~ diff --git a/docs/web/CMakeLists.txt b/docs/web/CMakeLists.txt index acc3e6a48e..a558debbf0 100644 --- a/docs/web/CMakeLists.txt +++ b/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) diff --git a/docs/web/MassonAffineFormulation.md b/docs/web/MassonAffineFormulation.md new file mode 100644 index 0000000000..72a7eafcef --- /dev/null +++ b/docs/web/MassonAffineFormulation.md @@ -0,0 +1,614 @@ +--- +title: Affine formulation for homogenization of a viscoplastic polycrystal +author: Antoine Martin +date: 30/01/2026 +lang: en-EN +link-citations: true +colorlinks: true +geometry: + - margin=2cm +papersize: a4 +bibliography: bibliography.bib +figPrefixTemplate: "$$i$$" +tblPrefixTemplate: "$$i$$" +secPrefixTemplate: "$$i$$" +eqnPrefixTemplate: "($$i$$)" +--- + +\newcommand{\tenseur}[1]{\underline{#1}} +\newcommand{\tenseurq}[1]{\underline{\mathbf{#1}}} +\newcommand{\tenseurt}[1]{\underline{\underline{\underline{#1}}}} +\newcommand{\tns}[1]{\underset{\tilde{}}{\mathbf{#1}}} +\newcommand{\transpose}[1]{#1^{\mathop{T}}} +\newcommand{\Ka}{\mathcal K} +\newcommand{\Sa}{\mathcal S} +\newcommand{\tsigma}{\tenseur{\sigma}} +\newcommand{\dbldot}{:} +\newcommand{\tepsilon}{\tenseur{\varepsilon}} +\newcommand{\te}{\tenseurq{e}} +\newcommand{\sigmaeq}{\sigma_{\mathrm{eq}}} +\newcommand{\epsiloneq}{\varepsilon_{\mathrm{eq}}} +\newcommand{\tE}{\tenseur E} +\newcommand{\tSigma}{\tenseur \Sigma} +\newcommand{\tepsilonto}{\underline{\epsilon}^{\mathrm{to}}} +\newcommand{\tepsilonel}{\underline{\epsilon}^{\mathrm{el}}} + +\newcommand{\trace}[1]{\mathrm{tr}\paren{#1}} +\newcommand{\Frac}[2]{{\displaystyle \frac{\displaystyle #1}{\displaystyle #2}}} +\newcommand{\deriv}[2]{{\displaystyle \frac{\displaystyle \partial #1}{\displaystyle \partial #2}}} +\newcommand{\derivdeux}[2]{{\displaystyle \frac{\displaystyle \partial^2 #1}{\displaystyle \partial #2^2}}} +\newcommand{\pderiv}[2]{{\dfrac{\partial #1}{\partial #2}}} +\newcommand{\pderivdeux}[2]{{\dfrac{\partial^2 #1}{\partial #2^2}}} +\newcommand{\derivderiv}[3]{{\displaystyle \frac{\displaystyle \partial^2 #1}{\displaystyle \partial #2\,\displaystyle \partial #3}}} +\newcommand{\dtot}{{\mathrm{d}}} +\newcommand{\paren}[1]{\left(#1\right)} +\newcommand{\nom}[1]{\textsc{#1}} +\renewcommand{\div}{\mathrm{div}} + + +We present here an implementation of the affine formulation [@masson_affine_2000] +for the homogenization of a viscoplastic polycrystal, example which is treated in [@bornert_second-order_2001] +but with Ponte-Castaneda second-order estimates. + +Here, the idea is to show that an implementation of that procedure based +on morphological tensors computed by FFT, +on a given geometry of polycrystal. + +This tutorial first presents the homogenization problem, recalls the methodology +of the affine formulation, presents different the possible implementations, +and show the details of the `mfront` file. + +# The viscoplastic polycrystal + +## Behaviour + +We consider a polycrystalline material, which means that each phase $r$ is associated to a crystal with +corresponding slip systems $\tenseur \mu_k^r$ ($1\leq k\leq K$). The strain rate in each crystal $r$ is given by +\[ + \begin{aligned} + \dot{\tepsilon}=\sum_k \dot{\gamma}_k^r\, \tenseur \mu_k^r\qquad\text{with}\quad \tenseur\mu_k^r = \dfrac12\left(\tenseur n_k^r\otimes\tenseur m_k^r + \tenseur m_k^r\otimes\tenseur n_k^r\right) + \end{aligned} +\] + where $\dot{\gamma}_k^r$ is the shear strain rate on the $k^{th}$ slip system of crystal $r$, and is given as a + function of the shear stress $\tau_k^r=\tsigma\dbldot\tenseur \mu_k^r$ by means of a potential $\psi_k^r$: +\[ + \begin{aligned} + \dot{\gamma}_k^r= \deriv{\psi_k^r}{\tau}\left(\tau_k^r\right) + \end{aligned} +\] + The expressions above show that on phase (or crystal) $r$: +\[ + \begin{aligned} + \dot{\tepsilon}=\deriv{\psi_r}{\tsigma}\left(\tsigma\right)\qquad\text{with}\quad\psi_r \left(\tsigma\right)=\sum_k \psi_k^r\left(\tau_k^r\right) + \end{aligned} +\] + Now, the behaviour of the polycrystal is governed by a potential $\psi$: +\[ + \begin{aligned} + \dot{\tepsilon}=\deriv{\psi_r}{\tsigma}\left(\tsigma\right)\qquad\psi(\tsigma)= \sum_{r=1}^{N}\chi_r\,\psi_r (\tsigma) + \end{aligned} +\] +where $N$ is the number of phases (or crystals) and $\chi_r$ is characteristic function of phase $r$. + + In all the sequel, we just note $\tepsilon$ for $\dot{\tepsilon}$. + +## The non-linear homogenization problem + +We impose a strain rate $\tE$ to the polycrystal and look for the solution $\tepsilon,\tsigma$ such that: +\[ +\begin{aligned} + &\div\,\tsigma=\tenseur 0\\ + &\tepsilon\in\Ka(\tE)\\ + &\tepsilon = \deriv{\psi}{\tsigma}\left(\tsigma\right) +\end{aligned} +\] + +where we introduced the space of kinematically admissible fields $\mathcal K(\tE)$, depending on the boundary conditions used (in the implementation, we work with periodic boundary conditions). + +# The affine formulation + +## The affine linearization + +The idea is to linearize the behaviour around a reference stress $\tsigma^r$: +\[ + \tepsilon_r(\tsigma)\approx\tenseurq M_r\left(\tsigma^r\right)\dbldot\tsigma+\tenseur e^r +\] + where +\[ + \tenseurq M_r\left(\tsigma^r\right)=\deriv{\tepsilon_r}{\tsigma}\left(\tsigma^r\right)=\derivdeux{\psi_r}{\tsigma}\left(\tsigma^r\right)\qquad\text{and}\quad\tenseur e^r = \deriv{\psi_r}{\tsigma}(\tsigma^r)-\tenseurq M_r\left(\tsigma^r\right)\dbldot\tsigma^r +\] + For a given set of reference stresses, this affine behaviour can be viewed as a so-called "thermoelastic comparison composite", and this composite + can be homogenized. This leads to a linear problem of the form +\[ +\begin{aligned} +&\div\, \left(\tenseurq L_r\dbldot\left(\tepsilon-\tenseur e^r\right)\right)=\tenseur 0\\ +&\tepsilon\in\Ka\left(\tE\right)\\ +&\tenseurq L_r=\left[\derivdeux{\psi_r}{\tsigma}\left(\tsigma^r\right)\right]^{-1}\qquad\text{and}\qquad\tenseur e^r=\deriv{\psi_r}{\tsigma}\left(\tsigma^r\right)-\tenseurq L_r^{-1}\dbldot\tsigma^r +\end{aligned} +\] + + Resolving this problem gives the average on phase $r$ of the strain solution $\tepsilon$ as +\[ + \begin{aligned} + \langle\tepsilon\rangle_r = \tenseurq A_r\left(\tsigma^1,...,\tsigma^N\right)\dbldot\tE + \sum_s \tenseurq B_{rs}\left(\tsigma^1,...,\tsigma^N\right)\dbldot\tenseur e^s\left(\tsigma^s\right) + \end{aligned} +\] + where $\tenseurq A_r\left(\tsigma^1,...,\tsigma^N\right)$ and $\tenseurq B_{rs}\left(\tsigma^1,...,\tsigma^N\right)$ can be obtained by a homogenization procedure (mean-field scheme, FFT...). The stresses are obtained by substracting the strain $\tenseur e^r$ and multiplying by the moduli $\tenseurq L_r=\left(\tenseurq M_r\right)^{-1}$, if it exists: + \[ + \begin{aligned} + \langle\tsigma\rangle_r = \tenseurq L_r\left(\tsigma^r\right)\dbldot\tenseurq A_r\left(\tsigma^1,...,\tsigma^N\right)\dbldot\tE + \sum_s \tenseurq L_r\left(\tsigma^r\right)\dbldot\tenseurq B_{rs}\left(\tsigma^1,...,\tsigma^N\right)\dbldot\tenseur e^s\left(\tsigma^s\right) - \tenseurq L_r\left(\tsigma^r\right)\dbldot\tenseur e^r\left(\tsigma^r\right) + \end{aligned} +\] + +The last question is the choice of the reference stresses $\tsigma^r$ ($1\leq r\leq N$). A discussion in [@masson_affine_2000] leads to the simple assumption that these reference stresses are the averages per phase of the stress in the thermoelastic composite: +\[ +\begin{aligned} +\tsigma^r = \langle\tsigma\rangle_r +\end{aligned} +\] +where $\tsigma$ here stands for the stress solution of the thermoelastic homogenization problem. + + +## Macroscopic stress and tangent operator + + The macroscopic stress hence can be obtained by the classical relation: +\[ + \begin{aligned} +\tSigma=\sum_rc_r\,\tsigma^r=\sum_rc_r\,\langle\tsigma\rangle_r=\tenseurq C^{\mathrm{eff}}\dbldot\tE+\tenseur\tau^{\mathrm{eff}} + \end{aligned} + \] + where + \[ +\begin{aligned} + \tenseurq C^{\mathrm{eff}}=\sum_rc_r\,\tenseurq L_r\dbldot\tenseurq A_r\qquad\text{and}\qquad\tenseur\tau^{\mathrm{eff}}=\sum_{r,s}c_r\,\tenseurq L_r\dbldot\tenseurq B_{rs}\dbldot\tenseur e^s-\sum_{r}c_r\,\tenseurq L_r\dbldot\tenseur e^r + \end{aligned} +\] + + The tangent operator is obtained by derivation: +\[ + \begin{aligned} + \dfrac{\mathrm{d}\tSigma}{\mathrm{d}\tE}=\sum_r c_r\,\deriv{\tsigma^r}{\tE} + \end{aligned} +\] + +In the expressions of the tangent operator, we can compute the term $\deriv{\tsigma^r}{\tE}$ by means of the jacobian matrix, also done and explained in the documentation of the [Implicit DSL](https://thelfer.github.io/tfel/web/implicit-dsl.html#computation-of-the-consistent-tangent-operator). We explain it below. + +# Summary and possible implementations + +## Summary + +The resolution consists in +\[ + \text{Find}\quad\tsigma^r\quad\text{such that}\quad\tsigma^r = \langle\tsigma\rangle_r +\] +where \(\tsigma\) is the stress field solution of a homogenization problem of the form: +\[ +\begin{aligned} +&\div\, \left(\tenseurq L_r\dbldot\left(\tepsilon-\tenseur e^r\right)\right)=\tenseur 0\\ +&\tepsilon\in\Ka\left(\tE\right)\\ +&\tenseurq L_r=\left[\derivdeux{\psi_r}{\tsigma}\left(\tsigma^r\right)\right]^{-1}\qquad\text{and}\qquad\tenseur e^r=\deriv{\psi_r}{\tsigma}\left(\tsigma^r\right)-\tenseurq L_r^{-1}\dbldot\tsigma^r +\end{aligned} +\] + +## Possible implementations + +The iterative resolution of the non-linear system can be summarized as + +\[ +\begin{aligned} +&\underset{\substack{\\ \\ \\ \Downarrow\\ \\ \\\text{analytic / finite difference}\\ \\ \\\text{\texttt{@BehaviourVariable} / directly provided}}}{\tenseurq L_r= \left[\pderivdeux{\psi_r}{\tsigma}\left(\tsigma^r\right)\right]^{-1},\quad\tenseur e^r=\deriv{\psi_r}{\tsigma}\left(\tsigma^r\right)-\tenseurq L_r^{-1}\dbldot\tsigma^r}\qquad\rightarrow\qquad \underset{\substack{\\ \\ \\ \Downarrow\\ \\ \\\text{mean-field scheme / morphological tensors}\\ \\ \\\texttt{tfel::material::homogenization}}}{\langle\tsigma\rangle_r = \tenseurq L_r\dbldot\tenseurq A_r\dbldot\tE + \sum_s \tenseurq L_r\dbldot\tenseurq B_{rs}\dbldot\tenseur e^s - \tenseurq L_r\dbldot\tenseur e^r}\\ +&\\ +&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\underset{\substack{\\ \\ \\ \Downarrow\\ \\ \\\text{Analytical Jacobian / Numerical Jacobian}\\ \\ \\\texttt{NewtonRaphson NumericalJacobian}}}{f_{\tsigma^r}=\tsigma^r-\langle\tsigma\rangle_r} +\end{aligned} +\] + +The first step of the resolution consists in computing the moduli $\tenseurq M_r$ and free strains $\tenseur e^r$. 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 local behaviour in an external file (in our example, the potential is to simple to do that). + +The second step is the computation of the tensors $\tenseurq A_r$ and $\tenseurq B_{rs}$. This can be achieved by means of mean-field schemes (the `namespace` `tfel::material::homogenization` provides the good tensors). This can also be done via morphological tensors, computed on a given geometry of polycrystal. We will present this second possibility below. + +We note that the equation relative to the third step will give our residue. + + +## Computation of tensors Ar and Brs + +### The Lippmann-Schwinger equation discretized + +In fact, the thermoelastic problem can be rewritten as a Lippmann-Schwinger equation (see [@willis_bounds_1977;@castaneda_effect_1995]). +Let us write $\tenseur\tau=\tsigma-\tenseurq L_0\dbldot\tepsilon$ with $\tsigma=\tenseurq L\dbldot\left(\tepsilon-\tenseur e\right)$. We have $\tenseur\tau=\tenseurq\delta\tenseurq L\dbldot\tepsilon-\tenseurq L\dbldot\tenseur e$ with $\tenseurq\delta\tenseurq L=\tenseurq L-\tenseurq L_0$. We have, because, $\div\,\tsigma=0$, +\[ +\begin{aligned} +\tepsilon=\tE-\tenseurq \Gamma_0\left(\tenseur\tau\right) +\end{aligned} +\] +where $\tenseurq \Gamma_0$ is the Green operator associated to the elasticity $\tenseurq L_0$ (this elasticity must be chosen by the user). It gives (due to the expression of $\tenseur\tau$) +\[ +\begin{aligned} +\tepsilon=\tE-\tenseurq \Gamma_0\left(\tenseurq\delta\tenseurq L\dbldot\tepsilon-\tenseurq L\dbldot\tenseur e\right) +\end{aligned} +\] +We will not resolve this equation exactly, but we will consider the average of the fields: +\[ +\begin{aligned} +\langle\tepsilon\rangle_r + \sum_s\tenseurq \Gamma_{rs}\dbldot\tenseurq\delta\tenseurq L_s\dbldot\langle\tepsilon\rangle_s= \tE +\sum_s\tenseurq \Gamma_{rs}\dbldot\tenseurq L_s\dbldot\tenseur e^s +\end{aligned} +\] +where +\[ +\begin{aligned} +\tenseurq \Gamma_{rs} = \sum_j\langle\,\tenseurq\Gamma_0(\chi_s\,\tenseur s_j)\otimes\tenseur s_j\,\rangle_r +\end{aligned} +\] +is what we call a morphological tensor or an interaction tensor, and can be computed by FFT or FEM before the `mfront` integration. This tensor depends on the Green operator $\tenseurq \Gamma_0$, relative to the elasticity $\tenseurq L_0$, and a basis of symmetric second-order tensors $(\tenseur s_1,...,\tenseur s_d)$ ($d=6$ in 3D). +Depending on the number of phases $N$, the computation of the $\tenseurq \Gamma_{rs}$ is more or less costly, but it is achieved +at the beginning, once and for all. + +After inversion of the linear system below, whose unknowns are the $\langle\tepsilon\rangle_s$, we obtain +\[ +\begin{aligned} +\langle\tepsilon\rangle_r= \sum_t\tenseurq T_{rt}\dbldot\left[\tE +\sum_s\tenseurq \Gamma_{ts}\dbldot\tenseurq L_s\dbldot\tenseur e^s\right] +\end{aligned} +\] +where $\tenseurq T_{rt}$ is the linear operator issued from the linear resolution. +and we directly obtain the expressions of $\tenseurq A_r$ and $\tenseurq B_{rs}$: +\[ +\begin{aligned} +&\tenseurq A_r= \sum_t\tenseurq T_{rt}\\ +&\tenseurq B_{rs}= \sum_t\tenseurq T_{rt}\dbldot\tenseurq \Gamma_{ts}\dbldot\tenseurq L_s +\end{aligned} +\] + +### Choice of reference medium + +We note that the reference medium can be chosen in this approach. A natural approach is to take a $\tenseurq L_0$ isotropic near $\tenseurq L^{\mathrm{hom}}$. To that extent we will make use of the fact that for an elasticity $\tenseurq C_1=r_0\tenseurq C_0$ the strain Green operator $\tenseurq \Gamma_1$ relative to $\tenseurq C_1$ is given by +\[ +\tenseurq \Gamma_1=\dfrac{1}{r_0}\tenseurq \Gamma_0 +\] +Hence, we will compute the morphological tensors $\tenseurq \Gamma_{rs}$ for a given elasticity $\tenseurq C_0$, and we can change the reference medium by multiplication of $\tenseurq C_0$ by $r_0$ and division of $\tenseurq \Gamma_0$ by $r_0$. + +Moreover, in our case, it makes sense to take a reference medium $\tenseurq C_0=2\mu_0\tenseurq K$. Afterwards, we will update the elasticity by taking $\tenseurq C_1=2\mu_{\mathrm{hom}}\tenseurq K$, where $\mu_{\mathrm{hom}}$ is equal to the shear modulus given by the projection of $\tenseurq C_{\mathrm{hom}}$ on $\tenseurq K$, where $\tenseurq C_{\mathrm{hom}}$ is the homogenized tensors obtained at previous step of the Newton-Raphson algorithm. + + +# Implementation in MFront + +All the files are available in the `MFrontGallery` project, [here](https://github.com/thelfer/MFrontGallery/tree/master/generic-behaviours/homogenization/), under the name `MassonAffineFormulation`. + +## Example used for the implementation + + The geometry of our polycrystal is generated with [`merope`](https://github.com/MarcJos/Merope/), + with a Voronoi tesselation. + This geometry is saved as an array and will be used only for the approach based on morphological + tensors. The image below shows an example of such a polycrystal. There are here 10 crystals. The + volume fractions can also be computed on this microstructure + + ![Geometry of the polycrystal](./img/polycrystal.png){width=50%} + + The morphological tensors were computed using a `python` script with `numpy.fft`. + The orientations of the slip systems will be given by angles, stored in the file `polycrystal.csv`. + This file also contains, at the end of each line, the volume fraction associated to a crystal. + They are in our case the same as the one computed on the generated microstructure. + + The potential $\psi_k^r$ will be of the form + \[ +\begin{aligned} + \psi_k^r \left(\tau\right)= \dfrac{\dot{\gamma}_0\tau_0^k}{n+1}\left(\dfrac{|\tau|}{\tau_0^k}\right)^{n+1} + \end{aligned} +\] + where the strain rate $\dot{\gamma}_0$ and the creep exponent $n\geq 1$ will be chosen identical for all $k,r$, and the resolved shear stress $\tau_0^k$ depends on $k$. + +The derivatives of the potential are + \[ +\begin{aligned} + &\tepsilon_r(\tsigma)=\deriv{\psi_r}{\tsigma}=\sum_k \deriv{\psi_k^r}{\tau}\left(\tau_k^r\right)\,\tenseur\mu_k^r=\sum_k\mathrm{sgn}(\tau)\,\dot{\gamma}_0\left(\dfrac{|\tau_k^r|}{\tau_0^{k}}\right)^{n}\tenseur\mu_k^r\qquad\text{with}\quad \tau_k^r = \tsigma\dbldot\tenseur \mu_k^r\\ + &\tenseurq M_r(\tsigma)=\derivdeux{\psi_r}{\tsigma}=\sum_k \derivdeux{\psi_k^r}{\tau}\left(\tau_k^r\right)\,\tenseurq N_k^r=\sum_k \dfrac{n\,\dot{\gamma}_0}{\tau_0^{k}}\left(\dfrac{|\tau_k^r|}{\tau_0^{k}}\right)^{n-1}\tenseurq N_k^r,\qquad\text{with}\quad\tenseurq N_k^r=\tenseur \mu_k^r\otimes\tenseur \mu_k^r + \end{aligned} +\] + and then defining $\overline{\tau}_k^r=\tsigma^r\dbldot\tenseur \mu_k^r$ where $\tsigma^r$ is the reference stress, we have + \[ +\begin{aligned} + &\tenseurq M_r(\tsigma^r)=\sum_k \dfrac{n\,\dot{\gamma}_0}{\tau_0^{k}}\left(\dfrac{|\overline{\tau}_k^r|}{\tau_0^{k}}\right)^{n-1}\tenseurq N_k^r\\ + &\tenseur e^r(\tsigma^r)=\sum_k\left(\mathrm{sgn}(\overline{\tau}_k^r)-n\right)\,\dot{\gamma}_0\left(\dfrac{|\overline{\tau}_k^r|}{\tau_0^{k}}\right)^{n}\tenseur\mu_k^r + \end{aligned} +\] +We see here that the modulus $\tenseurq M_r$ is not invertible. Hence, in the following, we add a small tensor to $\tenseurq M_r$ in order to make it invertible (see also appendix B.1 of [@bornert_second-order_2001]): +\[ +\tenseurq L_r=\left[\tenseurq M_r+\kappa\tenseurq I\right]^{-1} +\] + +## Details of implementation + +### Headers + +The beginning of the `mfront` file reads: + +~~~~ {#Begin .cpp .numberLines} +@DSL ImplicitII; +@Behaviour Affine_tensors; +@Author Martin Antoine; +@Date 06 / 03 / 26; +@Description{"Affine formulation for homogenization of a viscoplastic polycrystal, Masson et al. 2001., with morphological tensors"}; +@UseQt false; +@Algorithm NewtonRaphson; +@Epsilon 1e-14; +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Moreover, we need to include a header file for handling the rotations of the grains with their +slip systems: + + - PolyCrystalsSlidingSystems.hxx + +and a file which contains the morphological tensors. +This file is present in a repository, let's say `extra-headers/TFEL/Material/` +located at the same place as the `.mfront` file. + +~~~~ {#Begin .cpp .numberLines} +@TFELLibraries {"Material"}; +@Includes{ +#include "tensors.hxx" +#include "TFEL/Material/IsotropicModuli.hxx" +#include "TFEL/Material/PolyCrystalsSlidingSystems.hxx"} +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The compilation hence will be done like that: + +~~~~ {#Begin .sh .numberLines} +mfront -I ../extra-headers/TFEL/Material --obuild --interface=generic Affine_tensors.mfront +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Of cours the tensors are computed before via FFT as explained above +(see the folder mentioned above for the computation). + +### Sliding systems + +The sliding systems will be generated using these lines (the same +procedure as for the tutorial on the [Berveiller-Zaoui](ExplicitBerveillerZaouiPolyCrystals.html) scheme): + +~~~~ {#Begin .cpp .numberLines} +@ModellingHypothesis Tridimensional; +@OrthotropicBehaviour; +@CrystalStructure HCP; +@SlidingSystems{<1, 1, -2, 0>{1, -1, 0, 0}, + <-2, 1, 1, 3>{1, -1, 0, 1}, + <-2, 1, 1, 0>{0, 0, 0, 1}, + <1, 1, -2, 0>{1, -1, 0, 1}}; +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +but the rotations of the grains will be performed later, in the `Integrator` +code block: + +~~~~ {#Integrator .cpp .numberLines} +using ExtendedPolyCrystalsSlidingSystems = + ExtendedPolyCrystalsSlidingSystems, real>; +const auto& gs = + ExtendedPolyCrystalsSlidingSystems::getPolyCrystalsSlidingSystems("polycrystal.csv"); +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +### Properties and variables + +The state variables are the stresses $\tsigma^r$ on each phase. + +~~~~ {#Begin .cpp .numberLines} +@StateVariable Stensor sigma[Np]; +sigma.setEntryName("PhaseReferenceStress"); +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Note also that `Np` is the number of phases (or crystals). We define it +before like that: + +~~~~ {#Begin .cpp .numberLines} +@IntegerConstant Np = 10; +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Note that among several local variables, one of them, for +morphological tensors, is + +~~~~ {#Begin .cpp .numberLines} +@LocalVariable std::array,Np>,Np> GAMMA; +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This matrix will store the morphological tensors. Indeed, in the +`InitLocalVariables` code block, we can load these coefficients, +which are in fact situated in the header `tensors.hxx` + +~~~~ {#Begin .cpp .numberLines} +GAMMA=Gamma::get_tensor(); +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +### Affine modulus and free-strain + +There is no difficulty in computing the tangent modulus and free-strain, +in the `Integrator` code block: + +~~~~ {#Integrator .cpp .numberLines} +using namespace tfel::math; +Stensor4 LSC; +Stensor e[Np]; +Stensor dpsi_dsigma[Np]; +auto tau0k = tau1; + +for (int r=0;rint(Nss/12)){tau0k=tau1;} + else{tau0k=tau2;} + const auto Nkr = gs.mus[r][k]^gs.mus[r][k]; + const auto taukr = gs.mus[r][k] | (sigma[r]+dsigma[r]); + const auto puisn_1 = pow(max(abs(taukr),1e-5)/tau0k, nexp-1); + const auto fac= nexp*gamma0/tau0k*puisn_1; + M[r]+=fac*Nkr; + const auto puisn = puisn_1*abs(taukr)/tau0k; + const auto sgn= taukr< 0 ? -real(1) : real(1); + dpsi_dsigma[r]+=sgn*gamma0*puisn*gs.mus[r][k]; + } + e[r]=dpsi_dsigma[r]-M[r]*(sigma[r]+dsigma[r]); + M[r]=M[r]+kap*I; + L[r]=invert(M[r]); +} +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Here, on line 12, we set $\tau_0^k=\tau_1$ for all gliding systems, +except for $k=0$ and $k=1$, for which we set $\tau_0^k=\tau_2$. + +Note that here, a tensor `kap*I` is added to the compliance `M[r]` line 21. +In our case, we set `kap` as a `MaterialProperty` that can +evolve with time (see `mtest` file below). +We can hence consider $\tenseurq L_r$, the inverse of $\tenseurq M_r+\kappa\tenseurq I$. + +### Computation of Ar and Brs + +The computation of tensors Ar and Brs is as follows + + +~~~~ {#Integrator .cpp .numberLines} +tmatrix<6*Np,6*Np,real> MAT; +tmatrix<6*Np,6*Np,real> G; +tmatrix<6*Np,6,real> E; +for (int r=0;r(E,6*r,0)=I; + for (int s=0;s(MAT,6*r,6*s) =1/r0*GAMMA[r][s]*dL; + map_derivative(G,6*r,6*s) =1/r0*GAMMA[r][s]*L[s]; + } + map_derivative(MAT,6*r,6*r) +=I; +} + +TinyMatrixInvert<6*Np,real>::exe(MAT); +A = MAT*E; +tmatrix<6*Np,6*Np,real> B = MAT*G; +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Here we see that we multiply the morphological tensors by `1/r0` in order to adapt the reference medium. +The `r0` will be updated after, depending on the `Chom` obtained. +The `dL` also is impacted and written as `L[r]-r0*L0`. + +### Residues, jacobian and macroscopic stress + +Here is the computation of the residue, the jacobian and +of macroscopic stress: + +~~~~ {#Integrator .cpp .numberLines} +Chom=Stensor4::zero(); +for (int r=0;r(A,6*r,0); + Chom+=frac[r]*L[r]*Ar[r]; +} + +auto KGhom=tfel::material::computeKGModuli(Chom); +muhom=KGhom.mu; +r0=muhom/mu0; + +for (int r=0;r(B,6*r,6*s); + tau_eff+=frac[r]*L[r]*Brs*e[s]; + fsigma[r]-=L[r]*Brs*e[s]; + } + dfsigma_ddsigma(r,r)=I; +} + +sig=Chom*(eto+deto)+tau_eff; +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The jacobian here is approximative and given by `dfsigma_ddsigma(r,r)=I`. +We also see that we update the factor `r0` with the shear modulus `muhom`, +obtained by projection of `Chom` on the tensor `K`. + + +### Tangent operator + +The computation of the tangent operator necessitates the computation of $\deriv{\tsigma^r}{\tE}$. +To do that, we use a technique widely used in `MFront`. +The idea is to write the non-linear system to solve as: +\[ +\begin{aligned} +\tepsilon_r\left(\tsigma^r\right) - \tenseurq L_r\left(\tsigma^r\right)\dbldot\tenseurq A_r\left(\tsigma^1,...,\tsigma^N\right)\dbldot\tE - \sum_s \tenseurq B_{rs}\left(\tsigma^1,...,\tsigma^N\right)\dbldot\tenseur e^s\left(\tsigma^s\right)+\tenseurq L_r\left(\tsigma^r\right)\dbldot\tenseur e^r\left(\tsigma^r\right)= f_{\sigma^r}(\tE,\tsigma^1,...,\tsigma^N)=0 +\end{aligned} +\] +and by derivation of $f_{\sigma^r}$ w.r.t. $\tE$: +\[ +\begin{aligned} +-\tenseurq L_r\left(\tsigma^r\right)\dbldot\tenseurq A_r\left(\tsigma^1,...,\tsigma^N\right)+\sum_s\tenseurq J_{rs}\dbldot\deriv{\tsigma^s}{\tE}=0 +\end{aligned} +\] +where $\tenseurq J_{rs}=\deriv{f_{\sigma^r}}{\tsigma^s}$ stands for the sub-block $(r,s)$ of the Jacobian. Hence +\[ +\begin{aligned} +\deriv{\tsigma^r}{\tE}=\sum_s\mathbf {iJ}_{rs}\dbldot\tenseurq L_s\left(\tsigma^s\right)\dbldot\tenseurq A_s\left(\tsigma^1,...,\tsigma^N\right) +\end{aligned} +\] +where $\mathbf {iJ}$ is the inverse of the Jacobian. The implementation is + +~~~~ {#tangent .cpp .numberLines} +@TangentOperator{ + for (int r=0;r(A,6*r,0)=L[r]*Ar[r]; + } + tmatrix<6*Np,6*Np,real> iJ = jacobian; + TinyMatrixInvert<6*Np,real>::exe(iJ); + dsigma_deto=iJ*A; + Dt=Stensor4::zero(); + for (int r=0;r(dsigma_deto,6*r,0); + Dt+=frac[r]*dsigmar_deto; + } +} +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Note the use on line 7 of the function `tfel::math::map_derivative` +which allows to extract (but also to fill) a block from a `tmatrix`. + +## Results + +For the tests, we use `mtest`. However, note that the strain we impose +is in fact a strain rate. The idea for us is to keep this strain rate +constant and to make vary the parameter $n$ with time. +We hence set it as a function which varies exponentially with time +(this is just a choice). + +Besides, we choose $\tau_1=5 MPa$, and for $\tau_2$ +we try $\tau_2=5MPa$ and $\tau_2=1MPa$. + +The `.mtest` file is the following: + +~~~~ {#Affine_formulation .mtest .numberLines} +@ModellingHypothesis 'Tridimensional'; +@Behaviour 'src/libBehaviour.so' 'Affine_tensors'; +@ExternalStateVariable 'Temperature' {0 : 1000,1:1000}; + +@MaterialProperty 'nexp' "1+0.9*(exp(2.4*t)-1)"; +@MaterialProperty 'tau1' 5.; +@MaterialProperty 'tau2' 1.; + +@MaterialProperty 'kap' "0.05"; + +@ImposedStrain 'EXX' {0 : 1, 1 : 1}; +@ImposedStrain 'EYY' {0 : -1, 1 : -1}; +@ImposedStrain 'EZZ' {0 : 0, 1 : 0}; +@ImposedStrain 'EXY' {0 : 0, 1 : 0}; +@ImposedStrain 'EXZ' {0 : 0, 1 : 0}; +@ImposedStrain 'EYZ' {0 : 0, 1 : 0}; +@Times {0.,1 in 50}; + +@OutputFilePrecision 14; +~~~~~~~~~~~~~~ + +Note that `kap` on line 9 is related to the small compliance +we add to $\tenseurq M_r$ in order to make it invertible. + +![Axial macroscopic stress as a function of the parameter $n$](./img/Affine_formulation.png) + + +In this figure we plotted the results when $\tau_2$ is equal to $5 MPa$ +and when it is equal to $1MPa$. +We plot the implementation of the affine approach with morphological tensors (AFF). +FFT results are also plotted (FFT). +We see that the +morphological approach gives results in good agreement with FFT computation. + +# References + + diff --git a/docs/web/PonteCastaneda1992.md b/docs/web/PonteCastaneda1992.md index 2e21af6535..0e93324db5 100644 --- a/docs/web/PonteCastaneda1992.md +++ b/docs/web/PonteCastaneda1992.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(k_m,mu0); diff --git a/docs/web/Sachs.md b/docs/web/Sachs.md index cf7035258d..cfb216911d 100644 --- a/docs/web/Sachs.md +++ b/docs/web/Sachs.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, diff --git a/docs/web/Taylor.md b/docs/web/Taylor.md index 561fb98873..5ce135087c 100644 --- a/docs/web/Taylor.md +++ b/docs/web/Taylor.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`. diff --git a/docs/web/bibliography.bib b/docs/web/bibliography.bib index e2a213b2b5..19b3912a5c 100644 --- a/docs/web/bibliography.bib +++ b/docs/web/bibliography.bib @@ -12607,3 +12607,92 @@ @book{ponte_castaneda_nonlinear_1998 author = {Ponte Castaneda, Pedro and Suquet, Pierre}, year = {1998}, } + +@article{castaneda_exact_1996, + title = {Exact second-order estimates for the effective mechanical properties of nonlinear composite materials}, + volume = {44}, + copyright = {https://www.elsevier.com/tdm/userlicense/1.0/}, + issn = {00225096}, + url = {https://linkinghub.elsevier.com/retrieve/pii/0022509696000154}, + doi = {10.1016/0022-5096(96)00015-4}, + language = {en}, + number = {6}, + urldate = {2026-01-15}, + journal = {Journal of the Mechanics and Physics of Solids}, + author = {Castañeda, P.Ponte}, + month = jun, + year = {1996}, + pages = {827--862}, +} + +@article{bornert_second-order_2001, + title = {Second-order estimates for the effective behaviour of viscoplastic polycrystalline materials}, + volume = {49}, + abstract = {This paper is concerned with the application of the “second-order” nonlinear homogenisation procedure (Ponte Casta˜neda, J.Mech.Phys.Solids 44 (6) (1996) 827) to generate estimates of the self-consistent type for the e'ective behaviour of fcc and hcp viscoplastic polycrystals. The method has the distinctive property that it leads to estimates that are exact to second-order in the heterogeneity contrast, and which are expected to be more accurate, particularly when compared to rigorous bounds, than those resulting from earlier homogenisation schemes such as the Hill “incremental” method or its “total” formulation (Hutchinson) for pure power-law viscous materials.Special attention is paid to large grain anisotropy leading to correspondingly large heterogeneity contrast, and to highly nonlinear behaviour.Comparisons are also carried out with estimates derived from other more recent homogenisation schemes such as the “tangent” and “a7ne” methods.The results, illustrated for zirconium- and ice-type polycrystals, show that the second-order procedure o'ers the potential for signi9cantly improved results, at least relative to the Hill incremental formulation. ? 2001 Elsevier Science Ltd.All rights reserved.}, + language = {en}, + journal = {Journal of the Mechanics and Physics of Solids}, + author = {Bornert, M and Masson, R and Castaneda, P Ponte and Zaoui, A}, + year = {2001}, + pages = {2737--2764}, +} + +@article{masson_affine_2000, + title = {An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals}, + volume = {48}, + abstract = {Variational approaches for nonlinear elasticity show that Hill's incremental formulation for the prediction of the overall behaviour of heterogeneous materials yields estimates which are too sti and may even violate rigorous bounds. This paper aims at proposing an alternative `aÅne' formulation, based on a linear thermoelastic comparison medium, which could yield softer estimates. It is Ærst described for nonlinear elasticity and specified by making use of Hashin±Shtrikman estimates for the linear comparison composite; the associated aÅne self-consistent predictions are satisfactorily compared with incremental and tangent ones for power-law creeping polycrystals. Comparison is then made with the second-order procedure (Ponte Castanƒ eda, P., 1996. Exact second-order estimates for the e ective mechanical properties of nonlinear composite materials. J. Mech. Phys. Solids, 44 (6), 827±862) and some limitations of the aÅne method are pointed out; explicit comparisons between di erent procedures are performed for isotropic, two-phase materials. Finally, the aÅne formulation is extended to history-dependent behaviours; application to the self-consistent modelling of the elastoplastic behaviour of polycrystals shows that it o ers an improved alternative to Hill's incremental formulation. 7 2000 Elsevier Science Ltd. All rights reserved.}, + language = {en}, + journal = {Journal of the Mechanics and Physics of Solids}, + author = {Masson, R and Bornert, M and Suquet, P and Zaoui, A}, + year = {2000}, + pages = {1203--1227}, +} + +@article{castaneda_effect_1995, + title = {The effect of spatial distribution on the effective behavior of composite materials and cracked media}, + volume = {43}, + copyright = {https://www.elsevier.com/tdm/userlicense/1.0/}, + issn = {00225096}, + url = {https://linkinghub.elsevier.com/retrieve/pii/002250969500058Q}, + doi = {10.1016/0022-5096(95)00058-Q}, + abstract = {Estimates of the Hashin-Shtrikman type are developed for the overall moduli of composites consisting of a matrix containing one or more populations of inclusions, when the spatial correlations of inclusion locations take particular “ellipsoidal” forms. Inclusion shapes can be selected independently of the shapes adopted for the spatial correlations. The formulae that result are completely explicit and easy to use. They are likely to be useful, in particular, for composites that have undergone a prior macroscopically uniform large deformation. To the extent that the statistics that are assumed may not be realized exactly, the formulae provide approximations. Since, however, they are derived as variational approximations for composites with some explicit statistics that are realizable, they are free from some of the drawbacks of competitor approximations such as that of Mori and Tanaka (1973 Acta Metall. 21,571-574), which can generate tensors of effective moduli which fail to satisfy a necessary symmetry requirement. The new formulae are also the only ones known that take explicit account, at least approximately, of inclusion shape and spatial distribution independently.}, + language = {en}, + number = {12}, + urldate = {2026-01-20}, + journal = {Journal of the Mechanics and Physics of Solids}, + author = {Castañeda, P}, + month = dec, + year = {1995}, + pages = {1919--1951}, +} + +@article{willis_bounds_1977, + title = {Bounds and self-consistent estimates for the overall properties of anisotropic composites}, + volume = {25}, + copyright = {https://www.elsevier.com/tdm/userlicense/1.0/}, + issn = {00225096}, + url = {https://linkinghub.elsevier.com/retrieve/pii/0022509677900229}, + doi = {10.1016/0022-5096(77)90022-9}, + language = {en}, + number = {3}, + urldate = {2026-01-20}, + journal = {Journal of the Mechanics and Physics of Solids}, + author = {Willis, J.R.}, + month = jun, + year = {1977}, + pages = {185--202}, +} + +@article{idiart_thermodynamic_2025, + title = {Thermodynamic potentials for viscoelastic composites}, + volume = {194}, + issn = {00225096}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0022509624004022}, + doi = {10.1016/j.jmps.2024.105936}, + abstract = {Explicit expressions for the free-energy and dissipation densities of viscoelastic composites at fixed temperature are proposed. The composites are comprised of an arbitrary number of distinct constituents exhibiting linear Maxwellian rheologies and distributed randomly at a length scale that is much smaller than that over which applied loads vary significantly. Central to their derivation is the recognition that any viscous deformation field can be additively decomposed into an irrotational field and a solenoidal field in such a way that variational approximations available for elastic potentials become applicative to viscoelastic potentials. The thermodynamic potentials conform to a generalized standard model with a finite number of effective internal variables with explicit physical meaning. Specific approximations of the Hashin–Shtrikman and the Self-Consistent types are worked out in detail. Under particular circumstances, these approximations may turn out exact. Macroscopic stress–strain relations and intraphase statistics of the stress field up to second order are also provided.}, + language = {en}, + urldate = {2025-12-03}, + journal = {Journal of the Mechanics and Physics of Solids}, + author = {Idiart, Martín I.}, + year = {2025}, + pages = {105936}, +} diff --git a/docs/web/downloads/BetaRule.zip b/docs/web/downloads/BetaRule.zip deleted file mode 100644 index 23f8fd2c75..0000000000 Binary files a/docs/web/downloads/BetaRule.zip and /dev/null differ diff --git a/docs/web/downloads/PonteCastaneda1992.zip b/docs/web/downloads/PonteCastaneda1992.zip deleted file mode 100644 index 28720ef73a..0000000000 Binary files a/docs/web/downloads/PonteCastaneda1992.zip and /dev/null differ diff --git a/docs/web/gallery.md b/docs/web/gallery.md index 15cf4cf1b1..7de1746985 100644 --- a/docs/web/gallery.md +++ b/docs/web/gallery.md @@ -96,15 +96,20 @@ behaviour. This framework is described - The implementation of [Sachs scheme](Sachs.html) and [Taylor scheme](Taylor.html) show how to implement basic homogenized bounds with any behaviour laws on the phases using `@BehaviourVariable`. +- The implementation of Idiart's scheme [@idiart_thermodynamic_2025] for a + maxwellian viscoelastic composite with oriented ellipsoid is available [here](https://github.com/thelfer/MFrontGallery/blob/master/generic-behaviours/homogenization/Idiart_viscoelastic.mfront) - The description of the implementation of a polycrystal behaviour based on the Berveiller-Zaoui homogenisation scheme using an explicit scheme is available. [here](ExplicitBerveillerZaouiPolyCrystals.html) - [\(\beta\)-rule](BetaRule.html) can also be implemented with `@BehaviourVariable` and hence any behaviour law on each phase. -- A tutorial on implementation of Ponte-Castaneda bound (1992) - for a non-linear elastic composite is detaile [here](PonteCastaneda1992.html). +- A tutorial on implementation of Ponte-Castaneda bound [@castaneda_new_1992] + for a non-linear elastic composite is detailed [here](PonteCastaneda1992.html). +- A tutorial on implementation of the [affine formulation](MassonAffineFormulation.html) + from Masson et al. [@masson_affine_2000] for a viscoplastic polycrystal. +# References diff --git a/docs/web/img/Affine_formulation.png b/docs/web/img/Affine_formulation.png new file mode 100644 index 0000000000..15708278e5 Binary files /dev/null and b/docs/web/img/Affine_formulation.png differ diff --git a/docs/web/img/ParticulateMicrostructure.png b/docs/web/img/ParticulateMicrostructure.png new file mode 100644 index 0000000000..6787b37099 Binary files /dev/null and b/docs/web/img/ParticulateMicrostructure.png differ diff --git a/docs/web/img/distributions.png b/docs/web/img/distributions.png new file mode 100644 index 0000000000..49e4dbbff3 Binary files /dev/null and b/docs/web/img/distributions.png differ diff --git a/docs/web/img/distributions_PCW.png b/docs/web/img/distributions_PCW.png new file mode 100644 index 0000000000..0802b860c3 Binary files /dev/null and b/docs/web/img/distributions_PCW.png differ diff --git a/docs/web/img/ellipsoide_C0.png b/docs/web/img/ellipsoide_C0.png new file mode 100644 index 0000000000..5dc14be0f3 Binary files /dev/null and b/docs/web/img/ellipsoide_C0.png differ diff --git a/docs/web/img/polycrystal.png b/docs/web/img/polycrystal.png new file mode 100644 index 0000000000..d77bb37493 Binary files /dev/null and b/docs/web/img/polycrystal.png differ diff --git a/docs/web/tfel-material.md b/docs/web/tfel-material.md index a1af547587..8f308f7913 100644 --- a/docs/web/tfel-material.md +++ b/docs/web/tfel-material.md @@ -433,6 +433,8 @@ to the tolerance `eps`. The homogenization functions are part of the namespace `tfel::material::homogenization`. A specialization for elasticity is defined: `tfel::material::homogenization::elasticity`. +Note that the functionalities below are also available in +the `Python` module (see the doc [here](tfel-python.html#the-tfel.material.homogenization-module)). ## Eshelby, Hill and localisation tensors @@ -456,14 +458,17 @@ The expressions of Eshelby tensor can be found in [@torquato_2002] for the spheroidal inclusions and in [@eshelby_1957] for the general ellipsoid (three different semi-axes). +![Ellipsoidal inclusion](./img/ellipsoide_C0.png){width=35%} + Now if we consider an ellipsoid whose elasticity is \(\tenseurq C_i\), embedded in an infinite homogeneous medium whose elasticity is \(\tenseurq C_0\), submitted to a external uniform strain field at infinity \(\tenseur E\), the strain field within the ellipsoid is uniform and given by -\(\tenseur \varepsilon = \tenseurq A:\tenseur E\) +\(\tenseur \varepsilon = \tenseurq A_i:\tenseur E,\qquad\text{where}\quad + \tenseurq A_i = \left[\tenseurq I + \tenseurq P_0:\left(\tenseurq C_i -\tenseurq C_0\right)\right]^{-1}\) -where \(\tenseurq A \) is the strain localisation (or concentration) tensor. +where \(\tenseurq A_i \) is the strain localisation (or concentration) tensor. ### Computation in isotropic reference medium @@ -487,7 +492,7 @@ The second one computes the Hill tensor for an axisymmetrical ellipsoid (or sphe The user must provides the normal vector `n_a` for the axis, and `e` for the aspect ratio. The third line computes the Hill tensor of a more general ellipsoid whose semi-axis lengths are `a`,`b`,`c`. The axis `a` is related to direction given by `n_a` and `b` is related to the -direction given by `n_b`, which must be normal to `n_a`. +direction given by `n_b`, which must be normal to `n_a` (see the figure above). An `IsotropicModuli` can also be passed for the elasticity, as follows: @@ -567,6 +572,8 @@ const auto A_C = computePlaneStrainLocalisationTensor(IM0,C_i,n_a,a,b); The header `AnisotropicEshelbyTensor.hxx` introduces the computation of the Eshelby tensors and Hill tensors of general ellipsoids embedded in an anisotropic medium. +The anisotropic elasticity \(\tenseurq C_0\) must be +expressed in the global basis. These tensors can be computed as follows: @@ -603,6 +610,11 @@ not isotropic, it must be provided in the local basis defined by `n_a,n_b`. ## Homogenization schemes for biphasic media +Note that the functionalities below are also available in +the `Python` module (see the doc [here](tfel-python.html#homogenization-schemes-in-biphasic-media)). +See also [here](BiphasicLinearHomogenization.html) a tutorial on the computation of homogenized schemes for biphasic +particulate microstructures. + Different classical mean-field homogenization schemes are implemented for biphasic media. These schemes are introduced by the header `LinearHomogenizationSchemes.hxx`. They only deal with isotropic matrices and locally isotropic inclusions @@ -615,7 +627,7 @@ The available schemes are: - dilute scheme - Ponte Castaneda and Willis scheme -Each scheme is based on the average of the localisation tensor \(\tenseur A \) +Each scheme is based on the average of the localisation tensor \(\tenseurq A_i \) defined above. This average is computed assuming different distributions of ellipsoids. Hence different cases are considered: @@ -625,6 +637,8 @@ of ellipsoids. Hence different cases are considered: - transverse isotropic distribution of orientations (one axis \(\tenseur n_a\) of the ellipsoid is fixed, the others are uniformly distributed in the transverse plane) +![The three distributions of orientations considered here: oriented, isotropic, and transverse isotropic](./img/distributions.png){width=100%} + Hence we can compute the homogenized stiffness returned by the available schemes. For example, for the distribution of spheres: @@ -657,6 +671,16 @@ be created by the user. It is defined by two vectors \(\tenseur n_a,\tenseur n_b Distribution D = {.n_a = n_a, .a = a, .n_b = n_b, .b = b, .c = c}; ~~~~ +![The three distributions of orientations can be considered with PCW scheme, but a big ellipsoid defines the spatial distribution of inclusions](./img/distributions_PCW.png){width=100%} + +Indeed, in Ponte-Castaneda and Willis scheme, there is a difference between +the ellipsoid which defines the distribution of the inclusions, and the +ellipsoid which defines the shape of the inclusions (in the image above, we represent +the distribution `D` of inclusions by a big ellipsoid, with colored axes). +A bigger axis for `D` means that along this axis, the distribution of inclusions +is more diluted. A short axis means that, on the contrary, the distribution is denser +along this axis. + For the isotropic distribution of ellipsoids, we can do: ~~~~{.cpp} @@ -669,7 +693,7 @@ Here, the two first schemes return `KGModuli` objects, whereas `computeIsotropicPCWScheme` returns a `st2tost2` object. For this latter case, the ellipsoids have indeed a uniform isotropic distribution of orientations, but the user might use a non-isotropic -`Distribution D`. +`Distribution D` (it corresponds to the configuration at the center on the image). And finally, we can consider a transverse isotropic distribution of inclusions: @@ -680,7 +704,7 @@ const auto C_PCW = computeTransverseIsotropicPCWScheme(IM0,f,IMi,n_a,a,b ~~~~ Here, the three above schemes return `st2tost2` objects. -Because the functions are based on the average of the localisation tensor \(\tenseur A \) +Because the functions are based on the average of the localisation tensor \(\tenseurq A_i \) associated with each distribution, a `Base` function is also defined for each scheme, that only takes in argument the average of the localisation tensor `A_av`. We then can compute a homogenized stiffness with a very general averaged localisator: @@ -692,9 +716,6 @@ const auto C_PCW = computePCWScheme(E0,nu0,f,Ei,nui,A_av,D); ~~~~ Here, the three above schemes return `st2tost2` objects. -A tutorial on the computation of homogenized schemes for biphasic particulate microstructures -is available [here](BiphasicLinearHomogenization.html). - ## Homogenization bounds @@ -746,7 +767,8 @@ The number of phases is arbitrary, and the dimension is 2 or 3. ## Second moments of the strains -Some functions allow to compute the second moments of the strains +Some functions (defined in the header HomogenizationSecondMoments.hxx) +allow to compute the second moments of the strains if considering a Hashin-Shtrikman type microstructure. More precisely, we consider isotropic spherical inclusions embedded in an isotropic matrix, and we compute the following moments: @@ -778,6 +800,8 @@ of the `pair` corresponds to a phase (matrix or spheres). A `ParticulateMicrostructure` object can be created for homogenization of general matrix-inclusion microstructures. +Note that the functionalities below are also available in +the `Python` module (see the doc [here](tfel-python.html#homogenization-of-general-microstructures)). ### Description and construction of a microstructure @@ -788,6 +812,8 @@ via 2 template parameters: `ParticulateMicrostructure` with `N` the dimension. For the details, see the file 'MicrostructureDescription.hxx' which introduces the `class`. +![The `ParticulateMicrostructure class` is made of a matrix which embeds different distributions of inclusions](./img/ParticulateMicrostructure.png){width=50%} + A `ParticulateMicrostructure` consists on a matrix, in which are embedded several distributions of inclusions. The class has three (private) attributes: @@ -1019,7 +1045,7 @@ It is `12` by default. A last method of the `ParticulateMicrostructure` object allows to replace the matrix phase: -~~~~{.py} +~~~~{.cpp} micro_1.replaceMatrixPhase(C0); std::cout<< micro_1.get_matrix_elasticity()<< std::endl; std::cout<< micro_1.is_isotropic_matrix()<< std::endl; @@ -1032,7 +1058,7 @@ Note that the 4 `InclusionDistribution` classes are currently available in 3d on ### Computation of homogenization schemes -The file `MicrostructureLinearHomogenization.ixx` introduces the `HomogenizationScheme class` +The file `MicrostructureLinearHomogenization.ixx` introduces a `HomogenizationScheme` object which has three attributes: - `homogenized_stiffness` @@ -1059,7 +1085,7 @@ auto hmSC=computeSelfConsistent<3u,stress>(micro_1,1e-6,true); We note that `computeSelfConsistent` not only takes the microstructure as an argument, but also takes one real (`1e-6`) as -a parameter, which pilots the precision of the result. Indeed, at each iteration +a parameter, which pilots the accuracy of the result. Indeed, at each iteration of the self-consistent iterative algorithm, the function computes the relative difference between the new and the old homogenized stiffness. This relative difference must be smaller than the tolerance given as a parameter. @@ -1068,16 +1094,16 @@ the computation considers an isotropic matrix when computing the Hill tensors relative to the inclusions, at each iteration of the algorithm. Indeed, the homogenized stiffness may be non isotropic, so that the user can make the choice of isotropizing this homogenized stiffness at each iteration. -Otherwise, he can put `False`, so that a numerical integration (resulting +Otherwise, he can put `false`, so that a numerical integration (resulting in a slower computation) will be performed to compute the Hill tensors. Moreover, an integer parameter can be added after the boolean, that indicates the number of subdivisions in the numerical integration. This value is `12` by default: -~~~~{.py} +~~~~{.cpp} micro_2.addInclusionPhase(distrib_O); -hmSC_iso=computeSelfConsistent<3u,stress>(micro_2,10,True); -hmSC_aniso=computeSelfConsistent<3u,stress>(micro_2,10,False,10); +hmSC_iso=computeSelfConsistent<3u,stress>(micro_2,10,true); +hmSC_aniso=computeSelfConsistent<3u,stress>(micro_2,10,false,10); std::cout<< "SC iso: "<< hmSC_iso.homogenized_stiffness<< std::endl; std::cout<< "SC aniso: "<< hmSC_aniso.homogenized_stiffness<< std::endl; ~~~~ @@ -1094,7 +1120,7 @@ Moreover, if anisotropic, another parameter can be passed to specify the number of subdivisions in the numerical integration (this value is `12` by default): -~~~~{.py} +~~~~{.cpp} micro_1.replaceMatrixPhase(C0); micro_1.removeInclusionPhase(0); micro_1.addInclusionPhase(distrib_O); diff --git a/docs/web/tfel-python.md b/docs/web/tfel-python.md index 76339469a0..86a49a66a1 100644 --- a/docs/web/tfel-python.md +++ b/docs/web/tfel-python.md @@ -156,6 +156,8 @@ import tfel.material as tmat import tfel.math as tm ~~~~ +![Ellipsoidal inclusion](./img/ellipsoide_C0.png){width=35%} + ### Hill tensors The Hill tensors are available: @@ -247,7 +249,7 @@ Note that if the elasticity of the inclusion is not isotropic, an anisotropic elasticity `C_i` can be provided, assuming that this elasticiy is expressed in the same basis as the one defined by `n_a,n_b` (the local basis of the inclusion): -~~~~{.cpp} +~~~~{.py} A_aniso = hm.computeLocalisationTensor(IM0,C_i,n_a,a,n_b,b,c) ~~~~ @@ -280,7 +282,7 @@ The following schemes are available for biphasic media with - dilute scheme - Ponte Castaneda and Willis scheme -Here are some examples of computation: +Here are some examples of computation for the spherical inclusions: ~~~~{.py} young=1e9 @@ -298,7 +300,16 @@ KGDS_IM=hm.computeSphereDiluteScheme(IM,f,IMi) KGMT_IM=hm.computeSphereMoriTanakaScheme(IM,f,IMi) print(EnuDS.young,EnuDS.nu,KGDS_IM.kappa,KGDS_IM.mu) print(EnuMT.young,EnuMT.nu,KGMT_IM.kappa,KGMT_IM.mu) +~~~~ + +And we can also consider distribution of ellipsoidal inclusions, +with three kind of distributions of orientations. +![The three distributions of orientations considered here: oriented, isotropic, and transverse isotropic](./img/distributions.png){width=100%} + +Hence, here are the examples to compute the homogenized properties: + +~~~~{.py .numberLines} # Ellipsoidal inclusions a=10. b=1. @@ -334,6 +345,15 @@ print(C_O_MT) print(C_O_PCW) ~~~~ +![The three distributions of orientations can be considered with PCW scheme, but a big ellipsoid defines the spatial distribution of inclusions](./img/distributions_PCW.png){width=100%} + +In Ponte-Castaneda and Willis scheme (PCW), there is a difference between +the ellipsoid which defines the distribution of the inclusions, and the +ellipsoid which defines the shape of the inclusions (in the image above, we represent +the distribution `D` of inclusions by a big ellipsoid, with colored axes). +A bigger axis for `D` means that along this axis, the distribution of inclusions +is more diluted. A short axis means that, on the contrary, the distribution is denser +along this axis. The object `D` is defined above at line 12. ### Homogenization bounds @@ -406,6 +426,8 @@ in the namespace `tfel::material::homogenization::elasticity` for the construction and homogenization of general microstructures. The reader may want to consult this documentation [here](tfel-material.html#homogenization-of-general-microstructures). +![The `ParticulateMicrostructure` object is made of a matrix which embeds different distributions of inclusions](./img/ParticulateMicrostructure.png){width=50%} + The `ParticulateMicrostructure` object is defined and can be instantiated in various ways: diff --git a/include/TFEL/Material/IsotropicModuli.hxx b/include/TFEL/Material/IsotropicModuli.hxx index 98a147a510..0beb99f45f 100644 --- a/include/TFEL/Material/IsotropicModuli.hxx +++ b/include/TFEL/Material/IsotropicModuli.hxx @@ -61,10 +61,7 @@ namespace tfel::material { requires( tfel::math::checkUnitCompatibility()) struct IsotropicModuli { - IsotropicModuli() = default; - IsotropicModuli(const IsotropicModuli& IM) = default; - virtual IsotropicModuli& operator=( - const IsotropicModuli& IM) = default; + virtual ~IsotropicModuli() = default; virtual YoungNuModuli ToYoungNu() const = 0; virtual LambdaMuModuli ToLambdaMu() const = 0; @@ -78,7 +75,9 @@ namespace tfel::material { : public IsotropicModuli { StressType young; types::real nu; - + + YoungNuModuli() = default; + YoungNuModuli(const YoungNuModuli& ) = default; YoungNuModuli(const StressType& Young, const types::real& Nu) : IsotropicModuli(), young(Young), nu(Nu) {} @@ -112,6 +111,9 @@ namespace tfel::material { : public IsotropicModuli { StressType kappa; StressType mu; + + KGModuli() = default; + KGModuli(const KGModuli& ) = default; KGModuli(const StressType& Kappa, const StressType& Mu) : IsotropicModuli(), kappa(Kappa), mu(Mu) {} @@ -146,6 +148,9 @@ namespace tfel::material { : public IsotropicModuli { StressType lambda; StressType mu; + + LambdaMuModuli() = default; + LambdaMuModuli(const LambdaMuModuli& ) = default; LambdaMuModuli(const StressType& Lambda, const StressType& Mu) : IsotropicModuli(), lambda(Lambda), mu(Mu) {} diff --git a/include/TFEL/Material/MicrostructureLinearHomogenization.ixx b/include/TFEL/Material/MicrostructureLinearHomogenization.ixx index 6c540bc823..514abf6515 100644 --- a/include/TFEL/Material/MicrostructureLinearHomogenization.ixx +++ b/include/TFEL/Material/MicrostructureLinearHomogenization.ixx @@ -63,7 +63,7 @@ namespace tfel::material::homogenization::elasticity { KG0 = computeKGModuli(C0); } const auto tau0 = polarisations_[0]; - HomogenizationScheme h_s; + auto Chom = C0; auto tau_eff = tau0; @@ -86,9 +86,7 @@ namespace tfel::material::homogenization::elasticity { tau_eff += fi * tfel::math::transpose(Ai) * (taui - tau0); localisators.push_back(Ai); } - h_s.homogenized_stiffness = Chom; - h_s.effective_polarisation = tau_eff; - h_s.mean_strain_localisation_tensors = localisators; + HomogenizationScheme h_s={.homogenized_stiffness = Chom,.effective_polarisation = tau_eff,.mean_strain_localisation_tensors = localisators}; return h_s; }; @@ -114,7 +112,6 @@ namespace tfel::material::homogenization::elasticity { } const auto tau0 = polarisations_[0]; const auto f0 = micro.get_matrix_fraction(); - HomogenizationScheme h_s; auto Chom = C0; auto tau_eff = tau0; std::vector> localisators = {}; @@ -144,9 +141,7 @@ namespace tfel::material::homogenization::elasticity { Chom += fi * (Ci - C0) * Ai; tau_eff += fi * transpose(Ai) * (taui - tau0); } - h_s.homogenized_stiffness = Chom; - h_s.effective_polarisation = tau_eff; - h_s.mean_strain_localisation_tensors = localisators; + HomogenizationScheme h_s={.homogenized_stiffness = Chom,.effective_polarisation = tau_eff,.mean_strain_localisation_tensors = localisators}; return h_s; }; @@ -165,22 +160,21 @@ namespace tfel::material::homogenization::elasticity { const auto np = micro.get_number_of_phases(); const auto f0 = micro.get_matrix_fraction(); const auto C0 = micro.get_matrix_elasticity(); - + const auto polarisations_ = internals::initialize_polarisation(polarisations, np); - - HomogenizationScheme h_s; auto tau_eff = tfel::math::stensor::zero(); - std::vector> localisators = {}; + std::vector> localisators = {}; auto Chom = C0; auto Chom_ = C0; real rel_err = tolerance + 1; + while (rel_err > tolerance) { std::vector> localisators_try = {}; tfel::math::st2tost2 A0 = f0 * tfel::math::st2tost2::Id(); for (std::size_t i = 0; i < np - 1; i++) { - auto phasei = micro.get_inclusionPhase(i); + auto phasei=micro.get_inclusionPhase(i); auto fi = (*phasei).fraction; tfel::math::st2tost2 Ai; if (isotropic) { @@ -228,9 +222,7 @@ namespace tfel::material::homogenization::elasticity { auto Ai = localisators[i + 1]; tau_eff += fi * transpose(Ai) * taui; } - h_s.homogenized_stiffness = Chom; - h_s.effective_polarisation = tau_eff; - h_s.mean_strain_localisation_tensors = localisators; + HomogenizationScheme h_s={.homogenized_stiffness = Chom,.effective_polarisation = tau_eff,.mean_strain_localisation_tensors = localisators}; return h_s; } diff --git a/include/TFEL/Material/PolyCrystalsSlidingSystems.hxx b/include/TFEL/Material/PolyCrystalsSlidingSystems.hxx index 47e6636879..83a9aeb79d 100644 --- a/include/TFEL/Material/PolyCrystalsSlidingSystems.hxx +++ b/include/TFEL/Material/PolyCrystalsSlidingSystems.hxx @@ -35,7 +35,7 @@ namespace tfel::material { typedef NumType real; //! a simple alias typedef tfel::math::stensor<3u> StrainStensor; - //! return the uniq instance of the class + //! return the unique instance of the class static const PolyCrystalsSlidingSystems& getPolyCrystalsSlidingSystems( const std::string& = ""); //! tensor of directional senses, sorted by phases @@ -57,6 +57,34 @@ namespace tfel::material { */ PolyCrystalsSlidingSystems& operator=(const PolyCrystalsSlidingSystems&); }; + + /*! + * \brief the same class as PolyCrystalsSlidingSystems but with + * volume fraction associated to each grain (fourth item of each line of the file + * given as an argument to the method getPolyCrystalsSlidingSystems) + * \param[in] Np : number of phases + * \param[in] GS : class describing the sliding system of one phase + * \param[in] NumType : numeric type used + */ + template + struct ExtendedPolyCrystalsSlidingSystems { + static constexpr unsigned short Nss = GS::Nss; + typedef NumType real; + typedef tfel::math::stensor<3u, real> StrainStensor; + + static const ExtendedPolyCrystalsSlidingSystems &getPolyCrystalsSlidingSystems( + const std::string &); + tfel::math::vector> mus; + + //! \brief volume fractions per phases + tfel::math::vector volume_fractions; + + private: + ExtendedPolyCrystalsSlidingSystems(const std::string &); + ExtendedPolyCrystalsSlidingSystems(const ExtendedPolyCrystalsSlidingSystems &); + ExtendedPolyCrystalsSlidingSystems &operator=( + const ExtendedPolyCrystalsSlidingSystems &); + }; } // end of namespace tfel::material diff --git a/include/TFEL/Material/PolyCrystalsSlidingSystems.ixx b/include/TFEL/Material/PolyCrystalsSlidingSystems.ixx index 695f407e0a..052b05255e 100644 --- a/include/TFEL/Material/PolyCrystalsSlidingSystems.ixx +++ b/include/TFEL/Material/PolyCrystalsSlidingSystems.ixx @@ -87,6 +87,62 @@ namespace tfel::material { } } } + + template + const ExtendedPolyCrystalsSlidingSystems + &ExtendedPolyCrystalsSlidingSystems::getPolyCrystalsSlidingSystems( + const std::string &f) { + static ExtendedPolyCrystalsSlidingSystems gs(f); + return gs; + } + + template + ExtendedPolyCrystalsSlidingSystems:: + ExtendedPolyCrystalsSlidingSystems(const std::string &f) + : mus(Np), volume_fractions(Np) { + using namespace std; + using namespace tfel::math; + using namespace tfel::utilities; + // the sliding systems of one phase + const auto &gs = GS::getSlidingSystems(); + const real pi = real(4) * std::atan(real(1)); + // reading from the specified file + CxxTokenizer file(f); + file.stripComments(); + CxxTokenizer::const_iterator p = file.begin(); + const CxxTokenizer::const_iterator pe = file.end(); + for (unsigned short i = 0; i != Np; i++) { + this->mus[i].resize(Nss); + tmatrix<3u, 3u, real> drot; + const real psi = CxxTokenizer::readDouble(p, pe) * pi / 180.0; + const real the = CxxTokenizer::readDouble(p, pe) * pi / 180.0; + const real phi = CxxTokenizer::readDouble(p, pe) * pi / 180.0; + const real cospsi = std::cos(psi); + const real costhe = std::cos(the); + const real cosphi = std::cos(phi); + const real sinpsi = std::sin(psi); + const real sinthe = std::sin(the); + const real sinphi = std::sin(phi); + drot(0, 0) = cosphi * cospsi - sinphi * costhe * sinpsi; + drot(0, 1) = cosphi * sinpsi + sinphi * costhe * cospsi; + drot(0, 2) = sinphi * sinthe; + drot(1, 0) = -sinphi * cospsi - cosphi * costhe * sinpsi; + drot(1, 1) = -sinphi * sinpsi + cosphi * costhe * cospsi; + drot(1, 2) = cosphi * sinthe; + drot(2, 0) = sinthe * sinpsi; + drot(2, 1) = -sinthe * cospsi; + drot(2, 2) = costhe; + for (unsigned short j = 0; j != Nss; j++) { + auto &mu = this->mus[i][j]; + mu = gs.mus[j]; + mu.changeBasis(drot); + } + this->volume_fractions[i] = CxxTokenizer::readDouble(p, pe); + } + } + + + } // end of namespace tfel::material #endif /* TFEL_MATERIAL_POLYCRYSTALSSLIDINGSYSTEMS_IXX */ diff --git a/tests/Material/IsotropicModuli.cxx b/tests/Material/IsotropicModuli.cxx index 4ecd59b5d0..f4a1475185 100644 --- a/tests/Material/IsotropicModuli.cxx +++ b/tests/Material/IsotropicModuli.cxx @@ -57,6 +57,8 @@ struct IsotropicModuliTest final : public tfel::tests::TestCase { using namespace tfel::material; const auto tens = 3 * kappa * J + 2 * mu * K; const auto yes_it_is_isotropic = isIsotropic(tens, 2 * eps); + + TFEL_TESTS_ASSERT(yes_it_is_isotropic); @@ -64,6 +66,22 @@ struct IsotropicModuliTest final : public tfel::tests::TestCase { const auto Enu = KG.ToYoungNu(); const auto LG = KG.ToLambdaMu(); const auto Enu_bis = LG.ToYoungNu(); + + tfel::material::LambdaMuModuli lg_ter; + tfel::material::LambdaMuModuli lg_quad(LG); + lg_quad= LG; + lg_ter= LG; + + tfel::material::YoungNuModuli Enu_ter; + tfel::material::YoungNuModuli Enu_quad(Enu_bis); + Enu_quad= Enu_bis; + Enu_ter= Enu_bis; + + tfel::material::KGModuli kg_ter; + tfel::material::KGModuli kg_quad(KG); + kg_quad= KG; + kg_ter= KG; + const auto E = Enu.young; const auto nu = Enu.nu;