Geometric variational principles for computational homogenization

TL;DR

This paper introduces geometric variational principles for computational homogenization of elastic composites, enabling analysis under partial or uncertain material properties using FFT-based optimization methods.

Contribution

It develops a novel geometric variational framework for homogenization, incorporating Green's operators and gradient-based optimization, suitable for uncertain material data.

Findings

The proposed functionals effectively measure admissibility discrepancies.

Green's operators facilitate gradient computations in the optimization.

Numerical tests validate the approach against analytical solutions.

Abstract

The homogenization of periodic elastic composites is addressed through the reformulation of the local equations of the mechanical problem in a geometric functional setting. This relies on the definition of Hilbert spaces of kinematically and statically admissible tensor fields, whose orthogonality and duality properties are recalled. These are endowed with specific energetic scalar products that make use of a reference and uniform elasticity tensor. The corresponding strain and stress Green's operators are introduced and interpreted as orthogonal projection operators in the admissibility spaces. In this context and as an alternative to classical minimum energy principles, two geometric variational principles are investigated with the introduction of functionals that aim at measuring the discrepancy of arbitrary test fields to the kinematic, static or material admissibility conditions of the problem. By relaxing the corresponding local equations, this study aims in particular at laying the groundwork for the homogenization of composites whose constitutive properties are only partially known or uncertain. The local fields in the composite and their macroscopic responses are computed through the minimization of the proposed geometric functionals. To do so, their gradients are computed using the Green's operators and gradient-based optimization schemes are discussed. A FFT-based implementation of these schemes is proposed and they are assessed numerically on a canonical example for which analytical solutions are available.

Equations

$$\left\{\begin{aligned} &\text{(i)}&&{\bm{\varepsilon}}({\bm{x}})=\overline{{\bm{\varepsilon}}}+{\bm{\varepsilon}}^{*}({\bm{x}}),\quad{\bm{\varepsilon}}^{*}({\bm{x}})=\frac{1}{2}\big{(}{\bm{\nabla}}{\bm{u}}^{*}({\bm{x}})+{\bm{\nabla}}{{\bm{u}}^{*}}({\bm{x}})^{\!\top}\big{)},\quad{\bm{u}}^{*}\ \text{periodic on}\ \partial\mathcal{V},\\[2.84526pt] &\text{(ii)}&&{\bm{\sigma}}({\bm{x}})={\bm{L}}({\bm{x}}){\bm{\varepsilon}}({\bm{x}}),\\[5.69054pt] &\text{(iii)}&&\operatorname{div}{\bm{\sigma}}({\bm{x}})={\bm{0}},\quad{\bm{\sigma}}\cdot{\bm{n}}\ \text{anti-periodic on}\ \partial\mathcal{V},\end{aligned}\right.$$

$$\frac{1}{2}{{\bm{L}}}_{\operatorname{eff}}\,\overline{{\bm{\varepsilon}}}:\overline{{\bm{\varepsilon}}}=\min_{\begin{subarray}{c}{\bm{e}}\in\mathscrbf{E}\\ \langle{\bm{e}}\rangle=\overline{{\bm{\varepsilon}}}\end{subarray}}\mathcal{J}({\bm{e}})\quad\text{with}\quad\mathcal{J}({\bm{e}})=\frac{1}{2}\big{\langle}{\bm{L}}({\bm{x}}){\bm{e}}({\bm{x}}):{\bm{e}}({\bm{x}})\big{\rangle},$$

$$\big{\langle}f\big{\rangle}=\frac{1}{|\mathcal{V}|}\int_{\mathcal{V}}f({\bm{x}})\,\text{d}{\bm{x}}\quad\text{and}\quad{\bm{s}}:{\bm{e}}=\sum_{i,j=1}^{d}s_{ij}e_{ij}.$$

$$\mathbf{\epsilon }=\arg \underset{\begin{array}{c}\mathbf{e}\in \text{mathscrbf}E\\ ⟨\mathbf{e}⟩=\overline{\mathbf{\epsilon }}\end{array}}{\min }\mathcal{N}(\mathbf{e})\text{with}\mathcal{N}(\mathbf{e})=\text{Updelta}\text{Equil}(\mathbf{L}\mathbf{e})$$

$$\text{Updelta}\text{Equil}(\mathbf{s})=\frac{1}{2}\Vert {\mathbf{P}}_{\text{mathscrbf}{S}^{\perp }}\mathbf{s}{\Vert }^2,$$

$$(\mathbf{\sigma },\mathbf{\epsilon })=\arg \underset{\mathbf{s},\mathbf{e}}{\min }\mathcal{P}(\mathbf{s},\mathbf{e})\text{with}\mathcal{P}(\mathbf{s},\mathbf{e})=\text{Updelta}\text{Compat}(\mathbf{e})+\text{Updelta}\text{Const}(\mathbf{s},\mathbf{e})+\text{Updelta}\text{Equil}(\mathbf{s}).$$

$$\Updelta{\textrm{Compat}}({\bm{e}})=\frac{1}{2}\|\big{\langle}{\bm{e}}\big{\rangle}-\overline{{\bm{\varepsilon}}}+{\bm{P}}_{\!\mathscrbf{E}^{\bot}}{\bm{e}}\|^{2}.$$

$$\big{(}{\bm{\eta}},\tilde{{\bm{\eta}}}\big{)}_{\!\mathscrbf{H}_{e}}=\big{\langle}{\bm{L}}_{0}{\bm{\eta}}({\bm{x}}):\tilde{{\bm{\eta}}}({\bm{x}})\big{\rangle}.$$

$$\mathscrbf{E}_{0}=\left\{{\bm{e}}^{*}\in\mathscrbf{H}_{e}\text{ such that:}\ \exists\,{\bm{w}}^{*}\in{\bm{H}}^{1}_{\text{per}}\!\left(\mathcal{V}\right),\;{\bm{e}}^{*}({\bm{x}})=\frac{1}{2}\big{(}{\bm{\nabla w}}^{*}({\bm{x}})+{\bm{\nabla}}{{\bm{w}}^{*}}({\bm{x}})^{\!\top}\big{)}\text{ in}\ \mathcal{V}\right\}.$$

$$\text{mathscrbf}E=\text{mathscrbf}{U}_e\oplus \text{mathscrbf}{E}_0,$$

$$\text{mathscrbf}{H}_e=\text{mathscrbf}{U}_e\oplus \text{mathscrbf}{E}_0\oplus \text{mathscrbf}{E}^{\perp }.$$

$${\bm{\eta}}=\big{\langle}{\bm{\eta}}\big{\rangle}+{\bm{P}}_{\!\mathscrbf{E}_{0}}{\bm{\eta}}+{\bm{P}}_{\!\mathscrbf{E}^{\bot}}{\bm{\eta}}.$$

$$\big{(}{\bm{\eta}},\tilde{{\bm{\eta}}}\big{)}_{\!\mathscrbf{H}_{e}}=\big{(}{\bm{R}}_{s}{\bm{\eta}},\tilde{{\bm{\eta}}}\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}}\qquad\forall\tilde{{\bm{\eta}}}\in\mathscrbf{H}_{e},$$

$$\big{(}{\bm{\tau}},{\bm{\eta}}\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}}=\big{\langle}{\bm{\tau}}({\bm{x}}):{\bm{\eta}}({\bm{x}})\big{\rangle}.$$

$$\big{(}{\bm{\tau}},\tilde{{\bm{\tau}}}\big{)}_{\!\mathscrbf{H}_{s}}=\big{\langle}{\bm{\tau}}({\bm{x}}):{\bm{L}}_{0}^{-1}\tilde{{\bm{\tau}}}({\bm{x}})\big{\rangle},$$

$$\mathscrbf{S}_{0}=\Big{\{}{\bm{s}}\in\mathscrbf{H}_{s}\text{ such that:}\ \operatorname{div}{\bm{s}}({\bm{x}})={\bm{0}}\text{ in}\ \mathcal{V},\ {\bm{s}}\cdot{\bm{n}}\text{ anti-periodic on}\ \partial V,\ \big{\langle}{\bm{s}}({\bm{x}})\big{\rangle}={\bm{0}}\Big{\}}.$$

$$\text{mathscrbf}S=\text{mathscrbf}{U}_s\oplus \text{mathscrbf}{S}_0,$$

$$\text{mathscrbf}{H}_s=\text{mathscrbf}{U}_s\oplus \text{mathscrbf}{S}_0\oplus \text{mathscrbf}{S}^{\perp }$$

$${\bm{\tau}}=\big{\langle}{\bm{\tau}}\big{\rangle}+{\bm{P}}_{\!\mathscrbf{S}_{0}}{\bm{\tau}}+{\bm{P}}_{\!\mathscrbf{S}^{\bot}}{\bm{\tau}},$$

$$\mathbf{e}\in \text{mathscrbf}E$$

$$\mathbf{s}\in \text{mathscrbf}S$$

$$\text{mathscrbf}{E}_0=(\text{mathscrbf}{S}^{\perp }{)}^{\prime }\text{and}\text{mathscrbf}{S}_0=(\text{mathscrbf}{E}^{\perp }{)}^{\prime }.$$

$$\forall{\bm{s}}\in\mathscrbf{S},\;\forall{\bm{e}}\in\mathscrbf{E}:\qquad\big{(}{\bm{s}},{\bm{e}}\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}}=\big{\langle}{\bm{s}}({\bm{x}}):{\bm{e}}({\bm{x}})\big{\rangle}=\big{\langle}{\bm{s}}({\bm{x}})\big{\rangle}:\big{\langle}{\bm{e}}({\bm{x}})\big{\rangle}.$$

$$\left\{\begin{aligned} &\text{(i)}&&{\bm{\varepsilon}}\in\mathscrbf{E},\quad\big{\langle}{\bm{\varepsilon}}\big{\rangle}=\overline{{\bm{\varepsilon}}},\\[2.84526pt] &\text{(ii)}&&{\bm{\sigma}}({\bm{x}})={\bm{L}}({\bm{x}}){\bm{\varepsilon}}({\bm{x}})\text{ in }\mathcal{V},\\[5.69054pt] &\text{(iii)}&&{\bm{\sigma}}\in\mathscrbf{S}.\end{aligned}\right.$$

$$\text{Find }{\mathbf{e}}^{\ast }\in \text{mathscrbf}{E}_0\text{ such that }\mathbf{s}\stackrel{\mathrm{def}}{=}({\mathbf{L}}_0{\mathbf{e}}^{\ast }-\mathbf{\tau })\in \text{mathscrbf}S.$$

$${\mathbf{\Gamma }}_0:\mathbf{\tau }\mapsto {\mathbf{\Gamma }}_0\mathbf{\tau }={\mathbf{e}}^{\ast }\text{such that }{\mathbf{e}}^{\ast }\text{ is the solution of (refthermoelas)}.$$

$${\mathbf{\Delta }}_0:\mathbf{\eta }\mapsto {\mathbf{\Delta }}_0\mathbf{\eta }={\mathbf{s}}^{\ast }\text{ such that }{\mathbf{s}}^{\ast }\in \text{mathscrbf}{S}_0\text{ and }\mathbf{e}\stackrel{\mathrm{def}}{=}({\mathbf{L}}_0^{-1}{\mathbf{s}}^{\ast }-\mathbf{\eta })\in \text{mathscrbf}E.$$

$$\big{(}{\bm{\eta}},{\bm{\Gamma}}_{0}{\bm{L}}_{0}\tilde{{\bm{\eta}}}\big{)}_{\!\mathscrbf{H}_{e}}=\big{(}{\bm{\Gamma}}_{0}{\bm{L}}_{0}{\bm{\eta}},\tilde{{\bm{\eta}}}\big{)}_{\!\mathscrbf{H}_{e}},\qquad\forall{\bm{\eta}},\tilde{{\bm{\eta}}}\in\mathscrbf{H}_{e}.$$

$$\big{(}{\bm{\eta}},{\bm{\Gamma}}_{0}{\bm{L}}_{0}\tilde{{\bm{\eta}}}\big{)}_{\!\mathscrbf{H}_{e}}=\big{(}{\bm{L}}_{0}{\bm{\eta}},{\bm{\Gamma}}_{0}{\bm{L}}_{0}\tilde{{\bm{\eta}}}\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}}=\big{(}{\bm{\tau}},{\bm{\Gamma}}_{0}\tilde{{\bm{\tau}}}\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}},$$

$${\mathbf{P}}_{\text{mathscrbf}{E}_0}={\mathbf{\Gamma }}_0{\mathbf{L}}_0\text{and}{\mathbf{P}}_{\text{mathscrbf}{E}^{\perp }}={\mathbf{L}}_0^{-1}{\mathbf{\Delta }}_0.$$

$${\mathbf{P}}_{\text{mathscrbf}{S}_0}={\mathbf{\Delta }}_0{\mathbf{L}}_0^{-1}\text{and}{\mathbf{P}}_{\text{mathscrbf}{S}^{\perp }}={\mathbf{L}}_0{\mathbf{\Gamma }}_0.$$

Full text

Keywords:

Composite materials – Helmholtz decomposition – Green’s operators – Lippmann-Schwinger equation – Gradient-based algorithms

1 Introduction

1.1 Context

Focusing on linear composite materials, the early works [10] and [30, 31] have shown that the local fields satisfying the governing equations of the associated mechanical problem are also solutions of some linear integral equations, which are reminiscent of the well-known Lippmann-Schwinger equation. These formulations rely typically on the introduction of a homogeneous comparison material and of the corresponding Green’s operators for the strain or the stress fields. Remarkably, the Fourier transforms of the kernels of these integral operators are known in closed forms for different types of anisotropy of the reference medium, see [8, 22]. These bases have enabled the development of methods aiming at computing local fields, and their macroscopic responses as well, as the solutions to these integral equations, starting from the work [18, 19] where they are solved using a fixed-point iterative scheme in a FFT-based implementation. This method has developed in the field of computational homogenization with successful applications to a wide range of configurations and a concomitant improvement of the corresponding algorithms over the years, see [20] and the references therein. Some of these algorithms were developed without reference to variational principles [18, 21, 5, 15], while others [13, 4, 7] made explicit use of variational properties of the local fields.

A first link between such algorithms has been investigated in [4]. However, it is rather recently, see [7], that the link between the Lippmann-Schwinger equation and the gradient of the strain-based energy functional has been evidenced. In particular, it has been shown that the iterative scheme introduced in [19] can be interpreted as a gradient descent method with fixed step for this functional. The critical observation was that the gradient of this functional can be computed using the Green’s operator when the space of second-order tensor fields is endowed with an energetic scalar product defined by the reference medium. A similar change of metric was previously used in [13]. Doing so, it was then clear that the scheme of [19] and its variations can be obtained directly in the form of gradient-based algorithms according to minimum energy principles. Moreover, the avenues for improvement from the original scheme with fixed step became clear, namely by using optimal or fast gradients methods as done subsequently in [26]. Conjugate-gradient methods have also been investigated in a number of earlier studies, see [32, 4, 6]. However, the latter do not make use of the key property that the integral operator featured in the Lippmann-Schwinger equation is the gradient of the energy functional in a well-chosen Hilbert space and as such it is a self-adjoint operator. This fact has major implications for gradient descent methods, which results in subtle but fundamental differences between these algorithms.

The present study is structured around two key points:

Within the framework of classical energetic variational principles the gradients of given functionals can be computed using the available Green’s operators and efficient gradient-based minimization algorithms can be employed for fast and accurate computations of composites responses. 2. 2.

The computation of these gradients relies on endowing the space of second-order tensor fields with a geometric, i.e., Hilbertian, structure. Such a functional framework is relatively well-known, see [14], and traces back to [16, 28, 29]. What is crucial in the present study is the definition of well-chosen energetic scalar products so that the spaces of kinematically and statically tensors fields, together with their complementary subspaces, are linked by a number of orthogonality or duality properties, with the associated orthogonal projection operators being generated by the Green’s operators.

In this context, this study aims at blending these ideas together by formulating some geometrical variational principles that allow us to address the computational homogenization of composites from a new angle. Our objectives are detailed below, following a preliminary subsection to present the key elements of the problem.

1.2 Preliminaries

Consider a periodic composite material characterized by the unit-cell $\mathcal{V}\subset\mathbb{R}^{d}$ and the fourth-order elasticity tensor ${\bm{L}}({\bm{x}})$ with major and minor symmetries. The local strain and stress fields ${\bm{\varepsilon}}$ and ${\bm{\sigma}}$ solve the so-called local problem in $\mathcal{V}$ consisting of the compatibility equations, constitutive relations and equilibrium equations with periodic boundary conditions:

[TABLE]

with ${\bm{u}}^{*}$ being the fluctuation of the displacement field in $\mathcal{V}$, ${\bm{n}}$ the unit outward normal on $\partial\mathcal{V}$ and $\overline{{\bm{\varepsilon}}}$ an applied macroscopic strain. A mathematical definition of periodicity conditions is given in Appendix A.1.

The average strain $\overline{{\bm{\varepsilon}}}$ being prescribed then the effective elasticity tensor ${{\bm{L}}}_{\operatorname{eff}}$ is characterized by the energetic variational principle:

[TABLE]

where $\mathscrbf{E}$ denotes the space of second-order tensor fields that are admissible strains, which will be properly defined in Section 2.1, and with the averaging operator \big{\langle}\!\cdot\!\big{\rangle} over $\mathcal{V}$ and the doubly contracted product given by

[TABLE]

The actual strain field ${\bm{\varepsilon}}$ solution of the local problem (1) is the unique minimizer in (2). Note that the tensor ${\bm{L}}({\bm{x}})$ in (1.ii) is interpreted as a local operator hence the omission in such a relation of the doubly contracted product. We reserve the latter for products between second-order tensors.

1.3 Objective

In this article, our objective is three-fold:

First, we aim at revisiting from a geometric standpoint the classical energetic variational approaches such as (2), by relying on the geometric properties of the spaces of compatible strains and equilibrated stresses. In the minimum energy principle (2) the strain compatibility equation (1.i) is satisfied through the choice of the minimization space $\mathscrbf{E}$ and the constitutive relations (1.ii) are directly imposed in the definition of the cost functional $\mathcal{J}$. The equilibrium equation (1.iii) for the stress field, formally rewritten as the condition ${\bm{\sigma}}\in\mathscrbf{S}$ for the time being, with a proper definition of the space $\mathscrbf{S}$ given in the next section, is achieved through the minimization of $\mathcal{J}$ over $\mathscrbf{E}$. In this context, we show in Section 3.1 that we can adopt an alternative approach using a variational principle of the form:

[TABLE]

where “$\Updelta{\textrm{Equil}}$” stands for a measure of the static admissibility of the test field defined as ${\bm{s}}={\bm{L}}{\bm{e}}$. With $\mathscrbf{S}$ being the space of admissible stresses, achieving the condition ${\bm{s}}\in\mathscrbf{S}$ is equivalent to minimizing the norm of the projection of ${\bm{s}}$ onto the subspace orthogonal to $\mathscrbf{S}$. Upon introducing the corresponding orthogonal projection operator ${\bm{P}}_{\!\mathscrbf{S}^{\bot}}$ and a suitable norm, then one defines $\mathcal{N}$ with

[TABLE]

which justifies referring to (3) as a geometric variational principle. 2. 2.

Building on this idea, we introduce next a two-field geometric variational principle that allows to treat the constitutive relations (1.ii) on the same level as the strain and the stress admissibility equations. This is of particular interest in the situations where the constitutive model is partially or fully unknown as when dealing with inverse problems of material identification, see [2], in data-driven computational approaches [9], or when material uncertainties must be accounted for, see [23, 27] and the references therein. As shown in Section 3.2, a two-field variational approach can be adopted through an unconstrained minimization problem of the form

[TABLE]

In (4), the functional $\mathcal{P}$ features the term $\Updelta{\textrm{Const}}({\bm{s}},{\bm{e}})$ that is to be defined as a positive measure of the local error in constitutive relations between the test fields ${\bm{s}}$ and ${\bm{e}}$. The other two terms are geometric measures of the kinematic and static admissibilities of these fields: the third term is defined as in (3) while the first one is a measure of the strain compatibility condition, i.e., ${\bm{e}}\in\mathscrbf{E}$ with \big{\langle}{\bm{e}}\big{\rangle}=\overline{{\bm{\varepsilon}}}. By resorting to the projector ${\bm{P}}_{\!\mathscrbf{E}^{\bot}}$ onto the subspace orthogonal to the space $\mathscrbf{E}$ of admissible strains and using a appropriate norm, then $\mathcal{P}$ is defined with

[TABLE]

According to the definition of the functional $\mathcal{P}$, any pair of fields $({\bm{s}},{\bm{e}})$ can be represented as a point in a projection space where its coordinates correspond to the values of $\Updelta{\textrm{Compat}}({\bm{e}})$, $\Updelta{\textrm{Const}}({\bm{s}},{\bm{e}})$ and $\Updelta{\textrm{Equil}}({\bm{s}})$, see the schematics of Fig. 1. This allows in particular to visualize the evolution of the output of an iterative algorithm aiming at computing the solution to the minimization problem (4) and to compare different computation strategies. 3. 3.

The solutions to the geometric variational principles (3) and (4) can be computed using gradient-based minimization algorithms. In this perspective, we show that the gradients of the functionals $\mathcal{N}$ and $\mathcal{P}$ can be expressed by means of the periodic strain and stress Green’s operators, which will be defined in Section 2.2. In fact, those operators are shown to generate orthogonal projectors on the spaces of strain and stress tensor fields when the latter are endowed with well-chosen energetic scalar products. The obtained explicit forms of the gradients of the functionals $\mathcal{N}$ and $\mathcal{P}$ will also shed light on the proposed geometric variational principles by establishing some relationships with the classical minimum energy principles.

The article is organized as follows. In Section 2 the geometric functional framework is set, the Green’s operators are introduced and their properties as projection operators are investigated. Section 3 focuses on the introduction and the study of the geometric variational principles. Numerical implementation and examples are discussed in Section 4. Mathematical definitions and classical properties of the Green’s operators are deferred to the appendices.

2 Strain and stress orthogonal decompositions

2.1 Geometric setting

Let ${\bm{L}}_{0}$ denote a uniform (no spatial dependence) elastic stiffness tensor, i.e., a positive definite fourth-order tensor with major and minor symmetries. One defines as $\mathscrbf{H}_{e}$ the space of symmetric $\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$-tensor fields, with the notation of Appendix A.1, when equipped with the following energetic scalar product

[TABLE]

Doing so, $\mathscrbf{H}_{e}$ is a Hilbert space of strain tensor fields. Let $\mathscrbf{U}_{e}\subset\mathscrbf{H}_{e}$ be the linear subspace of uniform strain fields (typically averages of local fields) and $\mathscrbf{E}_{0}\subset\mathscrbf{H}_{e}$ the linear subspace of kinematically compatible strain fields, which derive from a periodic displacement field:

[TABLE]

By definition, every field ${\bm{e}}^{*}$ in $\mathscrbf{E}_{0}$ satisfies \big{\langle}{\bm{e}}^{*}({\bm{x}})\big{\rangle}={\bm{0}}. As $\mathscrbf{U}_{e}$ and $\mathscrbf{E}_{0}$ are closed subspaces of $\mathscrbf{H}_{e}$, they are Hilbert spaces for the scalar product (5). Moreover, defining the subspace $\mathscrbf{E}$ as

[TABLE]

then $\mathscrbf{H}_{e}$ admits the orthogonal decomposition

[TABLE]

Therefore, any field ${\bm{\eta}}\in\mathscrbf{H}_{e}$ can be decomposed as

[TABLE]

with \big{\langle}\!\cdot\!\big{\rangle}, ${\bm{P}}_{\!\mathscrbf{E}_{0}}$ and ${\bm{P}}_{\!\mathscrbf{E}^{\bot}}$ being the orthogonal projection operators onto $\mathscrbf{U}_{e}$, $\mathscrbf{E}_{0}$ and $\mathscrbf{E}^{\bot}$ respectively, which are mutually orthogonal for the energetic scalar product (5).

Its topological dual is denoted as $\mathscrbf{H}_{s}=\mathscrbf{H}^{\prime}_{e}$. According to the Riesz representation theorem, there exists an isomorphic mapping ${\bm{R}}_{s}:\mathscrbf{H}_{e}\to\mathscrbf{H}_{s}$ such that for every ${\bm{\eta}}\in\mathscrbf{H}_{e}$:

[TABLE]

where \big{(}\cdot,\cdot\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}} denotes the duality product between $\mathscrbf{H}_{s}$ and $\mathscrbf{H}_{e}$. The definition of the energetic scalar product (5) allows to identify any element in $\mathscrbf{H}_{s}$ to a $\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$-tensor field that is dimensionally consistent with a stress field. Moreover, the Riesz mapping reads ${\bm{R}}_{s}:{\bm{\eta}}\mapsto{\bm{L}}_{0}{\bm{\eta}}$ when the duality product is defined according to the principle of virtual work as

[TABLE]

Therefore, $\mathscrbf{H}_{s}$ is itself a Hilbert space with the induced energetic scalar product

[TABLE]

and, identifying $\mathscrbf{H}_{e}$ with its bidual, the inverse operator ${\bm{R}}_{e}={\bm{R}}_{s}^{-1}:\mathscrbf{H}_{s}\to\mathscrbf{H}_{e}$ with ${\bm{R}}_{e}{\bm{\tau}}={\bm{L}}_{0}^{-1}{\bm{\tau}}$ corresponds to the Riesz operator on $\mathscrbf{H}_{s}$.

As in (7), let $\mathscrbf{U}_{s}\subset\mathscrbf{H}_{s}$ be the linear subspace of uniform stress fields and $\mathscrbf{S}_{0}$ the linear subspace of self-equilibrated fields, i.e., divergence-free and mean-free, so that

[TABLE]

The subspaces $\mathscrbf{U}_{s}$ and $\mathscrbf{S}_{0}$ are closed in $\mathscrbf{H}_{s}$ so that, defining the subspace $\mathscrbf{S}$ as

[TABLE]

one has the following orthogonal decomposition

[TABLE]

relatively to the energetic scalar product (10). Thus, any field ${\bm{\tau}}\in\mathscrbf{H}_{s}$ can be decomposed as

[TABLE]

with \big{\langle}\!\cdot\!\big{\rangle}, ${\bm{P}}_{\!\mathscrbf{S}_{0}}$ and ${\bm{P}}_{\!\mathscrbf{S}^{\bot}}$ being three mutually orthogonal projectors onto $\mathscrbf{U}_{s}$, $\mathscrbf{S}_{0}$ and $\mathscrbf{S}^{\bot}$. In this setting, the following lemma holds:

Lemma 1** (Hill-type orthogonality properties).**

[TABLE]

The proof of this lemma can be found in [28] and makes use of the characterization of distributions which are compatible strain fields in the sense of distributions [25, 16]. It establishes that the space $\mathscrbf{E}$ (resp. $\mathscrbf{S}$) is the polar space of $\mathscrbf{S}_{0}$ (resp. $\mathscrbf{E}_{0}$). This implies, see [1], the following duality characterizations of these closed spaces:

[TABLE]

Similar orthogonality properties that make use of the standard $L^{2}$-scalar product on the tensor space $\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$ can be found in particular in [29] and [14]. Lastly, Lemma 1 yields the original result of Hill transposed to periodic boundary conditions, i.e.,

[TABLE]

With these definitions at hand, the local problem (1) can be rewritten in the condensed form:

[TABLE]

Remark 1**.**

From now on, we reserve the notation ${\bm{\varepsilon}},\,{\bm{\sigma}}$ for the actual solution to (15). Tensor fields that belong to the subspaces $\mathscrbf{E}$ and $\mathscrbf{S}$ are referred to as admissible strain and stress, and they are denoted as ${\bm{e}},\,{\bm{s}}$ respectively. The notations ${\bm{\eta}},\,{\bm{\tau}}$ are used for arbitrary fields in $\mathscrbf{H}_{e}$ and $\mathscrbf{H}_{s}$.

2.2 Green’s operators and orthogonal projectors

Two (periodic) Green’s operators, ${\bm{\Gamma}}_{0}$ for the strain field and ${\bm{\Delta}}_{0}$ for the stress field, can be associated with ${\bm{L}}_{0}$ and ${\bm{L}}_{0}^{-1}$. More specifically, for a given field ${\bm{\tau}}$ in $\mathscrbf{H}_{s}$, consider the following Eshelby problem:

[TABLE]

The problem (16) has a unique solution ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$. This allows to define the Green’s operators as follows:

Definition 1**.**

The periodic strain Green’s operator ${\bm{\Gamma}}_{0}:\mathscrbf{H}_{s}\to\mathscrbf{H}_{e}$ of the reference medium with stiffness ${\bm{L}}_{0}$ is defined as

[TABLE]

The periodic stress Green’s operator ${\bm{\Delta}}_{0}:\mathscrbf{H}_{e}\to\mathscrbf{H}_{s}$ is defined in a similar way as

[TABLE]

Classical properties of ${\bm{\Gamma}}_{0}$ (and similarly of ${\bm{\Delta}}_{0}$) are summarized in Appendix B. As discussed in the previous section, the strain and the stress have different dimensions and live in the dual Hilbert spaces $\mathscrbf{H}_{e}$ and $\mathscrbf{H}_{s}$ respectively, each being endowed with its own energetic scalar product (5) or (10). It will now be seen that these energetic scalar products generate two Helmholtz decompositions that are associated with (7) and (12). Such decompositions will be expressed in terms of the operators ${\bm{\Gamma}}_{0}{\bm{L}}_{0}$ and ${\bm{\Delta}}_{0}{\bm{L}}_{0}^{-1}$ so that our aim is now to study some useful properties of the latter.

Lemma 2**.**

Considering the Hilbert space $\mathscrbf{H}_{e}$ endowed with the energetic scalar product (5), then the operator ${\bm{\Gamma}}_{0}{\bm{L}}_{0}:\mathscrbf{H}_{e}\to\mathscrbf{H}_{e}$ is the orthogonal strain projector onto $\mathscrbf{E}_{0}$ as

${\bm{\Gamma}}_{0}{\bm{L}}_{0}$* is idempotent, i.e., ${\bm{\Gamma}}_{0}{\bm{L}}_{0}{\bm{\Gamma}}_{0}{\bm{L}}_{0}={\bm{\Gamma}}_{0}{\bm{L}}_{0}$.* 2. 2.

${\bm{\Gamma}}_{0}{\bm{L}}_{0}$* is self-adjoint, i.e.,*

[TABLE] 3. 3.

For all ${\bm{\eta}}\in\mathscrbf{H}_{e}$ then ${\bm{\Gamma}}_{0}{\bm{L}}_{0}{\bm{\eta}}\in\mathscrbf{E}_{0}$ and for all ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$ it holds ${\bm{\Gamma}}_{0}{\bm{L}}_{0}{\bm{e}}^{*}={\bm{e}}^{*}$.

Proof.

1. The idempotence of ${\bm{\Gamma}}_{0}{\bm{L}}_{0}$ is a direct consequence of (67) in Appendix B.

2. A straightforward calculation shows that

[TABLE]

where ${\bm{\tau}}={\bm{L}}_{0}{\bm{\eta}}$ and $\tilde{{\bm{\tau}}}={\bm{L}}_{0}\tilde{{\bm{\eta}}}$. Then (19) follows from the reciprocity identity satisfied by ${\bm{\Gamma}}_{0}$, see Lemma 3 in Appendix B.

3. The third item is a direct consequence of the definition of ${\bm{\Gamma}}_{0}$ and of Property 2 in [13]. This allows to conclude, see [3], that ${\bm{\Gamma}}_{0}{\bm{L}}_{0}$ is a projection operator from $\mathscrbf{H}_{e}$ onto the subspace $\mathscrbf{E}_{0}$, which is orthogonal for the energetic scalar product (5). ∎

A similar lemma can be proved for the operator ${\bm{\Delta}}_{0}{\bm{L}}_{0}^{-1}:\mathscrbf{H}_{s}\to\mathscrbf{H}_{s}$ using the duality principle of [14] and provided that $\mathscrbf{H}_{s}$ is endowed with the energetic scalar product (10). In particular, ${\bm{\Delta}}_{0}{\bm{L}}_{0}^{-1}$ is the orthogonal stress projector from $\mathscrbf{H}_{s}$ onto $\mathscrbf{S}_{0}$ for the energetic scalar product (10). As a consequence of Lemma 2, one arrives at the main result of this section:

Proposition 1**.**

Considering the orthogonal decomposition (8) such that {\bm{\eta}}=\big{\langle}{\bm{\eta}}\big{\rangle}+{\bm{P}}_{\!\mathscrbf{E}_{0}}{\bm{\eta}}+{\bm{P}}_{\!\mathscrbf{E}^{\bot}}{\bm{\eta}} for all ${\bm{\eta}}\in\mathscrbf{H}_{e}$, then the featured orthogonal projection operators can be expressed in terms of the Green’s operators as

[TABLE]

Similarly, the orthogonal projectors that enter the decomposition (13), i.e., {\bm{\tau}}=\big{\langle}{\bm{\tau}}\big{\rangle}+{\bm{P}}_{\!\mathscrbf{S}_{0}}{\bm{\tau}}+{\bm{P}}_{\!\mathscrbf{S}^{\bot}}{\bm{\tau}} for all ${\bm{\tau}}\in\mathscrbf{H}_{s}$, are given by

[TABLE]

Proof.

For all ${\bm{\eta}}\in\mathscrbf{H}_{e}$, Definition 1 entails that there exists ${\bm{e}}\in\mathscrbf{E}$ such that ${\bm{L}}_{0}{\bm{e}}={\bm{\Delta}}_{0}{\bm{\eta}}-{\bm{L}}_{0}{\bm{\eta}}$. Using (17) with ${\bm{\tau}}={\bm{L}}_{0}{\bm{e}}$ then, on the one hand, there exists ${\bm{s}}\in\mathscrbf{S}$ such that

[TABLE]

using the relation ${\bm{\Gamma}}_{0}{\bm{\Delta}}_{0}={\bm{0}}$, consequence of Lemma 3 in Appendix B. On the other hand, introducing the decomposition {\bm{e}}=\big{\langle}{\bm{e}}\big{\rangle}+{\bm{e}}^{*} with \big{\langle}{\bm{e}}\big{\rangle}\in\mathscrbf{U}_{e} and ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$, then one has

[TABLE]

Therefore, one obtains that

[TABLE]

Finally, according to (18) one has $\langle{\bm{e}}\rangle=-\langle{\bm{\eta}}\rangle$ so that ${\bm{s}}={\bm{L}}_{0}\langle{\bm{\eta}}\rangle$. Combining this last equation with (22) yields

[TABLE]

which after multiplication by ${\bm{L}}_{0}^{-1}$ yields the sought identity (8) by defining the orthogonal projections as \big{\langle}\!\cdot\!\big{\rangle}, ${\bm{P}}_{\!\mathscrbf{E}_{0}}={\bm{\Gamma}}_{0}{\bm{L}}_{0}$ and ${\bm{P}}_{\!\mathscrbf{E}^{\bot}}={\bm{L}}_{0}^{-1}{\bm{\Delta}}_{0}$. The fact that they are mutually orthogonal follows from the definitions of the Green’s operators as well as (66). The identity (13) with (21) is obtained by duality. ∎

Remark 2**.**

The strain and stress Hilbert spaces $\mathscrbf{H}_{e}$ and $\mathscrbf{H}_{s}$ differ by the physical dimension of their elements. Noticeably, one could avoid working with the two different spaces by considering the space of $\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$-tensor fields endowed with the standard $L^{2}$-scalar product. Based on Proposition 1, when (8) is multiplied by ${\bm{L}}_{0}^{1/2}$ or (13) by ${\bm{L}}_{0}^{-1/2}$ then a single orthogonal decomposition is obtained as

[TABLE]

In (23), the operators \big{\langle}\!\cdot\!\big{\rangle}, ${\bm{L}}_{0}^{1/2}{\bm{\Gamma}}_{0}{\bm{L}}_{0}^{1/2}$ and ${\bm{L}}_{0}^{-1/2}{\bm{\Delta}}_{0}{\bm{L}}_{0}^{-1/2}$ are three mutually orthogonal projectors in $\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$. However, according to this decomposition, a field ${\bm{\zeta}}\in\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$ would have neither the dimension of a strain nor of a stress but be consistent with either

[TABLE]

for ${\bm{\eta}}\in\mathscrbf{H}_{e}$ or ${\bm{\tau}}\in\mathscrbf{H}_{s}$. For this particular reason, we prefer to keep working with the spaces $\mathscrbf{H}_{e}$ and $\mathscrbf{H}_{s}$ that are dimensionally consistent with the mechanical problem considered.

3 Geometric variational principles

3.1 Strain-based variational principles

3.1.1 Minimum energy principle

With a slight abuse of notation, the energetic functional in (2) is redefined as $\mathcal{J}:\mathscrbf{E}_{0}\to\mathbb{R}$ with

[TABLE]

with $\overline{{\bm{\varepsilon}}}\in\mathscrbf{U}_{e}$ given. Doing so, the energetic variational principle (2) is equivalent to the minimization of the functional (24) with respect to ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$. As discussed in the introduction, this optimization constraint together with the definition of $\mathcal{J}$ itself enforce both conditions (15.i) and (15.ii). The remaining equilibrium equation (15.iii) for the stress field, i.e., ${\bm{\sigma}}\in\mathscrbf{S}$, is achieved through the minimization of $\mathcal{J}$ over $\mathscrbf{E}_{0}$ as the space $\mathscrbf{S}$ constitutes the polar space of the former. In this context, the gradient of $\mathcal{J}$ at ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$ can be computed. It is defined as the element of $\mathscrbf{E}_{0}$ such that

[TABLE]

where ${\bm{e}}=\overline{{\bm{\varepsilon}}}+{\bm{e}}^{*}$. This entails that \big{(}{\bm{\nabla}}\mathcal{J}({\bm{e}}^{*})-{\bm{L}}_{0}^{-1}{\bm{L}}{\bm{e}}\big{)}\in\mathscrbf{E}_{0}^{\bot}, or equivalently

[TABLE]

As ${\bm{\nabla}}\mathcal{J}({\bm{e}}^{*})\in\mathscrbf{E}_{0}$ and ${\bm{P}}_{\!\mathscrbf{E}_{0}}{\bm{L}}_{0}^{-1}{\bm{L}}{\bm{e}}={\bm{\Gamma}}_{0}{\bm{L}}{\bm{e}}$, it implies that

[TABLE]

Note that the identity (25) was first reported in [7] where the scalar product (5) is recovered through the change of a metric.

The necessary optimality conditions corresponding to the minimization of the functional $\mathcal{J}$ reads

[TABLE]

Using the properties of the Green’s operator ${\bm{\Gamma}}_{0}$ such that for all ${\bm{\eta}}\in\mathscrbf{H}_{e}$ one has

[TABLE]

with ${\bm{\delta L}}({\bm{x}})={\bm{L}}({\bm{x}})-{\bm{L}}_{0}$, and according to the decomposition ${\bm{\varepsilon}}=\overline{{\bm{\varepsilon}}}+{\bm{\varepsilon}}^{*}$, then (26) is equivalent to

[TABLE]

where ${\bm{I}}$ is the identity operator. This equation coincides with the (periodic) Lippmann-Schwinger equation used by [19] in linear composites.

Remark 3**.**

Consider a composite made of non-linear constituents with a local energy density $w({\bm{x}},{\bm{e}})$ which is convex with respect to ${\bm{e}}$ and bounded from above and below by two quadratic functions $\frac{1}{2}L{\bm{e}}:{\bm{e}}$ and $\frac{1}{2}\mu{\bm{e}}:{\bm{e}}$ respectively, such as \big{(}\frac{1}{2}L{\bm{e}}:{\bm{e}}-w({\bm{x}},{\bm{e}})\big{)} and \big{(}w({\bm{x}},{\bm{e}})-\frac{1}{2}\mu{\bm{e}}:{\bm{e}}\big{)} are convex. Then the response of the composite derives from an effective potential ${w}_{\operatorname{eff}}$ defined as

[TABLE]

In this context, it can be proved that

[TABLE]

in $\mathscrbf{E}_{0}$ endowed with the energetic scalar product (5).

3.1.2 Geometric alternative to the energetic principle

Rather than minimizing the energetic functional $\mathcal{J}$ in (24), one can adopt a geometric approach by minimizing the norm of its gradient (25). To do so, consider the geometric functional $\mathcal{N}:\mathscrbf{E}_{0}\to\mathbb{R}$ defined as

[TABLE]

where $\overline{{\bm{\varepsilon}}}\in\mathscrbf{U}_{e}$ is given. According to the Riesz mapping (9), it is seen that $\|{\bm{\eta}}\|_{\mathscrbf{H}_{e}}=\|{\bm{L}}_{0}{\bm{\eta}}\|_{\mathscrbf{H}_{s}}$ for all ${\bm{\eta}}\in\mathscrbf{H}_{e}$. Moreover, Proposition 1 shows that the orthogonal projector from $\mathscrbf{H}_{s}$ onto $\mathscrbf{S}^{\bot}$ for the scalar product (10) is given by ${\bm{P}}_{\!\mathscrbf{S}^{\bot}}={\bm{L}}_{0}{\bm{\Gamma}}_{0}$. Therefore, (28) can be recast as

[TABLE]

with the strain field solution ${\bm{\varepsilon}}^{*}$ being characterized by the variational problem

[TABLE]

The stationary value of $\mathcal{N}$ is zero, which is reached when ${\bm{L}}{\bm{e}}\in\mathscrbf{S}$ with ${\bm{e}}=\overline{{\bm{\varepsilon}}}+{\bm{e}}^{*}$. This constitutes the sought geometric variational principle where the stress admissibility condition is achieved through the minimization of the norm of the projection of the test field defined as ${\bm{\tau}}={\bm{L}}(\overline{{\bm{\varepsilon}}}+{\bm{e}}^{*})$ in the orthogonal space $\mathscrbf{S}^{\bot}$. This justifies the notation used in Section 1.3, i.e., \mathcal{N}({\bm{e}}^{*})=\Updelta{\textrm{Equil}}({\bm{L}}\big{(}\overline{{\bm{\varepsilon}}}+{\bm{e}}^{*})\big{)} with “$\Updelta{\textrm{Equil}}$” being a measure of the static admissibility of ${\bm{\tau}}={\bm{L}}(\overline{{\bm{\varepsilon}}}+{\bm{e}}^{*})$.

Lastly, note that for the strain-based geometric variational principle (29–30), the effective properties ${{\bm{L}}}_{\operatorname{eff}}$ are computed according to the energetic identity (2). In this context, the main result of this section is the following:

Proposition 2**.**

The gradient of $\mathcal{N}$ in $\mathscrbf{E}_{0}$ endowed with the energetic scalar product (5), and at ${\bm{e}}^{*}$, is the element of $\mathscrbf{E}_{0}$ defined as

[TABLE]

where ${\bm{\Gamma}}_{0}$ is the strain Green’s operator associated with the reference elastic medium with stiffness ${\bm{L}}_{0}$.

Proof.

The gradient of $\mathcal{N}$ in $\mathscrbf{E}_{0}$ endowed with the energetic scalar product (5) is defined as the element of $\mathscrbf{E}_{0}$ such that

[TABLE]

where ${\bm{e}}=\overline{{\bm{\varepsilon}}}+{\bm{e}}^{*}$. Making use of the properties of ${\bm{\Gamma}}_{0}$, then straightforward calculations lead to

[TABLE]

The above equation implies that \big{(}{\bm{\nabla}}\mathcal{N}({\bm{e}}^{*})-{\bm{L}}_{0}^{-1}{\bm{L}}{\bm{\Gamma}}_{0}{\bm{L}}{\bm{e}}\big{)}\in\mathscrbf{E}_{0}^{\bot} and therefore

[TABLE]

As ${\bm{\nabla}}\mathcal{N}({\bm{e}}^{*})\in\mathscrbf{E}_{0}$ and ${\bm{P}}_{\!\mathscrbf{E}_{0}}{\bm{L}}_{0}^{-1}{\bm{L}}{\bm{\Gamma}}_{0}{\bm{L}}{\bm{e}}={\bm{\Gamma}}_{0}{\bm{L}}{\bm{\Gamma}}_{0}{\bm{L}}{\bm{e}}$, this leads to

[TABLE]

∎

At this point one has shown that the strain field ${\bm{\varepsilon}}=\overline{{\bm{\varepsilon}}}+{\bm{\varepsilon}}^{*}$ solution of the local problem (1) with imposed macroscopic strain $\overline{{\bm{\varepsilon}}}\in\mathscrbf{U}_{e}$ can be characterized by either the energetic or the geometric variational principles, using the functionals (24) or (29) respectively, which pertain both to the fluctuating term ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$. To shed light on the connection between these two variational principles we consider that ${\bm{\varepsilon}}^{*}$ will be computed by an iterative gradient-based minimization scheme. The following proposition establishes a link between the two variational principles.

Proposition 3**.**

For any ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$, the vector $-{\bm{\nabla}}\mathcal{N}({\bm{e}}^{*})$ is a descent direction for the energetic functional $\mathcal{J}$ and, reciprocally, $-{\bm{\nabla}}\mathcal{J}({\bm{e}}^{*})$ is a descent direction for the geometric functional $\mathcal{N}$.

Proof.

A descent direction for a given functional, say $\mathcal{J}$, at ${\bm{e}}^{*}$ is defined as a tensor ${\bm{p}}\in\mathscrbf{E}_{0}$ such that \big{(}{\bm{p}},{\bm{\nabla}}\mathcal{J}({\bm{e}}^{*})\big{)}_{\!\mathscrbf{H}_{e}}<0. Owing to Proposition 2 and Equation (25), then for any ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$, one has

[TABLE]

with ${\bm{e}}=\overline{{\bm{\varepsilon}}}+{\bm{e}}^{*}$ and where we introduced $\tilde{{\bm{e}}}^{*}={\bm{\Gamma}}_{0}{\bm{L}}{\bm{e}}\in\mathscrbf{E}_{0}$. The properties of ${\bm{\Gamma}}_{0}$ entail

[TABLE]

The positive definiteness of the quadratic form \big{(}{\bm{L}}\tilde{{\bm{e}}}^{*},\tilde{{\bm{e}}}^{*}\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}} allows to conclude that

[TABLE]

∎

This result shows that a gradient descent algorithm based on either ${\bm{\nabla}}\mathcal{J}$ or ${\bm{\nabla}}\mathcal{N}$ would result in the simultaneous minimization of both functionals $\mathcal{J}$ and $\mathcal{N}$.

3.2 A two-field variational principle

As discussed in Section 1.3, there are situations where the constitutive relations are partially, or even fully, unknown. In this context, we now introduce a variational principle that allows leeway in treating these relations through a proper minimization, as it is done for the stress and strain admissibility conditions.

3.2.1 Minimum projections principle

First, the case of linear constituents is considered and the local problem (1) with prescribed macroscopic strain $\overline{{\bm{\varepsilon}}}\in\mathscrbf{U}_{e}$ is addressed in its condensed form (15). On the one hand, to deal with Eqn. (15.ii) in a variational setting we draw from the concept of error in constitutive relations, initially introduced in [11] for error estimation in the finite element method, and define the following functional locally for all ${\bm{x}}\in\mathcal{V}$:

[TABLE]

While expressed in this form this functional can be easily generalized to non-linear composites, see the next section, it can be conveniently rewritten in the linear case as

[TABLE]

which makes clear that $r({\bm{x}},{\bm{\tau}},{\bm{\eta}})\geq 0$ in $\mathcal{V}$, while $r({\bm{x}},{\bm{\tau}},{\bm{\eta}})=0$ locally if and only if ${\bm{\tau}}({\bm{x}})={\bm{L}}({\bm{x}}){\bm{\eta}}({\bm{x}})$. Note that we use the tensor ${\bm{L}}^{-1}$ to define the quadratic form (32), rather than the reference tensor ${\bm{L}}_{0}^{-1}$ that is used in the energetic scalar product (10). Our motivation to do so is to be able to generalize the formulation (32) to the case of non-linear composites.

On the other hand, the stress admissibility condition (15.iii), i.e., ${\bm{\tau}}\in\mathscrbf{S}$, is handled as in the geometric variational principle (29–30) through the minimization of the projection norm $\|{\bm{P}}_{\!\mathscrbf{S}^{\bot}}{\bm{\tau}}\|_{\mathscrbf{H}_{s}}$. Likewise, the strain admissibility condition (15.i), i.e., $({\bm{\eta}}-\overline{{\bm{\varepsilon}}})\in\mathscrbf{E}_{0}$, is achieved by minimizing the norm $\|{\bm{P}}_{\!\mathscrbf{E}_{0}^{\bot}}({\bm{\eta}}-\overline{{\bm{\varepsilon}}})\|_{\mathscrbf{H}_{e}}$ of the projection onto the subspace orthogonal to $\mathscrbf{E}_{0}$. According to the orthogonal decomposition (8), then $\mathscrbf{E}_{0}^{\bot}=\mathscrbf{U}_{e}\oplus\mathscrbf{E}^{\bot}$ and {\bm{P}}_{\!\mathscrbf{E}_{0}^{\bot}}({\bm{\eta}}-\overline{{\bm{\varepsilon}}})=\big{(}{\bm{I}}-{\bm{P}}_{\!\mathscrbf{E}_{0}}\big{)}{\bm{\eta}}-\overline{{\bm{\varepsilon}}}, so that it amounts to the minimization of the functional {\bm{\eta}}\mapsto\|\big{(}{\bm{I}}-{\bm{P}}_{\!\mathscrbf{E}_{0}}\big{)}{\bm{\eta}}-\overline{{\bm{\varepsilon}}}\|_{\mathscrbf{H}_{e}}.

In summary, considering that none of the equations in (15) is enforced exactly leads to the introduction of the cost functional $\mathcal{P}:\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}\to\mathbb{R}$ defined as

[TABLE]

with “$\Updelta{\textrm{Compat}}$”, “$\Updelta{\textrm{Const}}$” and “$\Updelta{\textrm{Equil}}$” being error measures in kinematic, material and static admissibilities respectively. Based on the above, these terms are defined as:

[TABLE]

The stress and strain fields pair $({\bm{\sigma}},{\bm{\varepsilon}})$ solution to (15) are then characterized by the variational problem:

[TABLE]

Note that the stationary value of $\mathcal{P}$ is zero with each of its additive subparts being zero too at the solution. In the perspective of implementing a gradient-based minimization scheme for (35) and as done previously, we now compute the gradient of the functional (33) in $\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}$ endowed with the cross scalar product defined by (10) and (5).

Proposition 4**.**

The partial gradients of $\mathcal{P}$ with respect to ${\bm{\tau}}$ and to ${\bm{\eta}}$ in $\mathscrbf{H}_{s}$ and $\mathscrbf{H}_{e}$ respectively, each being endowed with the associated energetic scalar product (10) and (5) are given by:

[TABLE]

Proof.

The partial gradient ${\bm{\nabla}}_{\!{\bm{\tau}}}\mathcal{P}({\bm{\tau}},{\bm{\eta}})$ of $\mathcal{P}$ with respect to ${\bm{\tau}}\in\mathscrbf{H}_{s}$ is defined as the element of $\mathscrbf{H}_{s}$ that satisfies

[TABLE]

According to (32), (33) and (34) one finds

[TABLE]

Owing to the properties of the Green’s operator ${\bm{\Gamma}}_{0}$, this identity can be recast as

[TABLE]

which proves that

[TABLE]

Similarly, one seeks ${\bm{\nabla}}_{\!{\bm{\eta}}}\mathcal{P}({\bm{\tau}},{\bm{\eta}})$ as the element of $\mathscrbf{H}_{e}$ satisfying:

[TABLE]

By definition, one has

[TABLE]

which can be rewritten as

[TABLE]

an identity which finally yields

[TABLE]

∎

Proposition 4 allows to characterize the solution to the variational problem (33–35). The first-order optimality condition ${\bm{\nabla}}_{\!{\bm{\tau}}}\mathcal{P}({\bm{\sigma}},{\bm{\varepsilon}})={\bm{0}}$ implies that ${\bm{L}}_{0}({\bm{L}}^{-1}{\bm{\sigma}}-{\bm{\varepsilon}})\in\mathscrbf{S}^{\bot}$. Therefore, according to the duality characterization (14) and the Riesz mapping (9), there exists ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$ such that

[TABLE]

Considering the second optimality condition ${\bm{\nabla}}_{\!{\bm{\eta}}}\mathcal{P}({\bm{\sigma}},{\bm{\varepsilon}})={\bm{0}}$ and applying to it the orthogonal projector ${\bm{P}}_{\!\mathscrbf{E}_{0}}$ entail ${\bm{\Gamma}}_{0}({\bm{L}}{\bm{\varepsilon}}-{\bm{\sigma}})={\bm{0}}$. As a consequence, see (66), there exists ${\bm{s}}\in\mathscrbf{S}$ such that

[TABLE]

Combining the equations (36) and (37) leads to ${\bm{L}}{\bm{e}}^{*}=-{\bm{s}}$. Applying the duality product with ${\bm{e}}^{*}$ to the previous equation implies

[TABLE]

where the last equality follows from Lemma 1. Then, the positive definiteness of the quadratic form \big{(}{\bm{L}}{\bm{e}}^{*},{\bm{e}}^{*}\big{)}_{\!\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}} yields ${\bm{e}}^{*}={\bm{0}}$. Hence, from (36) one obtains ${\bm{\sigma}}={\bm{L}}{\bm{\varepsilon}}$ which, inserted back in the optimality conditions, leads to ${\bm{P}}_{\!\mathscrbf{S}^{\bot}}{\bm{\sigma}}={\bm{0}}$ and $({\bm{I}}-{\bm{P}}_{\!\mathscrbf{E}_{0}}){\bm{\varepsilon}}=\overline{{\bm{\varepsilon}}}$. In summary, the optimality conditions associated with $\mathcal{P}$ lead to the following equations

[TABLE]

which are equivalent to the original elasticity problem (15) considered.

Finally, note that for the two-field geometric variational principle (33–35), making use of the energetic principle (2) and based on Hill’s lemma, the effective properties ${{\bm{L}}}_{\operatorname{eff}}$ are computed using the identity:

[TABLE]

3.2.2 Connection with minimum energy principles

Minimum energy principles for admissible fields.

As the two-field variational principle (33–35) has been motivated by geometric considerations, we would like to explore its connections with the conventional minimum energy principles in the case of linear constituents. The relationships between the strain-based principles have been explored in Proposition 3, making use of the properties of the energy functional $\mathcal{J}$ in (24). To do so for the two-field variational principle, we introduce the principle that is dual to (2), i.e., the following stress-based minimum energy principle under controlled overall strain:

[TABLE]

As in Section 3.1.2, the gradient of $\mathcal{J}_{c}$ at ${\bm{\sigma}}\in\mathscrbf{S}$ is defined as the element of $\mathscrbf{S}$ that satisfies

[TABLE]

Therefore \big{(}{\bm{\nabla}}\mathcal{J}_{c}({\bm{s}})-{\bm{L}}_{0}\big{(}{\bm{L}}^{-1}{\bm{s}}-\overline{{\bm{\varepsilon}}}\big{)}\big{)}\in\mathscrbf{S}^{\bot} which, making use of the projector onto $\mathscrbf{S}$, yields:

[TABLE]

Proposition 1 can be used to obtain an equivalent but more convenient form as

[TABLE]

Reminding that, according to the Riesz mapping, one has for all $({\bm{\tau}},{\bm{\eta}})\in\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}$:

[TABLE]

then combining the definitions (33) and (34) with (25), (40) and (41) yields the result below that establishes a link between the geometric functional $\mathcal{P}$ and the energetic ones $\mathcal{J}$ and $\mathcal{J}_{c}$ for tensor fields that are admissible.

Proposition 5**.**

Considering a kinematically admissible strain field ${\bm{e}}\in\mathscrbf{E}$ such that \big{\langle}{\bm{e}}\big{\rangle}=\overline{{\bm{\varepsilon}}} and a statically admissible stress field ${\bm{s}}\in\mathscrbf{S}$, then along the following four “trajectories” it holds:

[TABLE]

Comparison of minimization principles.

Proposition 5 sheds light on the relationships existing between the energy functionals and $\mathcal{P}$ when the latter is employed with an admissible strain and/or an admissible stress field while assuming that the constitutive relations might be satisfied too. As the proposed geometric variational principle (35) does not actually use any such constraint, we explore now its relationships with the minimum energy principles (2) and (39) in general.

For short-hand notations, given $({\bm{\tau}},{\bm{\eta}})\in\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}$ and making use of Proposition 1 we write:

[TABLE]

with the orthogonality in $\mathscrbf{H}_{s}$ and in $\mathscrbf{H}_{e}$ being understood in the sense of the corresponding energetic scalar products. Moreover, for the purpose of establishing the connection between the geometric and the energetic variational principles, let us recognize that the positive definite tensor ${\bm{L}}$ defines some energetic scalar products on the spaces of $\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$ strain and stress tensor fields, that we denote respectively as:

[TABLE]

In this context, the energetic scalar products (5) and (10) are temporarily denoted as

[TABLE]

to emphasize their dependence on the reference tensor ${\bm{L}}_{0}$ for the clarity of the exposition to come.

A tensor field in $\mathscrbf{H}_{e}$ can be seen as a strain field only if it belongs to the subspace $\mathscrbf{E}$, since it thus satisfies the kinematic admissibility condition. In order to extend the energetic functional $\mathcal{J}$ in (24) to any tensor field in $\mathscrbf{H}_{e}$ we evaluate the energy of its admissible component {{\bm{\eta}}}_{\raisebox{-1.5pt}{\scriptsize\mathscrbf{E}_{0}}}\in\mathscrbf{E}_{0}, based on (42), by

[TABLE]

where we make use of the norm associated with the scalar product (43) defined by ${\bm{L}}$. Likewise, for a generic field ${\bm{\tau}}\in\mathscrbf{H}_{s}$, the complementary stress-based mechanical energy under the controlled overall strain $\overline{{\bm{\varepsilon}}}$ is defined in terms of its statically admissible components {{\bm{\tau}}}_{\raisebox{-1.0pt}{\scriptsize\mathscrbf{S}}}\in\mathscrbf{S} as

[TABLE]

using (42) and the norm associated with (43). In this context, by a direct application of Lemma 1, the evaluation of the averaged value of the error in constitutive relations (32) for the admissible fields {{\bm{\tau}}}_{\raisebox{-1.0pt}{\scriptsize\mathscrbf{S}}} and (\overline{{\bm{\varepsilon}}}+{{\bm{\eta}}}_{\raisebox{-1.5pt}{\scriptsize\mathscrbf{E}_{0}}}) yields

[TABLE]

This shows that, for any admissible strain and stress fields {{\bm{\tau}}}_{\raisebox{-1.0pt}{\scriptsize\mathscrbf{S}}}\in\mathscrbf{S} and {{\bm{\eta}}}_{\raisebox{-1.5pt}{\scriptsize\mathscrbf{E}_{0}}}\in\mathscrbf{E}_{0}, the mean error in constitutive relations, which is positive, is equal to the associated total mechanical energy, i.e.,

[TABLE]

Now, our aim is to show that, for any fields $({\bm{\tau}},{\bm{\eta}})\in\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}$, the geometric functional $\mathcal{P}({\bm{\tau}},{\bm{\eta}})$ defined by (33) provides an upper bound on the total energy (45) of their admissible components ({{\bm{\tau}}}_{\raisebox{-1.0pt}{\scriptsize\mathscrbf{S}}},{{\bm{\eta}}}_{\raisebox{-1.5pt}{\scriptsize\mathscrbf{E}_{0}}})\in\mathscrbf{S}\times\mathscrbf{E}_{0}. To do so, let us first rewrite (45) as

[TABLE]

so that

[TABLE]

By Lemma 1 and Definition (43) one has

[TABLE]

Moreover, the triangular inequality and the definition of the error in constitutive relations entail

[TABLE]

Therefore, using (47) and (48) in (46) leads to the following inequality:

[TABLE]

Now, in the definition of the geometric functional (33) with (34) one can recognize that

[TABLE]

using the notation (44) to highlight that the norms are associated with the scalar products defined by the reference tensor ${\bm{L}}_{0}$. Moreover, as the tensor ${\bm{L}}$ is such that its components are real-valued functions in $L^{\infty}\!\left(\mathcal{V}\right)$, the norms defined by ${\bm{L}}$ and ${\bm{L}}_{0}$ are equivalent. Therefore, there exists a constant $C>0$ that depends only on ${\bm{L}}$ and ${\bm{L}}_{0}$, such that

[TABLE]

Using this inequality in (49), together with (45) and the definitions (33–34) of the geometric cost functional yields the following result.

Proposition 6**.**

There exists $c>0$, which depends only on ${\bm{L}}$ and ${\bm{L}}_{0}$, such that for all $({\bm{\tau}},{\bm{\eta}})\in\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}$:

[TABLE]

where ({{\bm{\tau}}}_{\raisebox{-1.0pt}{\scriptsize\mathscrbf{S}}},{{\bm{\eta}}}_{\raisebox{-1.5pt}{\scriptsize\mathscrbf{E}_{0}}})\in\mathscrbf{S}\times\mathscrbf{E}_{0} are the statically and kinematically admissible components of $({\bm{\tau}},{\bm{\eta}})$.

This proposition implies that the overall minimization, with respect to $({\bm{\tau}},{\bm{\eta}})\in\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}$, of the geometric functional $\mathcal{P}$ to its null stationary value results in the overall minimization of the total mechanical energy associated with these fields. Indeed, the pair $({\bm{\sigma}},{\bm{\varepsilon}})$ solution to (15) satisfies \mathcal{J}({\bm{\varepsilon}}-\big{\langle}{\bm{\varepsilon}}\big{\rangle})=-J_{c}({\bm{\sigma}}) with \big{\langle}{\bm{\varepsilon}}\big{\rangle}=\overline{{\bm{\varepsilon}}}, which is equivalent to having $\mathcal{P}({\bm{\sigma}},{\bm{\varepsilon}})=0$.

3.2.3 Non-linear composites

As it can be done for energetic variational principles, we show how to extend the geometric variational principle (33–35) to the case of non-linear constituents. To do so, consider a local energy density $w:\mathcal{V}\times\mathbb{R}^{d\times d}_{\operatorname{sym}}\to\mathbb{R}$, with $\mathbb{R}^{d\times d}_{\operatorname{sym}}$ being the space of symmetric second-order tensors. The potential $w$ is assumed to be convex in $\mathbb{R}^{d\times d}_{\operatorname{sym}}$ and weakly-coercive in the sense that for all ${\bm{\eta}},\tilde{{\bm{\eta}}}\in\mathbb{R}^{d\times d}_{\operatorname{sym}}$

[TABLE]

Its dual $w^{*}:\mathcal{V}\times\mathbb{R}^{d\times d}_{\operatorname{sym}}\to\mathbb{R}$, which is convex in $\mathbb{R}^{d\times d}_{\operatorname{sym}}$, is defined locally according to the classical Legendre-Fenchel transform:

[TABLE]

for all ${\bm{\tau}}\in\mathbb{R}^{d\times d}_{\operatorname{sym}}$. Note that, unlike in the linear case, the generic energy densities $w$ and $w^{*}$ are not necessarily quadratic.

The local error in constitutive relations functional (32) in its form (31) can then be generalized to non-linear constituents as was done in earlier studies, see [12] and the references therein. To do so, one defines:

[TABLE]

As in the linear case, by definition of the Legendre-Fenchel transform one has $r({\bm{x}},{\bm{\tau}},{\bm{\eta}})\geq 0$ in $\mathcal{V}$, while $r({\bm{x}},{\bm{\tau}},{\bm{\eta}})=0$ locally if and only if ${\bm{\tau}}({\bm{x}})=\partial_{{\bm{\eta}}}w({\bm{x}},{\bm{\eta}})$. In this context, one considers the geometric variational principle (33–35) with the term $\Updelta{\textrm{Const}}({\bm{\tau}},{\bm{\eta}})$ being now defined using (52).

The partial gradients of the cost functional (33) that makes use of (52) can then be computed. First, the gradient ${\bm{\nabla}}_{\!{\bm{\tau}}}\mathcal{P}({\bm{\tau}},{\bm{\eta}})$ is defined as the element of $\mathscrbf{H}_{s}$ that satisfies

[TABLE]

As this identity can be rewritten as

[TABLE]

it can be deduced that

[TABLE]

Likewise, the partial gradient ${\bm{\nabla}}_{\!{\bm{\eta}}}\mathcal{P}({\bm{\tau}},{\bm{\eta}})$ is defined as the element of $\mathscrbf{H}_{e}$ satisfying

[TABLE]

from which one obtains finally:

[TABLE]

The knowledge of the gradients (53) and (54) allows the implementation of gradient-based minimization algorithms. Moreover, as in the linear case, consider the associated optimality conditions and let ${\bm{\varepsilon}}$, ${\bm{\sigma}}$ denote the fields that satisfy ${\bm{\nabla}}_{\!{\bm{\tau}}}\mathcal{P}({\bm{\sigma}},{\bm{\varepsilon}})={\bm{0}}$ and ${\bm{\nabla}}_{\!{\bm{\eta}}}\mathcal{P}({\bm{\sigma}},{\bm{\varepsilon}})={\bm{0}}$. Based on the same arguments as those used in Section 3.2.1, these equations imply that there exist ${\bm{e}}^{*}\in\mathscrbf{E}_{0}$ and ${\bm{s}}\in\mathscrbf{S}$ such that

[TABLE]

Inverting the relation $\partial_{{\bm{\tau}}}w^{*}(\cdot,{\bm{\sigma}})={\bm{\varepsilon}}+{\bm{e}}^{*}$ with the help of the Legendre-Fenchel transform (51) leads to ${\bm{\sigma}}=\partial_{{\bm{\eta}}}w(\cdot,{\bm{\varepsilon}}+{\bm{e}}^{*})$, which inserted into the second equation in (55) implies

[TABLE]

On the one hand, applying the duality product with ${\bm{e}}^{*}$ to the above equation and using Lemma 1 entail

[TABLE]

On the other hand, owing to the convexity of $w$ in $\mathscrbf{H}_{e}$ it holds

[TABLE]

Therefore, (56) and (57) imply that

[TABLE]

which finally yields ${\bm{e}}^{*}({\bm{x}})={\bm{0}}$ as a consequence of the assumption (50). In turn, this implies that ${\bm{\sigma}}=\partial_{{\bm{\eta}}}w(\cdot,{\bm{\varepsilon}})$ and thus, from (53) and (54), one gets ${\bm{P}}_{\!\mathscrbf{S}^{\bot}}{\bm{\sigma}}={\bm{0}}$ together with $({\bm{I}}-{\bm{P}}_{\!\mathscrbf{E}_{0}}){\bm{\varepsilon}}=\overline{{\bm{\varepsilon}}}$. Therefore, the optimality conditions on $\mathcal{P}$ yield the equations:

[TABLE]

which are equivalent to the original elasticity problem (15) transposed to non-linear composites.

4 Numerical implementation

4.1 Iterative minimization schemes

4.1.1 Gradient-based algorithms

Throughout this article, a number of variational principles have been investigated, based on energetic or geometric considerations. In this context, the aim of this section is to discuss the numerical implementation of some iterative minimization schemes. As the formulations considered involve either one-field or two-field cost functionals (see (24, 29) and (33) respectively), we consider in a generic setting the following variational problem

[TABLE]

where the cost functional $\mathcal{T}$ is defined in a Hilbert space $\mathscrbf{H}$ that is either $\mathscrbf{H}_{e}$, $\mathscrbf{H}_{s}$ or $\mathscrbf{H}_{s}\times\mathscrbf{H}_{e}$, or a subset thereof, equipped with the corresponding energetic scalar product (5), (10) or their cross product. The search for a minimizer of the variational problem (58) can be performed by several descent algorithms. While our aim is not to optimize or discuss extensively the performances of such algorithms, reference can be made to [32, 4, 6, 7] for conjugate-gradient based implementations and to [26] for a fast gradient method for computational homogeneization. Rather, we focus here on simple implementations and discuss their properties.

A first class of descent methods is obtained by choosing the gradient as the direction of descent:

[TABLE]

One of the simplest descent algorithms along the gradient of $\mathcal{T}$ is obtained with a fixed step $\rho_{n}=-1$.

Remark 4**.**

It is insightful to revisit the original fixed-point algorithm of [18, 19] in the context of gradient-based minimization schemes. Indeed, as noted in [7], the basic scheme introduced in [19] can be interpreted as a gradient descent method with fixed step for the energetic functional $\mathcal{J}$ in (24).

Another descent method is the gradient method with optimal step, which involves the following line search at each step:

[TABLE]

which ensures that $\mathcal{T}$ decreases monotonically over iterations (a property which is not guaranteed by gradient methods with fixed step). In the case of linear constituents, the generic minimization problem (58) is quadratic so that the gradient of $\mathcal{T}$ in $\mathscrbf{H}$ takes the form

[TABLE]

where ${\bm{T}}$ is a linear operator from $\mathscrbf{H}$ into itself and ${\bm{t}}\in\mathscrbf{H}$, which may both be expressed in terms of the reference tensor ${\bm{L}}_{0}$ and the associated Green’s operators ${\bm{\Gamma}}_{0}$ and ${\bm{\Delta}}_{0}$. From (60), the optimal step $\rho_{n}$ can be found analytically at each iteration by writing that

[TABLE]

Note that the above expression could be further simplified given the explicit form of the linear operator ${\bm{T}}$ for specific cost functionals.

Alternatively, conjugate-gradient methods can be used to solve the variational problem (58). For linear materials, the Lippmann-Schwinger integral equation corresponding to the optimality conditions on $\mathcal{T}$ amounts in the linear system ${\bm{T}}\tilde{{\bm{\chi}}}={\bm{t}}$ with the notations of Equation (60). As in the variational framework considered ${\bm{T}}$ is a self-adjoint operator from $\mathscrbf{H}$ into itself for the chosen scalar product, then the conjugate-gradient method may be used to solve this system. It reads as follows:

Initialization: choose ${\bm{\chi}}_{0}$ and compute

[TABLE]

Then: do $n=0,1,\dots$ until convergence

[TABLE]

Although this scheme shares the same structure with those used in [32, 4, 6], there are subtle but fundamental differences between these algorithms. These differences are associated with the use of the energetic scalar products (5) and (10) together with the associated norms.

For non-linear materials, alternative conjugate-gradient methods can be used. The corresponding schemes maintain the previous structure but with ${\bm{r}}_{n}$ being defined as the steepest descent direction at each step, i.e., ${\bm{r}}_{n}=-{\bm{\nabla}}\mathcal{T}({\bm{\chi}}_{n})$, the parameter $\alpha_{n}$ being found by line search as in (59) and the parameter $\beta_{n}$ being chosen according to known formulae, which include the ones by Fletcher-Reeves, Polak-Ribiere, Hestenes-Stiefel and Dai-Yuan.

Remark 5**.**

The gradient-based minimization schemes described previously depend all on the reference tensor ${\bm{L}}_{0}$, through the choice of the scalar product endowing $\mathscrbf{H}$. When assessing the convergence rate of the corresponding algorithms, not all scalar products are equivalent. Therefore, and as is well-known, ${\bm{L}}_{0}$ can be chosen so as to optimize the convergence rate, see the discussion in [20]. An optimization that frees itself from information about the microstructure has been proposed in [19]. A quite natural question is whether it is possible to improve on the corresponding convergence rate when information on the microstructure is actually available and used. This issue is discussed further in [20] but it is beyond the scope of the present paper to investigate this question.

4.1.2 Stopping criteria

In the approach by minimization based on either one of the iterative methods described previously, the expression of the gradient of $\mathcal{T}$ provides a sensible stopping criterion at the iterate $n$ as

[TABLE]

where the tolerance $\delta$ has to be chosen. For illustration purposes, let consider the energetic cost functional $\mathcal{J}:\mathscrbf{E}_{0}\subset\mathscrbf{H}_{e}\to\mathbb{R}$ defined by (24), for which the criterion (61) can be recast as

[TABLE]

On the one hand, this criterion resembles the one suggested by [15] with two differences: (i) the Green’s operator used in the stopping criterion of [15] is slightly different from the actual operator ${\bm{\Gamma}}_{0}$ associated with ${\bm{L}}_{0}$, and (ii) [15] makes use of the $L^{2}$-norm of the operator, whereas (62) involves the norm on $\mathscrbf{H}_{e}$ that is induced by the energetic scalar product (5). On the other hand, the criterion (62) differs from the original stopping criterion of [18, 19], which is based on the divergence of the stress field ${\bm{s}}_{n}$, i.e.,

[TABLE]

where the $L^{2}$-norm of the divergence is computed in Fourier space using Parseval’s theorem. Note that in (62) the parameter $\delta$ can be chosen arbitrarily small, while by contrast $\delta^{\prime}$ in (63) cannot. Moreover, using propositions 2 and 3, a stopping criterion based on the norm $\|{\bm{\nabla}}\mathcal{N}({\bm{e}}^{*}_{n})\|_{\mathscrbf{H}_{e}}$ can also be used.

The stopping criterion (61) implies terminating the iterative minimization scheme once the gradient of the cost functional is sufficiently small. The value of the cost functional $\mathcal{T}$ itself is not relevant in this regard. In fact, for the energetic variational principles such as (2), the stationary values of the corresponding cost functionals are not known beforehand. These values can actually be expressed in terms of the effective tensor ${{\bm{L}}}_{\operatorname{eff}}$ which is rather computed a posteriori. On the contrary, for the geometric variational principles introduced in Section 3, the cost functionals considered are defined so that their stationary values are zero. This allows to use the following stopping criteria:

[TABLE]

where $\delta^{\prime\prime}$ can be chosen arbitrarily small at the continuous level. In the discretized versions of the proposed algorithms, the corresponding criteria must be adapted. Lastly, if one intends to compare the algorithm performances for different choices of the reference tensor ${\bm{L}}_{0}$ then these criteria must be normalized.

4.2 Sample example: the Obnosov problem

4.2.1 FFT-based implementation

The Green’s operator ${\bm{\Gamma}}_{0}$ considered is a volume integral operator defined over the unit cell $\mathcal{V}$. While all previous developments are independent of the boundary conditions provided that they are compatible with Lemma 1, for periodic media it is convenient to express this operator using the Fourier transform $\mathscr{F}$, see Appendix A.2, as

[TABLE]

with the symmetric fourth-order tensor $\hat{{\bm{\Gamma}}}_{0}({\bm{\xi}})$ being defined in closed-form in the Fourier space by

[TABLE]

Owing to the convolution theorem, Eqn. (64) is the Fourier transform of a convolution, hence ${\bm{\Gamma}}_{0}$ is a non-local integral operator in the real space. This is also the case for the stress Green’s tensor ${\bm{\Delta}}_{0}$. Note that, according to Proposition 1, the tensor ${\bm{\Delta}}_{0}$ can be fully expressed in terms of ${\bm{\Gamma}}_{0}$. The Fourier-based formulation (64) is at the foundations of the FFT-based computational homogenization methods for periodic media and it is used in the numerical examples of this section.

In a typical numerical implementation, the unit cell $\mathcal{V}$ is discretized using a regular grid of $N^{d}$ pixels or voxels, which is used to sample data and unknowns. The discrete Fourier transform is computed by means of the FFT algorithm, using all the discrete frequencies associated with the discretization, see [19]. The gradient-based minimization schemes of Section 4.1 are implemented and, in accordance with the general principle of the FFT methods, the algorithms make use of (64) in order to compute the action of the operator ${\bm{\Gamma}}_{0}$ locally in the Fourier space, while actions of operators such as ${\bm{L}}({\bm{x}})$ are computed locally in the real space. At convergence of these iterative schemes the accuracy of the solution is governed by the discretization, i.e., the grid size $N$.

4.2.2 Comparison between minimization schemes

In [24] Y. Obnosov derived the effective conductivity of a square array of square inclusions with volume fraction $0.25$, see Figure 2. As the corresponding problem is equivalent to antiplane elasticity, let ${\bm{L}}({\bm{x}})$ denote the second-order shear modulus matrix. Its effective value is isotropic according to:

[TABLE]

where the $L_{j}$ for $j=1,\,2$ are the isotropic shear moduli of the individual phases, the matrix being phase $2$. The contrast used in the subsequent simulations is $z=10^{2}$ with $L_{2}=1$ and the computational domain is an image of the unit-cell discretized into $512\!\times\!512$ pixels.

Computations were performed using the energetic cost functional (24) with the gradient method with fixed or optimal step and the conjugate-gradient method. Minimizations of the geometric functionals (29) and (33) are also implemented using the gradient method with optimal step and the conjugate-gradient method. For the strain-based schemes associated with $\mathcal{J}$ and $\mathcal{N}$, the corresponding reference medium was the optimal value derived in [13] without information on the microstructure, i.e., ${\bm{L}}_{0}=L_{0}\,{\bm{I}}$ with $L_{0}=\frac{z+1}{2}$. Moreover, they are initialized using ${\bm{\chi}}_{0}\equiv{\bm{e}}^{*}_{0}={\bm{0}}$ with ${\bm{\chi}}_{0}$ being a first-order tensor in the antiplane elasticity problem considered. The effective property is computed at each iterate using {{\bm{L}}}_{\operatorname{eff}}\,\overline{{\bm{\varepsilon}}}:\overline{{\bm{\varepsilon}}}=\big{\langle}{\bm{L}}({\bm{e}}^{*}_{n}+\overline{{\bm{\varepsilon}}}):({\bm{e}}^{*}_{n}+\overline{{\bm{\varepsilon}}})\big{\rangle} according to (2).

For the mixed schemes associated with $\mathcal{P}$ the reference medium was chosen as ${\bm{L}}_{0}=L_{0}\,{\bm{I}}$ with $L_{0}=\sqrt{z}$. Alternative choices might be relevant, see Remark 5, but we do not aim at investigating this issue further in the present study. The schemes are initialized using a pair of non-zero first-order tensors defined arbitrarily as ${\bm{\chi}}_{0}\equiv({\bm{\tau}}_{0},{\bm{\eta}}_{0})=({\bm{1}},\overline{{\bm{\varepsilon}}})$. Note that, the functionals considered being convex, the gradient-based schemes are guaranteed to converge independently of the chosen starting point. Lastly, the effective modulus is computed according to (38), i.e., using {{\bm{L}}}_{\operatorname{eff}}\,\overline{{\bm{\varepsilon}}}:\overline{{\bm{\varepsilon}}}=\big{\langle}{\bm{\tau}}_{n}\big{\rangle}:\big{\langle}{\bm{\eta}}_{n}\big{\rangle}.

For the configuration investigated here, the conjugate-gradient is the fastest method. Moreover, for the minimization of the energetic function $\mathcal{J}$ very few performance differences are found between the gradient-based schemes and the so-called basic scheme of [19]. In terms of computation of the effective property, the minimization of the geometric functionals $\mathcal{N}$ using the conjugate-gradient method yields comparable performances while the optimal-step implementation is associated with a slower convergence rate, see Figure 3a. Similar conclusions are found for the minimization of the geometric functional $\mathcal{P}$ as shown Fig. 3b. All methods compute the effective property with the same precision with a relative error of about $9\cdot 10^{-6}$ measured at convergence. The specific values of the error reached at convergence for such algorithms is specific to (i) the physical configuration considered, (ii) the method implementation and (iii) the discretization used, see the discussion in [20].

For the proposed geometric variational principles, Figure 4 highlights that the functionals $\mathcal{N}$ and $\mathcal{P}$ are minimized monotonically as expected. Moreover, the mechanical energy, as quantified by the energetic functionals $\mathcal{J}$ and $\mathcal{J}_{c}$, for the admissible components of the fields at each iterate decrease in accordance with the propositions 3 and 6 with $\mathcal{J}$ converging to an $\mathcal{O}(1)$ value that is equal to the effective energy, while $(\mathcal{J}+\mathcal{J}_{c})$ converges to zero. In Fig. 4a and 4b it can also be seen that these energy terms are minimized monotonically.

As a conclusion of this section, the comparison of the proposed approach with earlier methods, i.e., the minimization of the geometric functionals compared to the energetic ones, is not characterized by an improvements of the numerical performances, in terms of, e.g., accuracy or speed of convergence. Rather, what should be retained from this comparison is that the proposed method is conceptually different in that it treats the three equations of the local problem on an equal footing. In doing so, its advantage lies in the fact that it make it easy to deal with situations where the constitutive model is partially or fully unknown, see [2, 9, 23, 27].

4.2.3 Visualization of cost functional minimizations

For the Obnosov problem described previously, Figure 5 represents the solutions computed at each iterate for the gradient-based implementations of the proposed geometric variational principle (33–35) using a projection space quantifying the deviation from the validity of each of the equations of the problem (15), see Fig. 1. More precisely, for each iterate the solution is represented as a point with the following coordinates:

[TABLE]

with the functionals “$\Updelta{\textrm{Compat}}$”, “$\Updelta{\textrm{Const}}$” and “$\Updelta{\textrm{Equil}}$” being defined by (34).

For comparison, three classical methods are also implemented and their corresponding solutions represented in the same projection space in Figure 6. The iterates of the so-called basic scheme of [19] are shown in Fig. 6a using the coordinates

[TABLE]

since, by definition, this scheme satisfies by construction both the kinematic admissibility and the constitutive relations. The static admissibility is met only at convergence. The Eyre-Milton scheme of [5] and the augmented-Lagrangian method of [13] are also implemented using the common single-field formalism of [17] and their iterates are represented using the coordinates

[TABLE]

as they only enforce the constitutive relations by construction. Kinematic and static admissibilities are satisfied at convergence. The figures 5 and 6 allow to visualize the behavior of the iterated solution for each algorithm and to compare their convergence performances. For example, Figure 5 sheds a new light on the behavior of the proposed algorithms compared to the only information provided by the computation of the effective property, see Fig. 3b.

5 Conclusions

Considering an arbitrary periodic composite subjected to an applied macroscopic strain, the associated local mechanical problem is recast in a geometric formalism that relies on the definition of spaces of kinematically and statically admissible tensor fields. These are Hilbert spaces associated with two different energetic scalar products that are formulated using a reference and uniform elasticity tensor. The corresponding strain and stress Green’s operators are considered and their geometric properties are investigated. In particular, these operators are shown to generate orthogonal decompositions of second-order tensors fields. In this context, two geometric variational principles are proposed to compute the solution to the mechanical problem considered, as well as the effective properties of the composite. First, a strain-based variational principle is proposed through the introduction of a functional that measures the lack of static admissibility of an associated stress test field. The gradient of this functional is computed and the connection with the classical energetic variational principle is discussed. Building from this preliminary idea, a two-field variational principle is proposed. Its aim is to relax all of the equations of the local mechanical problem so as to reach the solution through an unconstrained minimization process. Doing so, one of our objectives is to enable the treatment of problems where the constitutive relations are partially unknown or uncertain. The functional partial gradients are computed and some connections with minimum energy principles are discussed. With these geometric functionals at hand, their minimization is addressed using gradient-based iterative minimization schemes. A gradient descent scheme with optimal step and the conjugate-gradient method are both implemented and their performances are illustrated on the prototypical Obnosov problem for which analytical solutions are available. Using a FFT-based implementation, the proposed schemes are confronted to the standard approaches that revolve around the classical minimum energy principles. Lastly, the geometric setting considered allows to visualize the evolutions of the iterated solutions in a 3D system of coordinates corresponding to the kinematic, static and material admissibility conditions, which is illustrated for a number of schemes.

This study lays the groundwork for the homogenization of composites whose constitutive properties are partially unknown and this will be the subject of future works. Moreover, there are open questions regarding the connections between the energetic and the geometric variational principles, in particular to assess whether and under which conditions the latter yield some minimization principles for the former. Lastly, in the context considered, let us mention as a perspective that there exist other algorithms pertaining to convex optimization, such as the Alternating Direction Method of Multipliers, which can be used to improve numerical performances.

Acknowledgements:

Fruitful discussions with H. Moulinec and J.-C. Michel are gratefully acknowledged. The Authors have received funding from Excellence Initiative of Aix-Marseille University - A*MIDEX, a French “Investissements d’Avenir” program in the framework of the Labex MEC.

Appendix A Mathematical definitions

A.1 Periodic fields

Consider a unit-cell $\mathcal{V}$ allowing to fill the space $\mathbb{R}^{d}$ by translation along $d$ vectors ${\bm{Y}}_{1},\dots,{\bm{Y}}_{d}$. The lattice $\mathcal{R}$ generated by these vectors is defined as

[TABLE]

Define spaces of periodic scalar functions, vector fields and tensor fields as:

[TABLE]

A.2 Fourier transforms

The Fourier transform $\hat{f}$ of $f$ is defined as:

[TABLE]

Let $\mathcal{R}^{*}$ denote the reciprocal lattice of $\mathcal{R}$ generated by the vectors

[TABLE]

where $(i,j,k)$ is a direct circular permutation. Then, according to Plancherel’s theorem:

[TABLE]

and therefore

[TABLE]

The original periodic function $f$ in $L^{2}_{\text{per}}\!\left(\mathcal{V}\right)$ can be reconstructed from its Fourier transform by

[TABLE]

Appendix B Properties of the Green’s operators

Lemma 3**.**

The strain Green’s operator ${\bm{\Gamma}}_{0}$ satisfies the following properties:

${\bm{\Gamma}}_{0}$* satisfies the reciprocity identity:*

[TABLE] 2. 2.

The kernel of ${\bm{\Gamma}}_{0}$ coincides with the subspace $\mathscrbf{S}$:

[TABLE] 3. 3.

${\bm{\Gamma}}_{0}$* is such that:*

[TABLE]

Properties (66) and (67) were proved in [13]. The additional property (65) derives from the identity

[TABLE]

To prove (68), note that, by definition of ${\bm{\Gamma}}_{0}$, a stress field ${\bm{s}}={\bm{L}}_{0}{\bm{\Gamma}}_{0}{\bm{\tau}}-{\bm{\tau}}$ in $\mathscrbf{S}$ can be associated with ${\bm{\tau}}$ through (17). By Lemma 1:

[TABLE]

which, according to the definition (9) of the Riesz mapping, proves (68). Note that similar properties can be proved for the stress Green’s operator ${\bm{\Delta}}_{0}$ owing to the duality principle, see [14].

Remark 6**.**

If one defines the Green’s operator from $\mathscrbf{L}^{2}_{\text{per}}\!\left(\mathcal{V}\right)$, endowed with the standard $L^{2}$-scalar product, into itself then the reciprocity identity (65) amounts in the self-adjointness of ${\bm{\Gamma}}_{0}$, a property which is known since [10].