The Heterogeneous Multiscale Finite Element Method for the Homogenization of Linear Elastic Solids and a Comparison with the FE$^2$ Method

TL;DR

This paper develops a two-scale finite element method for homogenizing linear elastic solids, extending FE-HMM to vector fields, and compares its performance and convergence with the FE$^2$ method.

Contribution

It introduces an FE-HMM formulation for vector-valued problems in linear elasticity and provides a detailed comparison with FE$^2$, including convergence analysis and superconvergence properties.

Findings

FE-HMM achieves optimal micro-macro mesh refinement strategies.

Superconvergence properties are identified for FE-HMM.

Numerical results show comparable accuracy with FE$^2$.

Abstract

The Heterogeneous Multiscale Finite Element Method (FE-HMM) is a two-scale FEM based on asymptotic homogenization for solving multiscale partial differential equations. It was introduced in [W. E and B. Engquist, \emph{Commun. Math. Sci.}, 1 (2003), 87--132]. The objective of the present work is an FE-HMM formulation for the homogenization of linear elastic solids in a geometrical linear frame, and doing so, for the first time, of a vector-valued field problem. A key ingredient of FE-HMM is that macrostiffness is estimated by stiffness sampling on heterogeneous microdomains in terms a of modified quadrature formula, which implies an equivalence of energy densities of the microscale with the macroscale. Beyond this coincidence with the Hill-Mandel macrohomogeneity condition, which is the cornerstone of the FE$^2$ method, we elaborate a conceptual comparison with the latter method. After developing an algorithmic framework we (i) assess the existing a priori convergence estimates for the micro- and macro-errors in various norms, (ii) verify optimal strategies in uniform micro-macro mesh refinements based on the estimates, (iii) analyze superconvergence properties of FE-HMM, and (iv) compare the numerical results of FE-HMM with those of FE$^2$.

Tables

Table 1. Table 1 : Algorithm for FE-HMM in linear elasticity.

Table 2. Table 2 : Algorithm for assessing superconvergence.

Table 3. Table 3 : Matrix-inclusion problem for beam bending. Comparison of the standard FE-HMM with the FE-HMM( 𝔸 0 , h superscript 𝔸 0 ℎ \mathbb{A}^{0,h} )/FE 2 solution.

	microelements	$64\times 64$
	macroelements	$50\times 10$	$100\times 20$	$200\times 40$	$400\times 80$
FE-HMM	${\Vert 𝒖\Vert }_{\infty }$	11.8018	11.8498	11.8630	11.8667
	${\Vert 𝒖\Vert }_A$	1080.36	1082.50	1083.07	1083.22
FE2 $\widehat{=}$	${\Vert 𝒖\Vert }_{\infty }$	11.8018	11.8498	11.8630	11.8667
FE-HMM(${𝔸}^{0,h}$)	deviation	4.09 $\cdot {10}^{-11}$	-6.72 $\cdot {10}^{-13}$	4.68 $\cdot {10}^{-13}$	3.58 $\cdot {10}^{-11}$
	${\Vert 𝒖\Vert }_A$	1080.36	1082.50	1083.07	1083.22
	deviation	2.18 $\cdot {10}^{-11}$	-3.38 $\cdot {10}^{-11}$	2.35 $\cdot {10}^{-10}$	1.96 $\cdot {10}^{-11}$

Table 4. Table 4 : Superconvergence test. Convergence of the global macro error in the energy-norm, ϵ italic-ϵ \epsilon =5 mm.

macroelements	$25\times 5$	$50\times 10$	$100\times 20$	$200\times 40$	$400\times 80$
microelements	$32\times 32$					order
${\Vert {𝐮}^0-{𝐮}^H\Vert }_A$	4.6020	2.3821	1.3401	0.7034	0.2751	0.9888

Table 5. Table 5 : Microstructure with an analytical solution of homogenization. Convergence of the standard FE-HMM solution in comparison to the FE-HMM( 𝔸 0 , h superscript 𝔸 0 ℎ \mathbb{A}^{0,h} )/FE 2 solution and the percentaged deviation from standard FE-HMM.

	microelements	$20\times 20$
	macroelements	$20\times 20$	$40\times 40$	$80\times 80$	$160\times 160$
FE-HMM	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	77.7902	78.3328	78.5746	78.6865
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	81.6746	81.7999	81.8388	81.8510
FE2 $\widehat{=}$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	77.7902	78.3328	78.5746	78.6865
FE-HMM(${𝔸}^{0,h}$)	deviation	2.82 $\cdot {10}^{-13}$	1.50 $\cdot {10}^{-13}$	2.20 $\cdot {10}^{-13}$	−1.64 $\cdot {10}^{-11}$
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	81.6746	81.7999	81.8388	81.8510
	deviation	1.38 $\cdot {10}^{-13}$	1.16 $\cdot {10}^{-13}$	1.26 $\cdot {10}^{-12}$	−8.78 $\cdot {10}^{-12}$

Table 6. Table 6 : Non-uniformly periodic tensor. Convergence of the standard FE-HMM solution in comparison to the FE 2 solution along with the perecentaged deviation.

	microelements	$20\times 20$
	macroelements	$20\times 20$	$40\times 40$	$80\times 80$	$160\times 160$
FE-HMM	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	20.0246	20.1980	20.2810	20.3214
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	4.15603	4.16176	4.16353	4.16408
FE2 $\widehat{=}$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	20.0246	20.1980	20.2810	20.3214
FE-HMM(${𝔸}^{0,h}$)	deviation	3.79 $\cdot {10}^{-11}$	4.47 $\cdot {10}^{-14}$	1.41 $\cdot {10}^{-14}$	-1.29 $\cdot {10}^{-13}$
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	4.15603	4.16176	4.16353	4.16408
	deviation	2.73 $\cdot {10}^{-14}$	2.81 $\cdot {10}^{-14}$	-8.33 $\cdot {10}^{-16}$	-2.70 $\cdot {10}^{-13}$

Table 7. Table B.7 : Matrix-inclusion problem for beam bending. Macroconvergence of FE-HMM keeping the microdiscretization fixed, variation of ϵ italic-ϵ \epsilon .

	macroelements	$50\times 10$	$100\times 20$	$200\times 40$	$400\times 80$	$800\times 160$
$ϵ$	microelements	$32\times 32$
$5$	${\Vert {𝒖}^H\Vert }_{\infty }$	11.7997	11.8477	11.8609	11.8646	11.8658
	${\Vert {𝒖}^H\Vert }_A$	1080.26	1082.41	1082.98	1083.13	1083.17
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2}$	703.5393	189.9063	52.2047	14.5294	3.9509
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1}$ in ${10}^{-3}$	179.1201	89.1681	48.5004	26.8540	13.9814
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_A$	34.6706	17.3223	9.4171	5.2122	2.7175
$0.5$	${\Vert 𝒖\Vert }_{\infty }$	11.7997	11.8477	11.8609	11.8646	11.8658
	${\Vert 𝒖\Vert }_A$	1080.26	1082.41	1082.98	1083.13	1083.17
$0.05$	${\Vert 𝒖\Vert }_{\infty }$	11.7997	11.8477	11.8609	11.8646	11.8658
	${\Vert 𝒖\Vert }_A$	1080.26	1082.41	1082.98	1083.13	1083.17

Table 8. Table B.8 : Matrix-inclusion problem for beam bending. Microconvergence of FE-HMM keeping the macrodiscretization fixed, variation of ϵ italic-ϵ \epsilon .

	macroelements	$200\times 40$
$ϵ$	microelements	$16\times 16$	$32\times 32$	$64\times 64$	$128\times 128$	$256\times 256$
$5$	${\Vert {𝒖}^H\Vert }_{\infty }$	11.8545	11.8609	11.8630	11.8636	11.8638
	${\Vert {𝒖}^H\Vert }_A$	1082.68	1082.98	1083.07	1083.10	1083.11
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2}$	102.4091	32.6310	10.0525	3.0164	0.8806
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1}$ in ${10}^{-4}$	48.0344	15.2866	4.6909	1.3997	0.4056
$0.5$	${\Vert 𝒖\Vert }_{\infty }$	11.8545	11.8609	11.8630	11.8636	11.8638
	${\Vert 𝒖\Vert }_A$	1082.68	1082.98	1083.07	1083.10	1083.11
$0.05$	${\Vert 𝒖\Vert }_{\infty }$	11.8545	11.8609	11.8630	11.8636	11.8638
	${\Vert 𝒖\Vert }_A$	1082.68	1082.98	1083.07	1083.10	1083.11

Table 9. Table B.9 : Matrix-inclusion problem for beam bending. Convergence of 𝔸 0 , h superscript 𝔸 0 ℎ \mathbb{A}^{0,h} components.

$A_{ij}^{0,h}$	microelements
$ij$	$16\times 16$	$32\times 32$	$64\times 64$	$128\times 128$	$256\times 256$
$11$	46721.57	46700.78	46694.06	46691.95	46691.31
$12$	11662.05	11666.10	11667.41	11667.79	11667.90
$22$	46721.57	46700.78	46694.06	46691.95	46691.31
$33$	17443.96	17436.62	17434.23	17433.48	17433.55

Table 10. Table B.10 : Microstructure with an analytical solution of homogenization. Macroconvergence of FE-HMM keeping the microdiscretization fixed, variation of ϵ italic-ϵ \epsilon .

	macroelements	$20\times 20$	$40\times 40$	$80\times 80$	$160\times 160$	$320\times 320$
$ϵ$	microelements	$40\times 40$
$0.005$	${\Vert {𝒖}^H\Vert }_{\infty }$ in ${10}^{-3}$	77.8920	78.4353	78.6774	78.7895	78.8427
	${\Vert {𝒖}^H\Vert }_A$ in ${10}^{-3}$	81.7244	81.8499	81.8888	81.9011	81.9050
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2}$ in ${10}^{-6}$	47.9037	15.2437	4.9547	1.6312	0.5347
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1}$ in ${10}^{-5}$	38.2124	20.8187	11.7254	6.6562	3.7066
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_A$ in ${10}^{-4}$	38.0001	20.7804	11.7107	6.6520	3.7103
$0.0005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	77.8920	78.4353	78.6774	78.7895	78.8427
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	81.7244	81.8499	81.8888	81.9011	81.9050
$0.00005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	77.8920	78.4353	78.6774	78.7895	78.8427
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	81.7244	81.8499	81.8888	81.9011	81.9050

Table 11. Table B.11 : Microstructure with an analytical solution of homogenization. Microconvergence of FE-HMM keeping the macrodiscretization fixed, variation of ϵ italic-ϵ \epsilon .

	macroelements	$40\times 40$
$ϵ$	microelements	$16\times 16$	$32\times 32$	$64\times 64$	$128\times 128$	$256\times 256$
$0.005$	${\Vert {𝒖}^H\Vert }_{\infty }$ in ${10}^{-3}$	78.3328	78.4353	78.4611	78.4675	78.4691
	${\Vert {𝒖}^H\Vert }_A$ in ${10}^{-3}$	81.7999	81.8499	81.8624	81.8656	81.8664
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2}$ in ${10}^{-7}$	199.8699	50.1244	12.5275	3.1182	0.7653
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1}$ in ${10}^{-7}$	267.2828	67.0309	16.7529	4.1700	1.0234
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_A$ in ${10}^{-7}$	422.0078	105.7584	26.4273	6.5777	1.6142
$0.0005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	78.3328	78.4353	78.4611	78.4675	78.4691
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	81.7999	81.8499	81.8624	81.8656	81.8664
$0.00005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	11.8545	78.4353	78.4611	78.4675	78.4691
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	81.7999	81.8499	81.8624	81.8656	81.8664

Table 12. Table B.12 : Microstructure with an analytical solution of homogenization. Convergence of the coefficients of 𝔸 0 , h superscript 𝔸 0 ℎ \mathbb{A}^{0,h} .

	microelements
$A_{ij}^{0,h}$	$20\times 20$	$40\times 40$	$80\times 80$	$160\times 160$	$320\times 320$
$11$	100.235	100.059	100.015	100.004	100.001
$22$	140.028	140.028	140.028	140.028	140.028

Table 13. Table B.13 : Optimal refinement strategies. Comparison of the L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} errors for the FE-HMM.

	macroelements	$16\times 16$	$64\times 64$	$144\times 144$	$256\times 256$	$576\times 576$
$H^1$	microelements	$4\times 4$	$8\times 8$	$12\times 12$	$16\times 16$	$24\times 24$
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2\left(\Omega \right)}$ in ${10}^{-5}$	46.1045	12.8707	5.6976	3.1917	1.4097
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1\left(\Omega \right)}$ in ${10}^{-5}$	76.6287	22.2750	10.5488	6.2132	2.8595
$L^2$	microelements	$16\times 16$	$64\times 64$	$144\times 144$	$256\times 256$	$576\times 576$
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2\left(\Omega \right)}$ in ${10}^{-7}$	994.5764	89.4812	22.7339	8.5862	1.9925
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1\left(\Omega \right)}$ in ${10}^{-5}$	47.9858	14.1055	7.2641	4.4981	2.1436

Table 14. Table B.14 : Non-uniformly periodic tensor. Macroconvergence of FE-HMM keeping the microdiscretization fixed, variation of ϵ italic-ϵ \epsilon .

	macroelements	$20\times 20$	$40\times 40$	$80\times 80$	$160\times 160$	$320\times 320$
$ϵ$	microelements	$40\times 40$
$0.005$	${\Vert {𝒖}^H\Vert }_{\infty }$ in ${10}^{-3}$	20.0360	20.2095	20.2925	20.3330	20.3530
	${\Vert {𝒖}^H\Vert }_A$ in ${10}^{-3}$	4.1572	4.1630	4.1647	4.1653	4.1654
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2}$ in ${10}^{-7}$	104.7913	32.6314	10.3133	3.2907	1.0449
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1}$ in ${10}^{-6}$	108.7447	58.0419	31.7047	17.4515	9.4527
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_A$ in ${10}^{-5}$	18.5865	9.9356	5.4413	3.0054	1.6343
$0.0005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	20.0360	20.2095	20.2925	20.3330	20.3530
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	4.1572	4.1630	4.1647	4.1653	4.1654
$0.00005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	20.0360	20.2095	20.2925	20.3330	20.3530
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	4.1572	4.1630	4.1647	4.1653	4.1654

Table 15. Table B.15 : Non-uniformly periodic tensor. Microconvergence of FE-HMM keeping the macrodiscretization fixed, variation of ϵ italic-ϵ \epsilon .

	macroelements	$40\times 40$
$ϵ$	microelements	$20\times 20$	$40\times 40$	$80\times 80$	$160\times 160$	$320\times 320$
$0.005$	${\Vert {𝒖}^H\Vert }_{\infty }$ in ${10}^{-3}$	20.1980	20.2095	20.2123	20.2131	20.2133
	${\Vert {𝒖}^H\Vert }_A$ in ${10}^{-3}$	4.16176	4.16295	4.16325	4.16333	4.16335
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L^2}$ in ${10}^{-8}$	235.0888	59.0407	14.7612	3.6745	0.9018
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H^1}$ in ${10}^{-8}$	290.3249	72.9140	18.2298	4.5379	1.1137
	${\Vert {𝒖}^0-{𝒖}^H\Vert }_A$ in ${10}^{-8}$	137.4293	35.1447	8.8671	2.2173	0.5454
$0.0005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	20.1980	20.2095	20.2123	20.2131	20.2133
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	4.16176	4.16295	4.16325	4.16333	4.16335
$0.00005$	${\Vert 𝒖\Vert }_{\infty }$ in ${10}^{-3}$	20.1980	20.2095	20.2123	20.2131	20.2133
	${\Vert 𝒖\Vert }_A$ in ${10}^{-3}$	4.16176	4.16295	4.16325	4.16333	4.16335

Table 16. Table B.16 : Non-uniformly periodic tensor. Convergence of the coefficients of 𝔸 0 , h superscript 𝔸 0 ℎ \mathbb{A}^{0,h} .

$A_{ij}^{0,h}$	microelements
$ij$	$20\times 20$	$40\times 40$	$80\times 80$	$160\times 160$	$320\times 320$
11	3.74326	3.74018	3.73940	3.73921	3.73916
12	0.94034	0.93979	0.93965	0.93961	0.93961
22	3.99944	3.99767	3.99723	3.99712	3.99709
33	1.40440	1.40325	1.40295	1.40288	1.40286

Table 17. Table B.17 : Escher’s periodic Bird and Fish tesselation. Macro-convergence of FE-HMM keeping the microdiscretization fixed.

Macro DOFs	$1000$	$4000$	$16000$	$64000$	$256000$
Micro DOFs / $h$	$1462$ / 1
${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L_2}\cdot {10}^{-2}$	0.2042	0.0554	0.0153	0.0043	0.0012
${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H_1}\cdot {10}^{-6}$	0.5304	0.2665	0.1459	0.0811	0.0423
${\Vert {𝒖}^0-{𝒖}^H\Vert }_A\cdot {10}^{-2}$	5.2007	2.6193	1.4314	0.7948	0.4154

Table 18. Table B.18 : Escher’s periodic Bird and Fish tesselation. Micro-convergence of FE-HMM keeping the macrodiscretization fixed.

Macro DOFs	$16000$
Micro DOFs / $h$	1462/1	2486/0.75	5502/0.5	21670/0.25	82856/0.125
${\Vert {𝒖}^0-{𝒖}^H\Vert }_{L_2}\cdot {10}^{-2}$	0.1438	0.0969	0.0491	0.0150	0.0036
${\Vert {𝒖}^0-{𝒖}^H\Vert }_{H_1}\cdot {10}^{-7}$	0.6395	0.4349	0.2202	0.0677	0.0166
${\Vert {𝒖}^0-{𝒖}^H\Vert }_A\cdot {10}^{-2}$	0.8506	0.5755	0.2925	0.0900	0.0221

Equations

$$\left.\begin{array}[]{rcl}-\dfrac{\partial}{\partial x_{j}}\left(\mathbb{A}^{\epsilon}_{ijlm}\dfrac{\partial u^{\epsilon}_{l}}{\partial x_{m}}\right)&=&f_{i}\qquad\mbox{in}\quad\mathcal{B}\\ u_{i}^{\epsilon}&=&\bar{u}_{i}\qquad\mbox{on}\quad\partial\mathcal{B}_{D}\\ \left(\mathbb{A}^{\epsilon}_{ijlm}\dfrac{\partial u^{\epsilon}_{l}}{\partial x_{m}}\right)\,n_{j}&=&\bar{t}_{i}\qquad\mbox{on}\quad\partial\mathcal{B}_{N}\\ \end{array}\right\}\;$$

$${\epsilon }_{ij}({\mathbf{u}}^{ϵ})$$

$$\mathbf{\epsilon }({\mathbf{u}}^{ϵ})=\mathbf{L}{\mathbf{u}}^{ϵ}.$$

$$\underset{\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}{B}_{ϵ}({\mathbf{u}}^{ϵ},\mathbf{v})}{\underbrace{{\int }_{\mathcal{B}}{\mathbb{A}}^{ϵ}(\mathbf{x})\mathbf{\epsilon }({\mathbf{u}}^{ϵ}):\mathbf{\epsilon }(\mathbf{v})dV}}=\underset{\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}F(\mathbf{v})}{\underbrace{{\int }_{\mathcal{B}}\mathbf{f}\cdot \mathbf{v}dV+{\int }_{\partial {\mathcal{B}}_N}\bar{\mathbf{t}}\cdot \mathbf{v}dA}}$$

$$\mathcal{V}=\{\mathbf{v};\mathbf{v}\in H^1(\mathcal{B}{)}^d,\mathbf{v}{|}_{\partial {\mathcal{B}}_D}=0\}.$$

$$\left.\begin{array}[]{rcl}-\dfrac{\partial}{\partial x_{j}}\left(\mathbb{A}^{0}_{ijlm}\dfrac{\partial u^{0}_{l}}{\partial x_{m}}\right)&=&\langle f_{i}\rangle\qquad\mbox{in}\quad\mathcal{B}\\ u_{i}^{0}&=&\langle\bar{u}_{i}\rangle_{\Gamma}\qquad\mbox{on}\quad\partial\mathcal{B}_{D}\\ \left(\mathbb{A}^{0}_{ijlm}\dfrac{\partial u^{0}_{l}}{\partial x_{m}}\right)\,n_{j}&=&\langle\bar{t}_{i}\rangle_{\Gamma}\qquad\mbox{on}\quad\partial\mathcal{B}_{N}\\ \end{array}\right\}\;$$

$$B_0({\mathbf{u}}^0,\mathbf{v})={\int }_{\mathcal{B}}{\mathbb{A}}^0\mathbf{\epsilon }({\mathbf{u}}^0):\mathbf{\epsilon }(\mathbf{v})dV={\int }_{\mathcal{B}}\mathbf{f}\cdot \mathbf{v}dV+{\int }_{\partial {\mathcal{B}}_N}\bar{\mathbf{t}}\cdot \mathbf{v}dA\forall \mathbf{v}\in \mathcal{V},$$

$${\mathcal{S}}_{\partial {\mathcal{B}}_D}^p(\mathcal{B},{\mathcal{T}}_H)=\left\{{\mathbf{u}}^H\in H^1(\mathcal{B}{)}^d;{\mathbf{u}}^H{|}_{\partial {\mathcal{B}}_D}=\bar{\mathbf{u}};{\mathbf{u}}^H{|}_K\in {\mathcal{P}}^p(K{)}^d,\forall K\in {\mathcal{T}}_H\right\},$$

$$B_H({\mathbf{u}}^H,{\mathbf{v}}^H)={\int }_{\mathcal{B}}\mathbf{f}\cdot {\mathbf{v}}^{H}dV+{\int }_{\partial {\mathcal{B}}_N}\bar{\mathbf{t}}\cdot {\mathbf{v}}^{H}dA\forall {\mathbf{v}}^H\in {\mathcal{S}}_{\partial {\mathcal{B}}_D}(\mathcal{B},{\mathcal{T}}_H),$$

$$B_H({\mathbf{u}}^H,{\mathbf{v}}^H)$$

$$\approx$$

$$\left.\begin{array}[]{rcl}\left(\bm{u}^{h}_{K_{\delta}}-\bm{u}^{H}_{lin,K_{\delta}}\right)&\in&\mathcal{S}^{q}(K_{\delta},\mathcal{T}_{h})\\[5.69054pt] B_{K_{\delta}}(\bm{u}^{h}_{K_{\delta}},\bm{w}^{h}_{K_{\delta}})&:=&\displaystyle{\int_{K_{\delta}}}\mathbb{A}^{\epsilon}(\bm{x})\,\bm{\varepsilon}(\bm{u}^{h}_{K_{\delta}}):\bm{\varepsilon}(\bm{w}^{h}_{K_{\delta}})\,dV=0\\ &&\forall\,\bm{w}^{h}_{K_{\delta}}\in\mathcal{S}^{q}(K_{\delta},\mathcal{T}_{h})\,,\end{array}\quad\right\}\;$$

$${\mathcal{S}}^q(K_{\delta },{\mathcal{T}}_h)=\{{\mathbf{w}}^h\in \mathcal{W}(K_{\delta });{\mathbf{w}}^h\in ({\mathcal{P}}^q(T){)}^d,T\in {\mathcal{T}}_h\}.$$

$${\mathbf{u}}_{lin,K_{\delta }}^H={\mathbf{u}}^H({\mathbf{x}}_{K_l})+(\mathbf{x}-{\mathbf{x}}_{K_l})\cdot \nabla {\mathbf{u}}^H({\mathbf{x}}_{K_l}).$$

$${\int }_{K_{\delta }}{\mathbb{A}}^{ϵ}(\mathbf{x})\mathbf{\epsilon }({\mathbf{\chi }}_{K_{\delta }}^h+{\mathbf{u}}_{lin,K_{\delta }}^H):\mathbf{\epsilon }({\mathbf{w}}_{K_{\delta }}^h)dV=0\forall {\mathbf{w}}_{K_{\delta }}^h\in {\mathcal{S}}^q(K_{\delta },{\mathcal{T}}_h).$$

$${\mathbf{u}}^H$$

$${\mathbf{u}}^h$$

$${\mathbf{k}}_{IJ}^{e,mac}=B_H^e({\mathbf{N}}_I^H,{\mathbf{N}}_J^H)=\sum _{l=1}^{N_{qp}}\frac{{\omega }_{K_l}}{|K_{\delta }|}{\int }_{K_{\delta }}(\mathbf{L}{\mathbf{u}}_{K_{\delta }}^{h(I)}{)}^T{\mathbb{A}}^{ϵ}(\mathbf{x})\mathbf{L}{\mathbf{u}}_{K_{\delta }}^{h(J)}dV.$$

$$\left.\begin{array}[]{rcl}\left(\bm{u}^{h(I,x_{i})}_{K_{\delta}}-{\bm{u}}^{H(I,x_{i})}_{lin,K_{\delta}}\right)&\in&\mathcal{S}^{q}(K_{\delta},\mathcal{T}_{h})\\[5.69054pt] B_{K_{\delta}}(\bm{u}^{h(I,x_{i})}_{K_{\delta}},\bm{w}^{h(I,x_{i})}_{K_{\delta}})&:=&\displaystyle{\int_{K_{\delta}}}(\bm{L}\bm{u}^{h(I,x_{i})}_{K_{\delta}})^{T}\mathbb{A}^{\epsilon}(\bm{x})\,\bm{L}\bm{w}^{h(I,x_{i})}_{K_{\delta}}dV=0\\ &&\forall\,\bm{w}^{h(I,x_{i})}_{K_{\delta}}\in\mathcal{S}^{q}_{per}(K_{\delta},\mathcal{T}_{h})\end{array}\quad\right\}\;.$$

$${\mathbf{u}}_{lin,K_{\delta }}^{H(I,x_i)}$$

$${\mathbf{u}}_{K_{\delta }}^{h(I,x_i)}$$

$${\mathbf{k}}_{IJ}^{e,mac}$$

$${\mathbf{k}}_K^{e,mac}$$

$$\text{mbox}with{\mathbf{T}}_{K_{\delta }}$$

$${\mathbf{k}}_K^{e,mac}$$

$${\mathbf{T}}_{K_{\delta }}$$

$${\mathbf{K}}_{K_{\delta }}^{mic}$$

$${\mathbf{k}}_{mn}^{e,mic}$$

$$\bm{K}^{mic}_{K_{\delta}}=\begin{picture}(8.0,7.0)\put(0.0,-4.0){ $\stackrel{{\scriptstyle\stackrel{{\scriptstyle\mbox{\tiny$num_{ele}$}}}{{\mbox{\Sf A}}}}}{{\mbox{\tiny e = 1}}}$}\end{picture}\;\,\bm{k}_{K_{\delta}}^{e,mic}\,$$

$${\int }_{K_{\delta }}(\mathbf{L}{\mathbf{u}}_{K_{\delta }}^{h(I,x_i)}{)}^T{\mathbb{A}}^{ϵ}\mathbf{L}{\mathbf{u}}_{K_{\delta }}^{h(I,x_i)}dV\to \text{mbox}min$$

$$\mathcal{L}({\mathbf{d}}^{h(I,x_i)},{\mathbf{\lambda }}^{(I,x_i)})=\frac{1}{2}{\left({\mathbf{d}}^{h(I,x_i)}\right)}^T{\mathbf{K}}_{K_{\delta }}^{mic}{\mathbf{d}}^{h(I,x_i)}+{\mathbf{\lambda }}^{(I,x_i)T}\mathbf{G}\left({\mathbf{d}}^{h(I,x_i)}-{\overline{\mathbf{d}}}^{H(I,x_i)}\right),$$

$$\bm{G}=\left[\begin{array}[]{ccccccc}b_{1}&0&0&\,\ldots\ldots&b_{M_{mic}}&0&0\\ 0&b_{1}&0&\,\ldots\ldots&0&b_{M_{mic}}&0\\ 0&0&b_{1}&\,\ldots\ldots&0&0&b_{M_{mic}}\\[5.69054pt] \lx@intercol\hfil\overline{\bm{G}}\hfil\lx@intercol\\ \end{array}\right]\,.$$

$${\mathbf{\lambda }}^{(I,x_i)}=\{{\lambda }_{0,x_1}^{(I,x_i)},{\lambda }_{0,x_2}^{(I,x_i)},{\lambda }_{0,x_3}^{(I,x_i)},{\lambda }_{1,x_1}^{(I,x_i)},{\lambda }_{1,x_2}^{(I,x_i)},{\lambda }_{1,x_3}^{(I,x_i)},\dots ,{\lambda }_{L,x_1}^{(I,x_i)},{\lambda }_{L,x_2}^{(I,x_i)},{\lambda }_{L,x_3}^{(I,x_i)}{\}}^T.$$

Full text

**The Heterogeneous Multiscale Finite Element Method

** for the Homogenization of Linear Elastic Solids

and a Comparison with the FE2 Method

Bernhard Eidel, Andreas Fischer

Heisenberg-Group, Institute of Mechanics, Department Mechanical Engineering

University Siegen, 57068 Siegen, Paul-Bonatz-Str. 9-11, Germany

e-mail: bernhard.eidel@uni-siegen.de, phone: +49 271 740 2224, fax: +49 271 740 2436

Abstract

The Heterogeneous Multiscale Finite Element Method (FE-HMM) is a two-scale FEM based on asymptotic homogenization for solving multiscale partial differential equations. It was introduced in [W. E and B. Engquist, Commun. Math. Sci., 1 (2003), 87–132]. The objective of the present work is an FE-HMM formulation for the homogenization of linear elastic solids in a geometrical linear frame, and doing so, of a vector-valued field problem. A key ingredient of FE-HMM is that macrostiffness is estimated by stiffness sampling on heterogeneous microdomains in terms of a modified quadrature formula, which implies an equivalence of energy densities of the microscale with the macroscale. Beyond this coincidence with the Hill-Mandel condition, which is the cornerstone of the FE2 method, we elaborate a conceptual comparison with the latter method. After developing an algorithmic framework we (i) assess the existing a priori convergence estimates for the micro- and macro-errors in various norms, (ii) verify optimal strategies in uniform micro-macro mesh refinements based on the estimates, (iii) analyze superconvergence properties of FE-HMM, and (iv) compare FE-HMM with FE2 by numerical results.

Keywords: Heterogeneous multiscale method; Finite element method; Homogenization; Macro-to-micro modeling; Estimates; Superconvergence

1 Introduction

Almost all matter is heterogeneous in its structural composition and so are the mechanical properties of solids heterogeneous. As a consequence, deformation and failure of solids and structures at large can be critically influenced by the heterogeneity at small scales, typically referred to as microscales. Homogenization theories provide a framework for the macroscopic modeling of microheterogeneous matter via averaging. Spatial homogenization problems are typically described by partial differential equations (PDEs). For their solution on arbitrary domains and for complex boundary conditions numerical methods must be used. Which among many numerical homogenization methods is best suited for particular problems is still a controversial issue and therefore a most active field of research. Little discussion is about the existing benefits and future promises of numerical homogenization; first, homogenization can reduce the computational costs or can give at all access to problems, which are prohibitive in full microresolution. Second, averaging out small-scale fluctuations of a homogenized solid can provide a better understanding of processes at large. Third, and complementary to the second aspect, the explicit sampling of properties and analysis of processes on representative microdomains allows for predictions on critical phenomena like the initiation of failure of structures at large. Hence, in two-scale numerical schemes one can achieve the one, the description of processes at macroscopic scales without abandoning the other, the added value of small-scale processes, and finally, analyze their interactions.

The Heterogeneous Multiscale Finite Element Method, FE-HMM, is a numerical homogenization method for problems described by PDEs with fast oscillations in the components of the stiffness tensor, the conductivity tensor, or alike. Mathematical homogenization based on asymptotic expansion provides a sound theoretical basis for FE-HMM. It is guiding to a considerable extent the particular formulation of the numerical method in order to preserve theoretical characteristics.

The physical problems treated so far in FE-HMM for spatial 111The most general Heterogeneous Multiscale Method (HMM) introduced in [20] also covers problems living on multiple time scales, see [45], [22]. two-scale problems are diffusion on rough surfaces, [9], stationary and instationary heat conduction, [7], Darcy-flow, [4], and transport problems, [5]. In conclusion, existing FE-HMM implementations and simulations so far have been restricted to scalar-valued field problems.

The progress of FE-HMM and the underlying Heterogeneous Multiscale Method (HMM) since the seminal papers of E, Engquist and coworkers, [20], [21], is described in the overview papers [22] and [10]. The focus in [4] is the fully discrete analysis of multiscale PDEs covering elliptic, parabolic and advection diffusion problems. For elliptic problems, a comparison of the FE-HMM with the multiscale finite element method (MsFEM) as introduced in [35], [36] is presented in

Abdulle presents in [3] the mathematical analysis of FE-HMM for the elliptic case of linear elasticity. The focus is on the fully discrete convergence analysis taking into account the discretization errors at both micro and macro levels. Simulation results are not considered in that reference. Abdulle [3] points out, that no convergence estimates of the overall numerical process taking into account the microscopic and the macroscopic discretization parameters have been given for micro-macro methods, even not for the linear case. Many of the cited micro-macro methods can be classified as FE2 methods, which are equally two-level finite element methods for numerical homogenization. The FE2 method has its origins in solid mechanics, [32], [47], [48], [49], [29], [28], [56], and has found considerable interest in academia and industry; as a versatile method FE2 has been used in non-linear problems of elasticity and inelasticity. For recent, comprehensive overviews of the FE2 method we refer to [31], [62] and [57]. In order to account for size-dependency observed in materials science, Kouznetsova et al. [42], [43] have introduced a second-order homogenization into FE2. Beyond solid mechanics, FE2 has continuously reached out for various physical branches including multifield problems with a coupling of mechanics with thermo-/electro-/magneto-fields, [55], [60], [38], [40] to name but a few.

The theoretical cornerstone of the FE2 method is the macro homogeneity condition or Hill-Mandel condition, [33], [34]. This condition postulates the equality of macroscopic stress power with corresponding stress power on the microscale. Similarly, FE-HMM is a methodology based on a two-scale energy equivalence since its very first, most general conception in the HMM framework [20]. Remarkably, the two methods have not yet been compared in more detail, and to the best of our knowledge, convergence properties for FE2 covering both the microscale and the macroscale error have neither been derived nor have been measured in numerical tests. This is remarkable, since already one decade has past since the work of Abdulle [3] in which the author points out this gap.

In view of the addressed issues and open gaps, the present work has the following main aims:

(i) We derive an FE-HMM formulation for linear elasticity in a geometrical linear frame and address aspects of its implementation. We point to the method’s roots in asymptotic homogenization, which are summarized in the appendix for ready reference. Next, we continue with the FE-HMM characteristics of nonstandard numerical quadrature on the macroscale, which then leads to a detailed algorithmic framework. Therein, the novel aspect is the FE-HMM implementation for a vector-valued field problem putting the method into the context of solid mechanics. The setting of linear elasticity considerably eases a thorough comparison with related methods like FE2, see point (ii) below. In its extension to the vector-valued case of field problems for elliptic PDEs the present work stands on the shoulders of the FE-HMM framework for scalar-valued cases proposed by [9] and [7], [53]. Moreover, the present work considerably takes profit from Abdulle’s analysis in [3] and aims to complement that reference.

(ii) We elaborate a conceptual, algorithmic and numerical comparison of FE-HMM with the FE2 method. The comparison will comprise conceptual aspects like the link of the underlying mathematical theory of asymptotic homogenization in FE-HMM with the physical postulate of energy density equivalence in FE2. It covers the macro-micro kinematical coupling concept and, vice-versa, the micro-to-macro data transfer in terms of stiffness and stress.

(iii) We assess the existing a priori error estimates in various norms ($L^{2}$-, $H^{1}$-, energy-norm) and for the coefficients of the homogenized elasticity tensor. We anticipate that the error analysis including a priori estimates are –for our understanding– the strongest result and most valuable contribution of FE-HMM to the field of numerical homogenization, and the weakest spot of FE2 to go without them. While the existence of error estimates in FE-HMM is remarkable, the content of the estimates is even more exciting, since they seem at a first glance to be at odds with familiar estimates from standard finite element methods. Here, we will analyze the inherent property of FE-HMM to exhibit superconvergence properties for the micro-FEM part in that the error in the $L^{2}$-norm and in the $H^{1}$-norm exhibit the same convergence order. Superconvergence in standard finite element methods, –if present at all– refers to the non-standard property of stress and strain to converge pointwise in the same order as displacements as the primary variables. Superconvergence in standard FEM according to Barlow [13] requires several premises in that it is restricted to particular element shapes and, additionally, is restricted to particular element sites. For the microscale FEM part of FE-HMM in contrast, superconvergence does not require any premises. An additional benefit of error estimates in FE-HMM is that they allow for optimal uniform macro-/micro mesh-refinement strategies – how to refine the micromesh for uniform macro mesh refinement, if the full convergence order shall be achieved but for minimal computational costs?

2 The Heterogeneous Multiscale Finite Element Method

2.1 Model Problem of Linear Elasticity

We consider a body $\mathcal{B}$, a bounded subset of $\mathbb{R}^{d}$, $d=2,3$, with boundary $\partial\mathcal{B}=\partial\mathcal{B}_{D}\cup\partial\mathcal{B}_{N}$ where the Dirichlet boundary $\partial\mathcal{B}_{D}$ and the Neumann boundary $\partial\mathcal{B}_{N}$ are disjoint sets. The closure of the body $\mathcal{B}$ is denoted by $\overline{\mathcal{B}}$. The body, which exhibits an inhomogeneous microstructure, is subject to body forces $\bm{f}$ and surface tractions $\bar{\bm{t}}$ and in static equilibrium.

2.1.1 The microproblem

The displacement $\bm{u}^{\epsilon}=(u_{1}^{\epsilon},\ldots,u_{d}^{\epsilon})$ of the body is given by the solution of the system

[TABLE]

The constitutive law is assumed to be linear elastic where $\mathbb{A}^{\epsilon}_{ijlm}$ is the fourth order elasticity tensor. Superscript $\epsilon$ throughout indicates the dependency of suchlike marked quantities on the heterogeneity of the elastic material. Note that the body forces $\bm{f}$ and the traction vectors $\bar{\bm{t}}$ are assumed not to depend on $\epsilon$. In (1)3, $\bm{n}=(n_{1},\ldots,n_{d})$ is the unit outward normal to $\partial\mathcal{B}$.

For $\mathbb{A}^{\epsilon}_{ijlm}$ the following symmetries hold $\mathbb{A}^{\epsilon}_{ijlm}=\mathbb{A}^{\epsilon}_{jilm}=\mathbb{A}^{\epsilon}_{lmij}$ for any $i,j,l,m=1,\ldots,d$. For the deformation kinematics geometrical linearity is assumed to hold with the linearized strain tensor $\bm{\varepsilon}$

[TABLE]

It can be written in compact format by means of the linear differential operator $\bm{L}$

[TABLE]

For a finite element formulation the strong form (1) is transformed into a variational or weak form. Multiplying the strong form by a test function $\bm{v}\in\mathcal{V}$, using the Green formula yields the following variational formulation:

Find $\bm{u}^{\epsilon}$ such that

[TABLE]

which has to hold for all $\bm{v}\in\mathcal{V}$, where $\mathcal{V}$ is the space of admissible displacements, i.e. virtual displacements that fulfill homogeneous Dirichlet boundary conditions

[TABLE]

The existence and uniqueness of the solution of problem (4) can be shown by use of the first Korn inequality and the Lax-Milgram theorem, [37], [16].

The direct numerical solution of (4) by a standard finite element formulation is prohibitive for small $\epsilon$, since a proper account of the microheterogeneity of characteristic length $\epsilon$ requires an even finer finite element resolution with typical element size $h$, hence $h\ll\epsilon$. This is the main reason for numerical homogenization based on sampling in small regions of confined size instead of an accurate account of microstructure’s heterogeneity everywhere.

The theoretical basis of FE-HMM is mathematical homogenization by asymptotic expansion, [14], [58], [11], [18]. To put things into perspective and for ready reference we provide the main results of asymptotic homogenization for linearized elasticity in the Appendix A.1. These contents are used at various places in the present paper to explain, how FE-HMM is constructed along the lines of asymptotic homogenization.

2.1.2 The macroproblem

The strong form of the macroscopic/homogenized boundary value problem (BVP) reads

[TABLE]

for a derivation see Sec. A.1. The macroscopic displacement is denoted by $u_{i}^{0}$ and $\mathbb{A}^{0}$ is the homogenized elasticity tensor. The term in brackets in (6)1,3 is the macroscopic stress obtained by a volume average over the microdomain, hence it equals $\langle{\bm{\sigma}}^{0}\rangle$, see Sec. A.1.

The values for the Dirichlet as well as Neumann boundary conditions in (6)2,3 are obtained by surface averages according to (A.32) of corresponding boundary conditions in (1)2,3. Similarly, $\langle f_{i}\rangle$ is the volume average of body forces in (1)1.

The solution of the homogenized problem is obtained from the variational form

[TABLE]

which follows from multiplying the strong form (6) by test functions $\bm{v}$ along with the application of Green’s formula. For notational convenience we skip in (7) and in the rest of the paper the averaging symbols $\langle\bullet\rangle$, $\langle\bullet\rangle_{\Gamma}$ for $\bm{f}$, $\bar{\bm{u}}$ and $\bar{\bm{t}}$ but keep in mind that these quantities follow from volume and surface averages, respectively.

There is an alternative route to derive (7); the asymptotic stress expansion (A.5) is plugged into he variational form (4). The resultant weak form expressions are arranged according to their orders $\mathcal{O}\left(\epsilon^{n}\right),n=1,-1,-2$. Each of them is multiplied by $\epsilon^{n}$ with corresponding $n$ along with the limit of $\epsilon\rightarrow 0^{+}$. Doing so the above macroscale variational form can be identified. Hence, it is the same process as carried out for the identification of the macroscopic balance of linear momentum in Sec. A.1.

Next we consider the piecewise linear continuous FEM in macro- and microspace, respectively. The domain $\mathcal{B}$ is a convex polygonal domain in order to avoid regularity issues.

We define a macro finite element space as

[TABLE]

where ${\mathcal{P}}^{p}$ is the space of (in the present work: linear, $p=1$) polynomials on the element $K$, ${\mathcal{T}}_{H}$ the (quasi-uniform) triangulation of $\mathcal{B}\,\subset\,\mathbb{R}^{d}$. Index/superscript $H$ denotes the characteristic element size, with $H\gg\epsilon$ for efficiency. The space $\mathcal{S}^{p}_{\partial\mathcal{B}_{D}}$ is a subspace of $\mathcal{V}$ defined in (5).

For the solution of (1) in the macrodomain we use the two-scale FEM framework of the FE-HMM as originally proposed in [20] and analyzed for elliptic PDEs in [23], and, with the focus on linear elasticity, in [3].

The macrosolution of the FE-HMM is given by the following variational form:

Find $\bm{u}^{H}\in\mathcal{S}_{\mathcal{B}_{D}}(\mathcal{B},\mathcal{T}_{H})$ such that

[TABLE]

which reads as a standard finite element formulation.

2.2 The modified macro bilinear form of FE-HMM

If the homogenized constitutive tensor $\mathbb{A}^{0}(\bm{x})$ is explicitly known, the bilinear form $B_{H}(\bm{u}^{H},\bm{v}^{H})$ can be calculated using standard numerical quadrature according to (10), where $\bm{x}_{K_{l}}$ and $\omega_{K_{l}}$ are the quadrature points and quadrature weights, respectively

[TABLE]

Since $\mathbb{A}^{0}(\bm{x})$ is typically not known for heterogeneous materials, the ansatz of FE-HMM is to approximate the virtual work expression at point $\bm{x}_{K_{l}}$ in the semidiscrete form (10) by another bilinear form using the known microheterogeneous elasticity tensor $\mathbb{A}^{\epsilon}$, see (11). According to this approximation, the solution $\bm{u}_{K_{\delta}}^{h}$ is obtained on microsampling domains $K_{\delta}=\bm{x}_{K_{l}}+\delta\,[-1/2,+1/2]^{d}$, $\delta\geq\epsilon$, which are each centered at the quadrature points $\bm{x}_{K_{l}}$ of $K$, $l=1,\ldots,N_{qp}$. For a visualization see Fig. 1. These microsampling domains with volume $|K_{\delta}|$ provide the additive contribution to the stiffness matrix of the macro finite element.

Remark 1  From a physical point of view, the approximation of (10) by (11) indicates that the FE-HMM is a numerical homogenization scheme that is based on the equality of the macroenergy density with the microenergy density.

From a mathematical point of view, the replacement of the pointwise integrand value by another bilinear form indicates that FE-HMM is a modified numerical quadrature formula. Since it is this setting that enables the fully discrete a priori error analysis, it is a key FE-HMM ingredient. The inserted product of test functions leads to a priori estimates of the convergence order of the microerror, which seem to be at odds with familiar results for FEM; for details see Sec. 3 along with its references to original work.

On Notation  The microdomains attached to the macro quadrature points $l$ in the macro element $K$ are distinguishable by index $l$. In order to avoid too heavy notation we will use $K_{\delta}$ instead of $K_{\delta_{l}}$ in the rest of the paper.

2.3 Variational formulation of the microproblem

It can be shown that the FE-HMM microproblem resembles the discrete version of the cell problem of asymptotic expansion, if it is formulated for each microdomain $K_{\delta}$ in $K$ with $l=1,\ldots,N_{qp}$, $K\in\mathcal{T}_{H}$ like this:

Find $\bm{u}^{h}_{K_{\delta}}$ such that the conditions for macro-micro coupling and for the micro bilinear form (12) are fulfilled:

[TABLE]

where the micro finite element space $\mathcal{S}^{q}(K_{\delta},\mathcal{T}_{h})$ is defined by

[TABLE]

In (13) $\mathcal{T}_{h}$ is a quasi-uniform discretization of the sampling domain $K_{\delta}$ with mesh size $h\ll\varepsilon$ resolving the finescale, $\mathcal{W}(K_{\delta})$ denotes the coupling condition or boundary conditions used for computing the microfunctions, and $\mathcal{P}^{q}$ is the space of polynomials on the element $T$. In the present work we restrict to linear shape functions, hence $q=1$ and consider periodic boundary conditions (PBCs), see (A.27).

The linearization of $\bm{u}^{H}$ in (12)1 is carried out at the quadrature point $\bm{x}_{K_{l}}$

[TABLE]

Remark 2  Note that (12) resembles for PBC the discrete version of the cell problem of asymptotic homogenization, (A.26), if $\bm{u}^{h}_{K_{\delta}}-\bm{u}^{H}_{lin,K_{\delta}}=\bm{\chi}^{h}_{K_{\delta}}$ and if $\bm{u}^{H}$ follows the linearization of (14). Then, replacing $\bm{u}^{h}_{K_{\delta}}$, (12)2 reads as

[TABLE]

The linearized $\bm{u}^{H}_{K_{\delta}}$ is the FE-HMM counterpart of the quantity $\bm{I}^{lm}$ in asymptotic homogenization; they both induce a homogeneous deformation in the microdomain/unit cell and (15) coincides with (A.26). Furthermore, we hint already here at the coincidence of FE-HMM with the FE2 method with respect to the above superposition of a homogeneous deformation with periodic fluctuations, which is frequently referred to as strain-driven, first order computational homogenization.

2.4 Bases for spaces

For the solution of (12) a basis $\{N_{I}^{H}\}_{I=1}^{M_{mac}}$ for the macro finite element space $\mathcal{S}^{1}_{0}(\mathcal{B},\mathcal{T}_{H})$ is employed in order to represent the macrosolution $\bm{u}^{H}$ of (9). Similarly, a basis $\{N_{i}^{h}\}_{i=1}^{M_{mic}}$ of the micro finite element space $\mathcal{S}^{1}_{0}(K_{\delta},\mathcal{T}_{h})$, (13), is introduced in order to represent the solution $\bm{u}^{h}$ of a microproblem. In (16), $M_{mac}$ denotes the number of nodes of the macrodomain, and $M_{mic}$ in (17) denotes the number of nodes of each microdomain. Hence, the macro- and the microsolution follow the representation

[TABLE]

where $\bm{d}_{I}^{H}$ is the displacement vector of macronode $I$, and $\bm{d}_{i}^{h}$ is the displacement vector for micronode $i$.

The node-based shape functions $\{N_{I/i}^{H/h}\}$ exhibit the standard properties (i) partition of unity: $\sum_{j\in\mathcal{N}}N_{j}(\bm{X}_{i})\equiv 1\quad\forall\,\bm{X}_{i}\in\mathcal{B}/K_{l}$ and (ii) compact support $N_{j}(\bm{X}_{j^{\prime}})=\delta_{jj^{\prime}}\quad\forall\,j,j^{\prime}\in\mathcal{N}$, where $\mathcal{N}$ is the set of finite element nodes.

2.5 Macrostiffness calculation

According to (11), FE-HMM can be seen as a bottom-up multiscale method in that the required information for estimating stiffness on the macrolevel is obtained from microsampling domains $K_{\delta},l=1,\ldots,N_{qp}$.

The macro bilinear form $B^{e}_{H}(\bm{u}^{H},\bm{v}^{H})$ is the virtual work of internal forces in a macro finite element $e$. The corresponding bilinear form in terms of the shape functions $B^{e}_{H}(\bm{N}_{I}^{H},\bm{N}_{J}^{H})$ extracts the stiffness matrix contribution $\bm{k}^{e,mac}_{IJ}$ for macronodes $I,J$, a $d\times d$ matrix. Consequently, we replace in (10) the displacements $\bm{u}^{H},\bm{v}^{H}$ by macro shape function matrices $\bm{N}_{I}^{H},\bm{N}_{J}^{H}$ and obtain222Note that matrix $\bm{N}^{H}_{I}$ is a $d\times d$ diagonal matrix with elements $N^{H}_{I}$; hence for $d=3$ it holds $\bm{N}^{H}_{I}=\mbox{diag}(N^{H}_{I},N^{H}_{I},N^{H}_{I})$. Similarly, matrix $\bm{u}^{h(I)}$ in (18) is a $d\times d$ diagonal matrix diag$(u_{x},u_{y},u_{z})$.

[TABLE]

In (18) $\bm{u}^{h(I)}_{K_{\delta}}$ is the counterpart of $\bm{u}_{K_{\delta}}^{h}$ in (12). It is the dimensionless solution of the microproblem on $K_{\delta}$, which is driven by the shape function $N_{I}^{H}$ at macronode $I$. In the following, we add $x_{i},i=1,\ldots,d$ to account for the vector-valued field problem of dimension $d$. Consequently, $\bm{u}^{h(I,x_{i})}_{K_{\delta}}$ is the microsolution driven by a macroelement unit-displacement state $\bm{u}^{H(I,x_{i})}_{lin,K_{\delta}}$ at node $I$ in $x_{i}$-direction.

For stiffness calculation, problem (12) is to be reformulated as follows: Find $\bm{u}^{h(I,x_{i})}_{K_{\delta}}$ on each $K_{\delta}$ such that the conditions for the kinematical coupling and for the micro bilinear form are fulfilled

[TABLE]

In order to realize the coupling of the macrodisplacement field ${\bm{u}}^{H(I,x_{i})}_{lin,K_{\delta}}$ with the yet unknown microdisplacement field $\bm{u}^{h(I,x_{i})}_{K_{\delta}}$, the two fields are expanded into the same basis $\{N_{i}^{h}\}_{i=1}^{M_{mic}}$ of $\mathcal{S}^{1}(K_{\delta},\mathcal{T}_{h})$,

[TABLE]

The solution of the microproblems for the minimizers $\bm{d}^{h(I,x_{i})}$ is presented in Sec. 2.6. We continue with the calculation of the macroelement stiffness matrix in (18) and obtain

[TABLE]

where $\bm{d}^{h(I)}=\left(\,\bm{d}^{h(I,x_{1})}|\bm{d}^{h(I,x_{2})}|\bm{d}^{h(I,x_{3})}\,\right)$ for $d=3$. The assembly of $\bm{k}^{e,mac}_{IJ}$ yielding $\bm{k}^{e,mac}$ implies the arrangement of $\bm{d}^{h(I)}$ in different columns for $I=1,\ldots,N_{node}$, which yields the transformation matrix $\bm{T}_{K_{\delta}}$. Doing so we obtain

[TABLE]

Equations (23) and (24) give insight into important characteristics of $\bm{T}_{K_{l}}$. First, the matrix dimensions

[TABLE]

reveal that $\bm{T}_{K_{\delta}}$ ”compresses”333Following the terminology of HMM, $\bm{T}_{K_{\delta}}$ can be called a stiffness compression operator, see [20], Sec.2.1. the total micro stiffness matrix $\bm{K}^{mic}_{K_{\delta}}$ to the dimensions of the macroelement stiffness matrix $\bm{k}^{e,mac}_{K}$. For that reason $\bm{T}_{K_{\delta}}$ is a micro-macro stiffness transfer operator along with a model reduction/coarse-graining. Second, the transformation matrix $\bm{T}_{K_{\delta}}$ is built up by the column vectors $\bm{d}^{h(I,x_{i})}$, $x_{i}\,|\,i=1,\ldots,d$, $I=1,\ldots,N_{node}$.

The macroelement shape functions each represent a unit displacement state for each macro node $I$ in each direction of space $x_{i}\,|\,i=1,\ldots,d$. They drive the microproblem in terms of the corresponding nodal values ${\bm{d}}^{H(I,x_{i})}_{m},m=1,\ldots,M_{mic}$ via (20) in each microdomain to evaluate the macroelement stiffness $\bm{k}^{e,mac}_{IJ}$. Each unit displacement state in $x_{i}$-direction induces in ${\bm{d}}^{H(I,x_{i})}$ nonzero components only in $x_{i}$, for example ${\bm{d}}^{H(I,x_{i})}|_{i=2}=\left[0,{d}^{H(I,x_{2})}_{1,x_{2}},0,\ldots,0,{d}^{H(I,x_{2})}_{M_{mic},x_{2}},0\right]^{T}$.

Figure 2 (left) visualizes a unit displacement state, the corresponding uniform deformation of microdomains and (right) the microsolution obtained for the additional constraint of PBC444The problem deals with a stiff inclusion in a softer matrix, analyzed in detail in Sec.5.1..

2.6 Solution of the microproblems

While the element stiffness matrix of a macro finite element is computed by the modified quadrature formula expanding the collocation points to sampling domains, the element stiffness matrix $\bm{k}^{e,mic}_{mn}$ of a micro finite element is computed by standard Gauss-Legendre quadrature in the $\xi$-, $\eta$-, $\zeta$-parameter space on the unit cube (of side length 2)

[TABLE]

where $\bm{B}^{e}_{m}=\bm{L}\,\bm{N}_{m}$ and where $\bm{J}(\xi_{k},\eta_{k},\zeta_{k})$ is the Jacobian for an isoparametric finite element formulation.

The assembly of the element microstiffness matrices on the sampling domain $K_{\delta}$ yields the total microstiffness matrix

[TABLE]

required for the solution of a microproblem.

The problem (19) can be rewritten as a minimization problem

[TABLE]

over all functions $\bm{u}_{K_{\delta}}^{h(I,x_{i})}\in\mathcal{S}^{1}(K_{\delta},\mathcal{T}_{h})$ satisfying $\bm{u}_{K_{\delta}}^{h(I,x_{i})}-{\bm{u}}_{K_{\delta}}^{H(I,x_{i})}\in\mathcal{S}^{1}_{per}(K_{\delta},\mathcal{T}_{h})$.

For the solution of the microproblem (28), we choose the method of Lagrange multipliers and thus arrive at a saddlepoint problem, see for example [15]. It is convenient for our purposes to express the strain energy by the micro domain stiffness matrix and the microdisplacement vector.

Doing so, the problem reads: Minimize for each node on a macro finite element $I=1,\ldots,N_{node}$ and for each direction of space $x_{i}\,|\,i=1,\ldots,d$ the functional

[TABLE]

where $\bm{G}$ contains the kinematical coupling conditions and reads for $d=3$

[TABLE]

In $\bm{G}$, the first $d$ rows contain the normalization condition for the fluctuations, the macro-micro coupling condition, to which the first $d$ Lagrange multipliers are associated; for details see (33)–(2.6). The submatrix $\overline{\bm{G}}$ contains the $L$ non-redundant periodic coupling conditions for adjacent nodes at opposite edges/faces in 2d/3d on the microdomain boundary with the associated $L\cdot d$ Lagrange multipliers. The vector of Lagrange multipliers $\bm{\lambda}^{(I,x_{i})}\in\mathbb{R}^{(1+L)\cdot d}$ reads for $d=3$

[TABLE]

The first variations of $\mathcal{L}$ with respect to its variables give the stationarity condition

[TABLE]

This set of linear equations is to be solved for $\bm{d}^{h(I,x_{i})}$ and $\bm{\lambda}^{(I,x_{i})}$. Since the coefficient matrix in (32) is constant for all macro unit displacement states $(I,x_{i})$, the solution vectors as well as the right-hand sides for $I=1,\ldots,N_{node}$ and $x_{i},i=1,\ldots,d$ are augmented to full matrices, hence, $\bm{d}^{h(I,x_{i})}\rightarrow\bm{T}$, $\bm{\lambda}^{(I,x_{i})}\rightarrow\bm{\Lambda}$, ${\bm{d}}^{H(I,x_{i})}\rightarrow\bm{d}^{H}$.

The coefficients $b_{i}$, $i=1,\ldots,M_{mic}$ in (30) follow from the kinematical macro-micro coupling condition, (19)1, which can be written as

[TABLE]

Since (33) is the discrete counterpart of the normalization condition (A.21) for the periodic fluctuations, the particular choice of the constant is inconsequential for the microsolution. Here we choose $\bm{c}=\bm{0}$.

With (21), it holds

[TABLE]

Using the method of Lagrange multipliers, the PBC are not strongly imposed on the FE space but rather weakly enforced. In [9] it is found that despite of this weak coupling of scales, the error estimates of [23] still hold, where in the latter reference, the PBC are strongly imposed on the FE space.

2.7 The homogenized elasticity tensor and its efficient computation

The direct stiffness sampling on microdomains circumvents the necessity to calculate the homogenized elasticity tensor $\mathbb{A}^{0}$ and therefore the necessity of an explicit macroscopic constitutive model. Notwithstanding, calculating $\mathbb{A}^{0}$ in terms of its numerical approximation $\mathbb{A}^{0,h}$ offers a second route to numerical homogenzation. In its derivation we follow [3] and propose an efficient way of its computation.

A reference simplex $\hat{K}$ with its local vectorial basis $\hat{N}^{H}_{0,j}=(1-\xi_{1}-...-\xi_{d})\otimes e_{j},\hat{N}^{H}_{i,j}=\xi_{i}\otimes e_{j},i,j=1,...,d$ is considered. With $\hat{u}_{i,j}^{h}$ being the solution to (2.3) such that $\hat{u}_{i,j}^{h}-\hat{N}_{i,j}^{H}\in\mathcal{S}_{\text{per}}(K_{\delta},\mathcal{T}_{h})$ we can compute the numerical homogenized elasticity tensor as follows

[TABLE]

With

[TABLE]

(37) yields the coefficients of the homogenized elasticity tensor

[TABLE]

The error introduced by approximating $\mathbb{A}_{ijlm}^{0}$ by $\mathbb{A}_{ijlm}^{0,h}$ follows $|\mathbb{A}_{ijlm}^{0,h}-\mathbb{A}_{ijlm}^{0}|\leq C(h/\epsilon)^{2q}$, see [4], Sec. 3.2.2.

The homogenized elasticity tensor is calculated by inserting selected macro shape functions of a reference element in the variational formulation of the microproblem, which leads to a macroscopic unit strain state as shown in (38). By doing so the homogenized elasticity tensor can be calculated coefficientwise by evaluating the corresponding strains on the microdomain.

The homogenized elasticity tensor $\mathbb{A}^{0,h}$ can be calculated alongside the FE-HMM calculation most efficiently from the already available results for the calculation of the macroelement stiffness matrix. For that purpose the microdisplacements are used, which follow from the unit displacement states of a macroelement.

While the notation above is limited to simplex reference elements, the application of a unit strain state can be easily expanded to other elements. The unit displacement states have to be combined and multiplied by a corresponding factor, such that they yield the corresponding unit strain states on the macrodomain.

Figure 3 exemplarily visualizes this procedure for a quadrilateral. Two unit displacement states of the nodes on the right macroelement side are superimposed. The resultant constant strain state is multiplied with the element length and results in a unitary strain state.

Remark 3  The described computation of the homogenized elasticity tensor is only valid for linear elasticity. For nonlinear elasticity, the microstrains have to be calculated from the microdisplacements, which result from the macro unit displacement state.

2.8 The macrosolution

The global macroquantities, the stiffness matrix $\bm{K}^{mac}$, the force vector $\bm{F}$ and the displacement vector $\bm{D}$ follow from a standard assembly of the corresponding macroelement matrices/vectors

[TABLE]

Prescribed displacements $\bm{D}=\bm{D}_{d}$ on Dirichlet boundaries $\partial\mathcal{B}_{D}$ are eliminated by static condensation

[TABLE]

The first set of equations in (41) gives the desired macro solution $\bm{K}_{11}\,\bm{D}=\bm{F}_{I}-\bm{K}_{12}\bm{D}_{d}$.

2.9 Postprocessing

After the solution of the macroproblem, the microproblem can be finally solved. Now, the microproblems in a macro finite element are driven by the true, linearized macrodisplacements, $\bm{u}^{H}_{lin,K}$, which results in the microsolution vector $\bm{u}^{h}$. The postprocessing is carried out in the following steps

[TABLE]

The true resultant microdisplacements are obtained by the superposition of the macrosolution with the microsolution at the centre of the microdomain $\bm{x}_{K_{l}}$ according to (43). From the macrosolution $\bm{u}^{H}$ we obtain the macrostrains $\bm{\varepsilon}(\bm{u}^{H})$, (44), and from the resultant microsolution $\bm{u}^{Hh}$ the microstrains $\bm{\varepsilon}(\bm{u}^{Hh})$, (45).

Next, micro stresses $\bm{\sigma}(\bm{u}^{Hh})$, (46), and the macro stresses $\bm{\sigma}(\bm{u}^{H})$, (47), are calculated.

The second equality in (46) is due to (43). The equality of macrostress according to (47) with (48) follows from the fact that macrostress is defined as the volume average of microstress in asymptotic homogenization, and that the homogenized elasticity tensor is identified in this averaging process, see Sec. A.1.

2.10 Algorithm

For convenience, we summarize the FE-HMM algorithm for the geometrically linear, linear elastic case in Tab. 1. Step 1 can be understood as a stiffness preprocessor using already the microsolver. Step 2 is referred to as the macrosolver, step 3 uses again the microsolver for the postprocessing.

3 Convergence and a priori error estimates

For the fully discrete case [3] provides convergence results. In the periodic case, the macrosolution $\bm{u}^{H}$ approximates the homogenized solution $\bm{u}^{0}$, see Theorem 4.1 in [3]. Since $\bm{u}^{0}$ does not capture the small scale oscillations of $\bm{u}^{\epsilon}$, its numerical approximation $\bm{u}^{H}$ equally can not. But similar to the first order correction according to (A.38) and the corresponding enriched approximation (A.39), the reconstructed microsolution $\bm{u}^{Hh}$ according to (43) approximates the true fine scale solution $\bm{u}^{\epsilon}$ in the $H^{1}$-norm, see Theorem 4.2. in [3].

3.1 A priori estimates

In view of reliable solution data for engineering decisions, a key property of numerical methods is the existence of error estimates with guaranteed bounds.

FE-HMM exhibits a priori as well as a posteriori estimates for various types of partial differential equations. For elliptic PDEs we refer to [23] and [54], for the (elliptic) case of linear elasticity in a geometrical linear setting to [3]. A posteriori error indicators along with adaptive mesh refinement have been presented in [8] and [6].

We restrict to a priori estimates for the elliptic problem of linear elastic solids in a geometrical linear frame provided in [3]. In Sec. 5 we assess these estimates for uniform discretizations.

The total error of the FE-HMM-method can be decomposed into three parts

[TABLE]

where $e_{mac}$, $e_{mod}$, $e_{mic}$ are the macro error, the modeling error, and the micro error.

Here, $\bm{u}^{0}$ is the solution of the homogenized equation, (6), $\bm{u}^{H}$ is the FE-HMM solution, $\bm{u}^{0,H}$ is the standard FEM solution of problem (7), therefore obtained through exact $\mathbb{A}^{0}$; and $\widetilde{\bm{u}}^{H}$ is the FE-HMM solution obtained through exact microfunctions (in $W(K_{\delta})$).

For sufficiently regular problems, the following a priori estimates hold in the $L^{2}$-norm, the $H^{1}$-norm and the energy-norm:

[TABLE]

Remark 4

(i) Abdulle [4] underlines that for the fully discrete analysis of FE-HMM the setting of the modified quadrature formula according to (11) is essential. The reason is that the error analysis crucially relies on conditions (ellipticity condition and approximation condition, first studied by Ciarlet and Raviart, see [17]) which ensure that a FEM with numerical quadrature converges with the same rate to the exact solution than the same FEM with exact integration.

(ii) Equations (50) and (51) provide relations of the macro error and the micro error. They enable strategies to achieve the optimal convergence order for minimal computational costs in micro-macro uniform discretizations. Details will be discussed in Sec. 5 dealing with numerical examples.

3.2 Superconvergence

The convergence order of $2q$ of the FE-HMM micro error according to (50) is already in the $L^{2}$-norm at odds with familiar results of standard FEM. The fact however, that the micro error in the $H^{1}$-norm scales in the same order as in the $L^{2}$-norm is even more unusual, since it indicates the same order of a derivative as the primary quantity. This phenomenon is called superconvergence. We contrast this type of standard FE-HMM superconvergence with the (non-standard) superconvergence in single-scale FEM. For the latter Barlow [13] has shown that superconvergence exists only at particular sites of finite elements of rectangular shape and, moreover, that the number and locus of these points depend on the polynomial order of the shape functions. Zienkiewicz and Zhu presented in [65], [66] a procedure for the transfer of the superconvergence property from superconvergent, inner element points to element nodes referred to as ”superconvergent patch recovery” (SPR). Based on these superconvergent nodal values the same authors constructed a refinement indicator for adaptive remeshing.

Here, we analyze superconvergence for the numerical homogenization by FE-HMM.

3.2.1 Macro FEM and micro FEM

The macro error of FE-HMM in the $L^{2}$-norm (here: of displacements) scales in the order $\mathcal{O}(H^{p+1})$. The corresponding macro error in the $H^{1}$-norm and in the energy norm (hence of strain and stress) scale in the order $\mathcal{O}(H^{p})$, one order below the displacements.

Following the analysis of [13], for $p=1$, superconvergence of strain and stress –hence order $p+1=2$– can be expected in the element center of macro finite elements, if they exhibit rectangular shapes, which implies that the error in the $H^{1}$-norm and in the energy norm exhibit the same convergence order as the error in the $L^{2}$-norm.

Since the micro errors of FE-HMM in the $H^{1}$-norm and in the energy norm exhibit the same convergence order as the error in the $L^{2}$-norm, $\mathcal{O}((h/\varepsilon)^{2q})$, superconvergence is standard for the micropart of FE-HMM, which therefore is not restricted to rectangular element shapes and not to Barlow’s nominally superconvergent points.

3.2.2 The superconvergent patch recovery for the macro FEM

In order to verify superconvergence of the macro FEM in the FE-HMM context, the rationale of the SPR according to [65] is applied; strain and stress are calculated at superconvergent element sites, for $p=1$ in the center of an element. Next, these values are transferred by a least-square procedure to the finite element node in the direct neighborhood. Elements having such a node in common are referred to as the patch in the superconvergent recovery procedure. For a visualization see Fig. 4.

Stresses on the patch are prescribed component-wise by

[TABLE]

with

[TABLE]

Vector $\bm{P}$ contains polynomial terms of the ansatz functions, here for $d=2$ and the case of linear quadrilaterals. For the determination of the unknown vector $\bm{a}$ the function

[TABLE]

has to be minimized. Therein, $(x_{i},\,y_{i})$ are the coordinates of the superconvergent points, $n$ is the number of superconvergent points of the total patch and $\sigma_{h}(x_{i},\,y_{i})$ are the stresses in these superconvergent points. Minimization of $F(\bm{a})$ implies that $\bm{a}$ fulfills the condition

[TABLE]

which can be solved for $\bm{a}$

[TABLE]

with

[TABLE]

Stresses in the central node of the patch can be recovered by inserting its nodal coordinates $(x_{N},\,y_{N})$ into the $\bm{P}$-Vector in (53).

Summarizing, the numerical assessment of the superconvergence property of stresses and strains can be calculated as described in Tab. 2

4 Conceptual Comparison of FE-HMM with the FE2 Method

The aim of the present section is a conceptual comparison of the FE-HMM with the FE2 method, which reveals similarities, but also differences in some noticable aspects.

4.1 The FE2-Method

To put things into perspective we briefly reiterate the building blocks of the FE2 method.

The basic idea of computational homogenization in the FE2 method is to solve at each integration point of the macroproblem the associated microproblem. For the case of the first-order strain-driven computational homogenization the macroproblem provides input to the problem at the microscale in terms of macroscopic deformation and, vice-versa, the microproblem provides stiffness information in terms of the macroscopic tangent and stress. The latter is obtained from averaging microstresses, which implies that a constitutive law does not exist on the macroscale.

The theoretical cornerstone of the FE2 method is the Hill-Mandel condition, which postulates that the stress power on the macrolevel must be equal to the average microscopic stress power in a so-called Representative Volume Element, RVE

[TABLE]

where macrostress $\overline{\bm{\sigma}}$ is calculated as the volumetric mean of microstresses $\bm{\sigma}$ according to

[TABLE]

Boundary conditions for the RVE, which are consistent to the Hill-Mandel postulate, are the constraint condition $\dot{\bm{\varepsilon}}:=\dot{\overline{\bm{\varepsilon}}}$ on the whole RVE (Voigt condition), the constraint condition of constant stress $\bm{\sigma}=\overline{\bm{\sigma}}$ on the whole RVE (Reuss condition). Moreover, linear Dirichlet and linear Neumann conditions each fulfill the condition as well. Finally, PBC are consistent with the postulate. They are visualized in Fig. 6. Macrostrain $\overline{\bm{\varepsilon}}$ yields a homogeneous deformation on the RVE, which is superimposed by periodic fluctutations $\widetilde{\bm{w}}$.

For a proper definition of PBC, the boundary of the RVE is divided into pairwise periodic parts $\partial\text{RVE}^{+}$ and $\partial\text{RVE}^{-}$ with a corresponding sign convention for position vectors $\bm{x}^{+}$ and $\bm{x}^{-}$ and normal vectors $\bm{n}^{+}=-\bm{n}^{-}$. For pairwise periodicity of the displacements, i.e. $\widetilde{\bm{w}}^{+}=\widetilde{\bm{w}}^{-}=\widetilde{\bm{w}}$ along with

[TABLE]

yields

[TABLE]

which implies the additional condition $\widetilde{\bm{t}}^{+}=-\widetilde{\bm{t}}^{-}$. Thus, the PBC satisfying the macrohomogeneity condition can be given as

[TABLE]

The homogenized elasticity tensor $\overline{\mathbb{A}}$ is calculated for each macro quadrature point. If static condensation is applied, $\overline{\mathbb{A}}$ can be calculated according to

[TABLE]

with

[TABLE]

The corresponding system of linear equations which has to be solved to obtain $\overline{\mathbb{A}}$ is

[TABLE]

The number of right hand sides follows from the dimension of $\overline{\bm{L}}$, hence depends on ${\bm{B}}$ acccording to (65)2,3. The number is 3 for $d=2$, and 6 for $d=3$.

4.2 Comparison

4.2.1 Homogeneity condition. Energy equivalence

While the FE2-method is built on the Hill-Mandel condition (59), FE-HMM is built on asymptotic expansion. It is a straightforward excercise to show that for PBC, asymptotic homogenization theory fulfills the Hill-Mandel condition. With the strain influence function already introduced in (A.23) a direct calculation of microstrain from macrostrain is obtained

[TABLE]

With (67) we obtain

[TABLE]

For the evaluation of the second integral in (4.2.1) we have used integration by parts

[TABLE]

with $\bm{n}$ being the unit normal to the boundary. Here we exploited periodicity, namely that $\chi_{i}^{kl}$ is equal on the opposite site of the unit cell, whereas surface tractions $t_{i}=\sigma_{ij}n_{j}$ are equal and of opposite sign on opposite sites of the unit cell. The second integral in (69) vanishes, since it holds ${\partial\sigma_{ij}}/{\partial y_{j}}=0$ for equilibrium.

Equivalence of energy densities in FE-HMM is enforced in the frame of a modified quadrature rule for the calculation of the macrostiffness matrix by stiffness sampling on microdomains, see (11). Already in its very first, most general proposal, the Heterogeneous Multiscale Method (HMM) for variational problems was built on energy equivalence in terms of a fine-scale to coarse-scale stiffness transfer, see Sec. 2.2 in [20]. For that reason we can call the FE-HMM energy equivalence ”E-macrohomogeneity condition” for the initials of energy and the authors.

4.2.2 Kinematics of coupling conditions

Both FE-HMM and FE2 attach the microdomains to the macrodomain at the quadrature points of macro finite elements as visualized in Fig. 1. The explicit kinematical coupling follow for FE2 from the energy equivalence condition (59), and for FE-HMM from the cell problem of asymptotic homogenization. Periodic boundary conditions meet these requirements such that the macro-micro kinematical coupling effectively coincides for FE2 and FE-HMM; in FE2 strain at the macro quadrature point $\overline{\bm{\varepsilon}}$ induces a homogeneous deformation in the attached RVE with the corresponding displacement field $\bm{u}=\overline{\bm{\varepsilon}}\cdot\bm{X}$. In FE-HMM the homogeneous deformation in the microdomain is induced by the linearized macrodisplacement field $\bm{u}^{H}_{lin}$, which is superimposed by periodic fluctuations following from constrained energy minimization in the RVE and microdomain, respectively. In the above nomenclature, the periodic fluctuations are referred to as $\widetilde{\bm{w}}$ in FE2, in FE-HMM as the difference $\bm{u}^{h}-\bm{u}^{H}_{lin}$, and in asymptotic homogenization as $\bm{\chi}^{lm}$.

The enforcement of coupling conditions have been realized in the FE2 method either by static condensation/elimination, [47], by the method of Lagrange multipliers, [50], or by the Penalty method, [63]. In FE-HMM, the coupling was realized by Lagrange multipliers in [9], [7], [27]. For a discussion of the role of different types of boundary conditions on the RVE we refer to [64], [44].

4.2.3 Micro-to-macro coupling/bottom-up data transfer

In strain-driven two-level finite element frameworks, the microscale must pass over stress and stiffness information to the macroscale. For the particular case of linear elasticity, the deformation-dependent tangential material stiffness boils down to the constant elasticity tensor.

Stiffness transfer in FE2 is carried out by the effective fourth-order constitutive tensor $\overline{\mathbb{A}}$, although its computation may differ. It can be carried out by numerical perturbation as proposed in [28], by Lagrange multiplier techniques, [50], or by static condensation, [47]. The effective constitutive tensor is passed over to the corresponding macro quadrature point, where it enters the macro element stiffness matrix.

For FE-HMM in practice, the bottom-up stiffness transfer can be carried out by two different means, either by (i) a direct stiffness matrix transfer mediated by micro-to-macro transformation matrices or by (ii) the homogenized constitutive (here: elasticity) tensor, where the latter coincides with the FE2 method. We anticipate that the two variants of stiffness transfer show quantative agreement, see Sec. 5. Notwithstanding, it is only the first approach of the stiffness matrix transfer which enables due to its underlying modified quadrature formula the derivation of unified a priori estimates and therefore is a key, indispensible FE-HMM ingredient, see Sec. 4.2.5.

4.2.4 Data-driven versus model-driven method

Closely related to the bottom-up data transfer is the discrimination between data- and model-driven methods. Within the FE2 method, a macroscopic constitutive model is not defined since not required. Instead, the macroscale is fed by homogenized microdata in terms of stress and stiffness, see Sec.4.2.3. In this respect FE2 is a data-driven method, in view of the absent macro constitutive equations it can even be called an ”equation-free” method.

Similarly, the sole constitutive input for FE-HMM is $\mathbb{A}^{\epsilon}$, which lives on the microscale. It is passed over to the macroscale via the transformed micro stiffness matrix, and macrostresses are calculated by averaging microstresses. Of course, the macroscale solver is a finite element method but the physical model is not declared.

4.2.5 A priori error estimates

FE-HMM is endowed with a priori estimates for the macro error and the micro error in the $L^{1}$-, $H^{2}$- and the energy-norm as well as with estimates for the coefficients of the homogenized elasticity tensor, see Sec. 3.1.

While the a priori estimates of the macro error in FE-HMM directly follow from the approximation condition for standard FEM, [17], the a priori estimate for the micro error and its propagation to the macroscale crucially relies on the FE-HMM setting as a modified quadrature rule and on the method’s roots in asymptotic homogenization, see [2], [3]. FE2 lacks these roots and, as a consequence, no unified a priori estimate covering both the macro and the micro error exists. Moreover, there are, to the best of our knowledge, even no numerical convergence tests available for FE2.

Unified a priori estimates allow in FE-HMM for micro-macro uniform refinement strategies, which are optimal in that the full convergence order can be achieved for a minimum of computational costs. This is a considerable benefit of FE-HMM compared with FE2. In numerical tests related to scalar-valued field problems it was shown that the a priori estimates are quite sharp, [4], [7]. For the present case of linear elasticity, the a priori estimates and derived optimal micro-macro uniform refinement strategies will be assessed in Sec. 5.

4.2.6 Numerical costs

The method of Lagrange multipliers for stiffness calculation in FE-HMM requires to solve the microproblem for $N_{node}\times d$ right hand sides. Since the microstiffness matrix is independent of the unit displacement states, it is a linear set of equations with $N_{node}\times d$ different right hand sides thus reducing the computational complexity. In FE2 using the static condensation method, the number of right hand sides in the microproblem is 3 for $d=2$, and 6 for $d=3$, hence cheaper than the FE-HMM approach using Lagrange multipliers. Of course, the solution of the microproblems in FE2 by the method of Lagrange multipliers is equally more expensive than static condensation, see [50]. Moreover, if static condensation is used in FE-HMM, the numerical costs are the same as in FE2. If FE-HMM takes the route using the homogenized elasticity tensor $\mathbb{A}^{0,h}$, the numerical effort is exactly the same as in FE2.

In either case, linear problems, here for both geometrical linearity and material linearity, are solved not only by FE-HMM but also by FE2 in one single step, no iterations are required, see [3], Sec. 1; for that reason the algorithm in Tab. 1 similarly applies to the FE2 method.

Beyond parallelization several approaches can reduce the computational complexity of two-level FEM. Methods based on Fast-Fourier Transforms (FFT) as introduced in [52] use a direct, point-wise discretization of the Lippmann-Schwinger equation, with improvements in [47] and [59]. In an effort to reduce the computation time by a modified modeling [61] introduce into FE2 the rationale of a statistically similar representative volume element (SSRVE) in order to replace the true microstructure in its full geometrical complexity by a surrogate, which resembles the original one by geometrical features as analyzed by different geometrical similarity measures.

4.2.7 Limitations

The classical FE2 method has proven to be equally suitable for linear and non-linear problems. It has been applied to non-linear elasticity, to various inelastic constitutive models, and to multifield problems with a coupling of mechanics with thermo-/electro-/magneto-fields. The versatility of FE2 is its big plus; the extension of FE-HMM to finite deformation non-linear elasticity is a non-trivial task, most notably to provide the mathematical analysis in terms of estimates for that case.

In the frame of first-order computational homogenization methods FE2 and FE-HMM exhibit inherent limitations. Inhomogeneous deformation modes such as bending cannot be properly represented on the microscale. As a consequence, the microscale is not appropriate to capture geometrical size effects in material behavior. This issue was overcome by second order computational homogenization as introduced into FE2 by [42], [43]. For a thorough discussion of the limitations and future challenges we refer to [30]. Alternative to the second-order homogenization to capture the size effect, one can introduce surface/interface energies into homogenization which are more intuitive and can be physically motivated as described in [19] and [39].

4.2.8 Relation to other methods

There is a remarkable connection of FE2 and FE-HMM to the fully nonlocal, cluster-based Quasicontinuum Method (CQC), which was introduced by [41]. In CQC, nanosampling domains in the shape of spherical clusters are used to approximate the sum over all lattice sites of a crystal by a weighted sub-sum. For a modification to an energy-based approach of this concept endowing the method with a variational structure see [24], [25]. The nanosampling domains of CQC are strongly coupled to the displacement field of the finite elements. Additionally, the clusters are coupled by the nonlocal interactions of the atoms in the sampling domain with other atoms within the cutoff-radius. The concept of scale-separation is a necessary theoretical requirement in FE2 and FE-HMM. Quite in contrast, CQC as an atomistic-continuum coupling method aims at a seamless transition between the scales. For CQC, the sole constitutive input originates from the fine scale (potential of atomic interactions), and consequently, no explicit coarse-scale constitutive model exists. It is remarkable that in this particular characteristics, CQC is close to FE2 and FE-HMM, but is in contrast to other atomistic-continuum coupling methods, which equally require a coarse-scale/macro model, [26].

In conclusion, innovative quadrature formula (FE-HMM) and summation rules (CQC), see [12], have been advancing multiscale modeling to a considerable extent.

5 Numerical Examples

This section assesses the numerical performance of the FE-HMM method. The results and corresponding errors are measured in the following norms

[TABLE]

Four benchmark problems subject to plane strain conditions will be analyzed,

1. (i)

a matrix-inclusion problem, 2. (ii)

a microstructure with an analytical solution for homogenization, 3. (iii)

a non-uniformly periodic microstructure, 4. (iv)

a uniformly periodic tessellation inspired by M.C. Escher.

Key aspects are

1. (A)

the convergence order of the simulation results, where the errors are calculated by means of an accurate reference solution (”overkill”-solution) on very fine meshes. Doing so, the a priori estimates will be assessed. 2. (B)

Based on the unified error estimates covering the macro error and the micro error, tests will be carried out on how full convergence orders can be achieved for minimal costs in uniform micro-macro mesh refinements. 3. (C)

Superconvergence employing the superconvergent patch recovery will be tested. 4. (D)

Simulation results of FE-HMM are compared with those of FE2, thus complementing the conceptual comparison of Sec. 4.2. 5. (E)

Finally, the reconstruction of the microsolution from the macrosolution shall demonstrate the capability of postprocessing small scale features, which can be used to investigate local phenomena like stress concentrations, which eventually may initiate inelastic deformations and failure mechanisms.

To improve the readability but still properly document the present simulation results, we provide the numerical convergence data in tabular form as a supplement and restrict here to the more telling convergence diagrams.

The elasticity tensor will be given in standard Voigt notation. Moreover, we replace the $x_{i}$, $i=1,\ldots,d$ coordinate system by a $x$-,$y$-,$z$-coordinate system.

5.1 Matrix-inclusion problem for beam-bending

In the first numerical example we consider the microstructure of a stiff inclusion in a soft matrix, hence a problem with uniformly periodic coefficients, $\delta=\epsilon$, such that periodic boundary conditions are applied to the microdomains. The Young’s moduli of the inclusion $E_{i}=100000$ $\text{N/mm}^{2}$ and the matrix phase $E_{m}=40000$ $\text{N/mm}^{2}$ exhibit the contrast of $E_{i}/E_{m}=2.5$, for the Poisson’s ratio it holds $\nu=0.2$. The volume ratio of the inclusion phase is $V_{i}/V_{tot}=1/16$.

The beam as displayed in Fig. 7 exhibits length $l_{1}=5000\leavevmode\nobreak\ \text{mm}$, height $l_{2}=1000\leavevmode\nobreak\ \text{mm}$, and thickness $t=100\leavevmode\nobreak\ \text{mm}$. The external line load is $f=1\leavevmode\nobreak\ \text{N/mm}^{2}$. The side length of the square microcell is $\epsilon=5\leavevmode\nobreak\ \text{mm}$, hence $\epsilon/l_{1}=5/5000\ll 1$ and $\epsilon/l_{2}=5/1000\ll 1$.

5.1.1 Macro- and microconvergence

The simulation results are obtained along with a reference solution with 3200 $\times$ 640 macroelements and $32\times 32$ microelements. The results turn out to be invariant with respect to a variation of $\epsilon$. Figure 8 (left) shows the total error as a function of the macrodiscretization while keeping the microdiscretization fixed; the mean convergence orders in the $L^{2}$-norm of almost 2, and in the $H^{1}$- and in the energy-norm of almost 1 confirm the a priori estimates of Sec. 3.1.

Similarly, the convergence of the micro error in right diagram of Fig. 8 is in reasonable agreement with the theoretical predictions. The deviation (1.72 versus 2) is due to the magnitude of the stiffness jump at the inclusion-matrix interface, which lowers the regularity of the microproblem. This assertion is underpinned by the local convergence distribution of the error in the energy-norm as shown in the right of Fig. 9; at the matrix-inclusion interphase, most notably at the corners, the convergence order is considerably reduced. The reference solution for microerror calculation is obtained using $200\times 40$ macroelements and $2048\times 2048$ microelements.

5.1.2 Homogenized elasticity tensor

Next, we assess the convergence estimate for the components of the approximate homogenized elasticity tensor $\mathbb{A}^{0,h}$ where the reference solution $\mathbb{A}^{0}$ is obtained on a microdomain discretized by $2048\times 2048$ elements. The convergence as a function of $h/\epsilon$ is shown in the left diagram of Fig. 9. The mean slope varies between 1.70 for coefficient $\mathbb{A}_{33}$ and 1.84 for coefficient $\mathbb{A}_{12}$ and well agrees with the convergence of the microerrors in the $L^{2}$-norm and the $H^{1}$-norm. The deviation from the theoretical order 2 is caused by the aforementioned reduced regularity of the microproblem for its material stiffness contrast.

5.1.3 Comparison of FE-HMM with FE2

In practice, the use of $\mathbb{A}^{0,h}$ in the FE-HMM simulation is an alternative route to the standard FE-HMM solution method. The difference between the solutions is quantitatively analyzed in Tab. 3; the normed percental deviation is throughout less than $5\times 10^{-10}$ and independent of the macrodiscretization.

Since the FE2 method draws in the homogenization process on the fourth order constitutive tensor, here $\mathbb{A}^{0,h}$, it coincides with the FE-HMM byproduct of the homogenized elasticity tensor. The simulation results for FE-HMM along with the micro-macro stiffness matrix transfer show quantitative agreement with the FE2, which underpin the identity of the methods by numerical means as well.

5.1.4 Postprocessing on selected microdomains

Figure 10 displays the von-Mises stress distributions on two different microdomains, in the lower left and in the upper right of the beam. The stiff inclusion attracts stress, while the stress in the soft matrix is considerably smaller. Right at the interface of the inclusion with the matrix, a stress jump can be observed; the fact that the stiffer inclusion carries over large stresses shielding the softer matrix in terms of the matrix’ strong local reduction of stress is referred to as stress-shielding.

Refined information is obtained by the contour plots of relevant stress components and strain components for the two points considered in the lower left of the beam, Fig. 11, and for the upper right of the beam, Fig. 12.

Remark 5  The $H^{1}$-norm depends on the chosen physical unit of the primary quantity; the smaller the (length) unit, the larger the weight of the primary quantity compared to the derivative. In the limiting case the $H^{1}$-norm equals the $L^{2}$-norm, and the macro error in the $H^{1}$-norm will show the same order as in the $L^{2}$-norm – in contrast to the estimates. In order to avoid this falsifying effect, we use in the first example of Sec. 5.1 the length unit of kilometer instead of millimeter for evaluating the $H^{1}$-norm.

5.2 Assessment of superconvergence in the FE-HMM

Since the present example is uniformly periodic, the homogenized elasticity tensor is constant everywhere including the superconvergent element center.

5.2.1 Superconvergence for the macro FEM

The error in the energy-norm is computed at mesh nodes of the coarsest discretization, here using a $25\times 5$ mesh. Then, a hierarchical, uniform mesh refinement is carried out making 4 elements out of one of the coarsest mesh. The reference solution is obtained for $800\times 160$ elements. The microdomain constantly exhibits a discretization of $32\times 32$ elements, with $\epsilon=5\,\text{mm}$.

Figure 13 displays the local convergence order of the error in the energy-norm as a contour plot in the macrodomain. Obviously, superconvergence shows up in the macro error of FE-HMM in terms of order two (instead of one according to the error estimates of standard FE-HMM) almost everywhere. A reduced order is observed at the sites of load application and reaction forces; the variations beyond and below order 2 are most pronounced at the clamped end.

Figure 14 displays the convergence of the error of the stress components $\sigma_{xx}$ and $\sigma_{xy}$ as a function of the macroelement size. Component $\sigma_{yy}$ was discarded for its virtually vanishing magnitude. The convergence order of 2.09 demonstrates superconvergence of stresses at that particular node.

The integral error in the energy-norm is shown in Tab. 4. Superconvergence is not obtained. Instead we see the standard convergence order of one. The reason for the lost superconvergence is that for the integral error all nodes of the macrodomain are considered. The nodes at the re-entrant corners of the clamped edge exhibit only orders in the range of 0.5, see Fig. 13, which spoils the overall convergence order.

5.3 Problem with an analytical solution of the homogenized tensor

The BVP of the present example is a square plate of sidelength $l=1\,m$ and thickness $t=0.1\,m$, which is clamped at one edge and loaded by a constant line load of $q_{0}=1$ $[\text{N/mm}^{2}]$ at the opposite edge, see Fig. 15.

The microheterogeneous elasticity tensor is uniformly periodic

[TABLE]

with components in $[N/m^{2}]$. It was used in [7] for heat conduction in a slightly different format.

In contrast to the matrix-inclusion problem, the microheterogeneous tensor exhibits an analytical solution, such that FE-HMM can here be assessed by an exact solution without drawing to an overkill solution.

According to [18] and [8] the exact homogenized tensor $\mathbb{A}^{0}$ and the exact volumetric mean $\langle\mathbb{A}^{\epsilon}\rangle$ is given by

[TABLE]

The comparison of ${A}^{0}$ with $\langle\mathbb{A}^{\epsilon}\rangle$ shows that the volumetric mean cannot account for the induced anisotropy due to the laminate character of the microstructure; component $\langle\mathbb{A}^{\epsilon}_{11}\rangle$ considerably deviates from the exact solution.

5.3.1 Macro- and microconvergence

The simulation results for the macroconvergence are obtained based on a reference solution using 2560 $\times$ 2560 macroelements and $32\times 32$ microelements. The results are invariant with respect to $\epsilon$.

Figure 16 (left) displays the total error as a function of the macrodiscretization for a constant microdiscretization. The observed convergence order in the $L^{2}$-norm of 1.62 instead of 2, and in the $H^{1}$-norm as well as the energy-norm of almost 0.84 instead of 1 are in reasonable agreement with the a priori estimates. The deviation from the theoretical order can be traced back to notch effects at the clamped end.

The convergence of the micro error in the diagram of Fig. 16 (centre) for various $h/\epsilon$ keeping $H$ fixed confirms the theoretical predictions in quantitative agreement. The reference microsolution is obtained using $50\times 50$ macroelements and $2048\times 2048$ microelements.

5.3.2 Homogenized elasticity tensor

The obtained results by FE-HMM can capture the exact values of the above elasticity tensors in good agreement. A microdiscretization using $160\times 160$ elements yields

[TABLE]

The convergence of the components of the approximate homogenized elasticity tensor $\mathbb{A}^{0,h}$ is displayed in the right diagram of Fig. 16. The slope indicates order 2 in agreement with the theory, $|\mathbb{A}^{0,h}-\mathbb{A}^{0}|\leq C(h/\epsilon)^{2}$. The reference solution $\mathbb{A}^{0}$ is obtained on a microdomain with a discretization of 2560 $\times$ 2560 elements. Note that component $\mathbb{A}_{22}^{0,h}$ does not change in the leading 9 digits for meshes finer than 40 $\times$ 40 elements, the reason, why the diagram in Fig. 16 is restricted to component $\mathbb{A}_{11}^{0,h}$.

5.3.3 Comparison of FE-HMM with FE2

Table 5 shows, that FE-HMM and FE2 exhibit virtually the same results.

5.3.4 Optimal refinement strategies

The a priori error estimates ($\delta/\epsilon\in\mathbb{N}$ along with periodic coupling) enable an optimal uniform refinement strategy for uniform macro- ($H$) and micro- ($h$) meshes, see [4]. For $p=q=1$ i.e. bilinear shape functions for the macro- as well as on the micro-FE, we obtain

[TABLE]

We denote by $L$ the number of micro elements in each direction of space of the microdomain. Similarly, $M$ is the number of macro elements per direction. Consequently, the micro element size is $h=\epsilon/L$ and it holds $h_{L}:=h/\epsilon=1/L$ for the micromesh, and $H_{M}=1/M$ for the macromesh.

Denoting by $N_{mac}$ the number of macro DOF and by $N_{mic}$ the number of micro DOF, the above rates of convergence show that the following choices are the best uniform refinement strategies that is full order is achieved for minimal computational costs.

[TABLE]

We compute the solutions of the problem for uniform macro mesh refinement $H_{M}=1/8,1/12,1/16,\ldots,1/1024$. In the diagrams of Fig. 17 (right) and Fig. 18 (right) the micromesh is kept fixed for each of the distinct solid lines and is successively refined from one line to the other. Optimal refinements clearly follow the ratio $h_{L}=H_{M}$ for the error in the $L^{2}$-norm, Fig. 18, and $h_{L}=\sqrt{H_{L}}$ for the error in the $H^{1}$-norm, Fig. 17. The results demonstrate the sharpness of the a priori bounds. Similar results for piecewise bilinear FE are reported in [4].

5.4 Non-uniformly periodic tensor

The present example serves the purpose to analyze FE-HMM for a non-uniformly periodic microstructure, which is given by a Young’s modulus of

[TABLE]

as displayed for $\epsilon=0.005$ m in Fig. 19 on the domain $[0,1]$ m $\times$ $[0,1]$ m, see [7]. The coordinate system is centered in the domain, axes align with the boundaries.

The macro BVP is almost the same as in the previous example of Sec. 5.3. Only the line load at the free edge is adopted to $q_{0}=0.01$ N/m, since the Young’s modulus is significantly smaller than in the previous example.

5.4.1 Macro- and microconvergence

The convergence of the macro error for a fixed microdiscretization is displayed in the left diagram of Fig. 20. The corresponding reference solution is obtained by $2560\times 2560$ macroelements and $40\times 40$ microelements. The convergence order is below the theoretical order in the $L^{2}$-norm (1.66 instead of 2), and in the $H^{1}$- as well as in the energy-norm (0.88 instead of 1), which is due to the reduced regularity of the macro BVP, compare the results with those in Sec. 5.3.

The convergence of the micro error for a fixed macrodiscretization is displayed in the centre diagram of Fig. 20. Full order of convergence is obtained for the micro error in the $L^{2}$-norm and in the $H^{1}$-norm, respectively. The reference solution for microerror calculation is obtained by $40\times 40$ macroelements and $2048\times 2048$ microelements.

5.4.2 Homogenized elasticity tensor

Next, we assess the convergence estimate for the components of the approximate homogenized elasticity tensor $\mathbb{A}^{0,h}$. The convergence as a function of the microdiscretization is shown in the right diagram of Fig. 20. The slope indicates order 2 for all components of the elasticity tensor in the plane strain case in perfect agreement with the theory, [4].

5.4.3 Comparison of the FE-HMM variants, and with the FE2 method

The deviation of the FE-HMM solution based on micro-macro stiffness transfer in comparison with the second version using the homogenized elasticity tensor is throughout less than $4\cdot 10^{-11}$ and thus in the range of machine precision, see Tab. 6. Since the latter FE-HMM version coincides with FE2, the coincidence of FE-HMM and FE2 is underpinned again by numerical means in the present non-uniformly periodic setting.

5.5 Escher’s Bird and Fish tessellation

While Lord Kelvin was concerned with the problem, how to pack equal-sized objects together to fill space with minimal surface area, the graphic artist M.C. Escher (1898–1972) was devoted to the question, how to fill the plane by tessellations555Tessellations are arrangements of closed shapes that completely cover the plane without overlapping and without leaving gaps. with most intriguing beauty.

While from 1887 until 1993 the tetrakaidecahedron was continuously but erroneously believed666The so-called Weaire-Phelan structure was shown in 1993 to solve the Kelvin problem with even smaller surface area. to be the right answer of Kelvin to his minimization problem, the answer to Escher’s question is simple in that it is clearly non-unique at any time, as beauty is in the eye of the beholder777*”At moments of great enthusiasm it seems to me that no one in the world has ever made something this beautiful and important.”* – M.C. Escher.

It is a well-known mathematical result that of all the regular polygons, only the triangle, square, and hexagon can be used for a tessellation. Escher combined them applying reflections, glide reflections, translations, and rotations to augment his playground for creating a larger variety of tessellations.

From Escher’s oeuvre in the field which he called ”regular divisions of the plane” we devote our analysis to the ”Bird/Fish (1941 B)” tessellation, an artistic artifact that is uniformly periodic and consists of two phases, Fig. 21. Here we choose $\delta=\varepsilon=28.5$ mm. The macroscopic BVP is the same as in Sec. 5.1, and the Poisson’s ratio and Young’s moduli exhibit the same values, thus showing a contrast of the bright face to the dark face of $E_{b}/E_{d}=2.5$.

5.5.1 Macro and micro convergence

Figure 22 shows in its left diagram the convergence of the normed error for macroscopic mesh refinement for a reference solution of 4 096 000 DOFs (3200$\times$640 macroelements) and with 1462 DOFs on each microdomain. Similarly, Fig. 22 shows in its right diagram the convergence of the normed errors for microscopic mesh refinement for a reference solution of 16 000 DOFs on the macrodomain (200$\times$40 macroelements) and 312 506 DOFs on each microdomain.

The convergence rates for $p=q=1$ reasonably agree with the theoretical values. The deviations are in the reduced regularity of the problem sets on the two scales.

At the end, we abstain from homogenizing (the elasticities of) bird and fish and discover the microscale simulation results displayed in Fig. 23.

6 Summary

The novel contributions and main results shall be summarized.

Based on the seminal contributions of E, Engquist and Abdulle we have presented a finite element an FE-HMM formulation for linear elasticity in a geometrical linear frame. For that case of a vector-valued field problem we have generalized the FE-HMM micro-to-macro stiffness transfer by deriving corresponding transformation matrices which follow from macro element unit displacement states and result in energy minimizers. By virtue of the present setting in solid mechanics, a clear mechanical interpretation of FE-HMM as an energy-based homogenization scheme is facilitated. Moreover, this setting opens the door to a conceptual and numerical comparison with other two-scale finite element homogenization concepts. 2. 2.

A conceptual comparison of the FE-HMM with the FE2 method was presented; they both have in common the equivalence of micro- with macroenergy densities; it was shown that for periodic boundary conditions the macro-homogeneity condition of Hill-Mandel as the theoretical cornerstone of FE2 can be derived from asymptotic homogenization as the theoretical foundation of FE-HMM, which is the deepest common root of the two methods. FE2 in its micro-to-macro stiffness transfer draws on the homogenized fourth order constitutive tensor while keeping the macroscopic bilinear form unaltered. FE-HMM however crucially relies on a modified quadrature formula in the macroscopic bilinear form, which results in a micro-macro stiffness matrix transfer, a setting which enables the a priori error analysis. For that reason the stiffness matrix transfer is a key ingredient of FE-HMM. The homogenized elasticity tensor –coinciding with the one of FE2– must be seen as a byproduct.

The existence of a priori estimates is the particular strength and key novel contribution of FE-HMM to the field of numerical homogenization in general, and in particular to two-scale finite element methods like FE2. In the mathematical cornerstones of FE-HMM Therein are the asymptotic homogenization framework and the setting of numerical quadrature. Here, FE-HMM and its protagonists, W. E, B. Engquist, A. Abdulle and coworkers have brought mathematical rigor into the fully discrete analysis of two-scale finite element methods for homogenization. For FE2 methods, in spite of their considerable progress in a multitude of applications, a thorough mathematical analysis is missing. 3. 3.

In representative benchmark problems FE-HMM has been assessed for periodic boundary conditions imposed on the microdomains. Here we checked in various norms the a priori error estimates as given by [23] for the general case of elliptic PDEs and for linear elasticity by [3]. For sufficient regularity of the BVP, the estimates have been verified for both the macro error and the micro error and for the coefficients of the homogenized elasticity tensor.

Remarkably, superconvergence holds for the micro error in FE-HMM in that stress (and equally the error in the energy-norm and in the $H^{1}$-norm) converges in the same order as displacements (and equally the error in the $L^{2}$-norm). In contrast to the notion of superconvergence in standard, single-scale finite element methods, which holds only at particular points for particular element shapes, the verified estimates for the micro error are generally valid. We have demonstrated by the application of the Superconvergent Patch Recovery (SPR) technique, that for the macro error in FE-HMM standard superconvergence holds at the so-called Barlow points.

Moreover, the existing predictions for FE-HMM concerning the optimal uniform micro-macro refinement strategies on how to obtain full order for minimal computational costs are corroborated by the numerical results in the present work. 4. 4.

Furthermore, the simulations have shown quantitative agreement of the two FE-HMM variants with FE2. Of course, for the coincidence of the two methods in practice, FE-HMM findings equally apply for the FE2 method with respect to (i) the micro and macro convergence orders including superconvergence and with respect to (ii) the derived optimal uniform micro-macro mesh refinement strategies.

In conclusion, the present work shall contribute to bridge the still existing gap between the mathematics of multiscale modeling with FE-HMM as an excellent example, and, on the other hand side, the engineering community, in which the popular FE2 method has been growing into different branches of solid mechanics by virtue of its versatility. Here, FE-HMM as a non-standard discretization method, which is endowed with a fully-fledged mathematical structure at least for the linear case, advances numerical homogenization more than by an $\epsilon$ towards reliability as a key ingredient of scientific computations.

Acknowledgements. The first author acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) within the Heisenberg program (grant no. EI 453$/$2-1).

Appendix A Mathematical basis of FE-HMM

A.1 Asymptotic homogenization

This part presents the derivation and main results of asymptotic homogenization for linearized elasticity. It is a main reference point to underpin the mathematical basis of FE-HMM. For mathematical homogenization by asymptotic expansion we refer to [14], [58], [11], [18].

In the two-scale mathematical homogenization theory, various fields are assumed to depend on two coordinates, $\bm{x}$ as the coarse- or macroscale position vector, and $\bm{y}$ as the fine- or microscale position vector. These two vectors are related by $\bm{y}=\bm{x}/\epsilon$ with $0<\epsilon\ll 1$.

The unit cell domain $Y$ is typically chosen a rectangular domain, in which the microstructure is assumed to be locally periodic ($Y$-periodic). Even if the microstructure is not periodic, the response functions $\bm{r}^{\varepsilon}$ are in either case assumed to be locally periodic.

In addition to the symmetries introduced in Sec. 2.1.1 we assume for $\mathbb{A}^{\epsilon}_{ijlm}$

[TABLE]

The derivative of a response function $\bm{r}^{\epsilon}(\bm{x})$ with respect to the macroscale variable $x_{i}$ can be calculated by the chain rule according to

[TABLE]

For the case of linear elasticity we consider the yet unknown response functions $\bm{u}^{\epsilon}$, $\varepsilon_{ij}^{\epsilon}$, $\sigma_{ij}^{\epsilon}$ and their expansions into series

[TABLE]

The strain expansion (A.4) is obtained by inserting the expansion for $u_{i}^{\epsilon}(\bm{x})$ into (2) with an account of (A.2) and with the following nomenclature for the strain components of various ”orders”

[TABLE]

where the abbreviation $u^{k}_{(ij)y}=1/2(\partial u_{i}^{k}/\partial y_{j}+\partial u_{j}^{k}/\partial y_{i})$ (and analogously for $u^{k}_{(ij)x}$) is employed. For each scale an elasticity relation is assumed to hold

[TABLE]

Inserting (A.5) into the balance of linear momentum (1)1 yields

[TABLE]

from which, according to the order, different balance equations follow (by multiplying (A.8) by $\epsilon^{n}$ and considering the limit $\epsilon\rightarrow 0^{+}$, each for $n=2,1,0$)

[TABLE]

For the Dirichlet and Neumann boundary conditions (1)2,3 we obtain by (A.3) and (A.5)

[TABLE]

The boundary conditions for various orders of $\epsilon$ are identified from (A.12) and (A.13) according to

[TABLE]

[TABLE]

It can be shown, that $\partial u^{0}_{i}/\partial y_{j}=0$, which implies $\bm{u}^{0}=\bm{u}(\bm{x})$ as already anticipated in (A.3), and moreover, that $\sigma_{ij}^{-1}=0$. Thus, the first term in the expansion of $\bm{u}^{\epsilon}$, i.e. $\bm{u}^{0}$, is the coarse-scale/macroscale displacement field, and according to (A.14)1 and (A.15)1 the Dirichlet as well as Neumann boundary conditions (A.12), (A.13) refer to the macroproblem.

For $\sigma_{ij}^{-1}=0$ it follows from the $\mathcal{O}\left(\epsilon^{-1}\right)$ equilibrium, (A.10), that

[TABLE]

which yields along with the kinematical relations (A.6)2 and the elastic constitutive law (A.7) for $s=0$

[TABLE]

A separation of variables $\bm{x}$ and $\bm{y}$ in (A.17) is achieved by the ansatz

[TABLE]

where $\chi_{k}^{lm}(\bm{y})$ is the so-called first-order displacement influence function, which is locally periodic ($Y$-periodic) and symmetric $\chi_{k}^{lm}(\bm{y})=\chi_{k}^{ml}(\bm{y})$. Inserting (A.18) into (A.17) along with the identity $u^{0}_{(kh)x}=\mathbb{I}_{khlm}u^{0}_{(lm)x}$ where $\mathbb{I}_{khlm}=1/2\,(\delta_{kl}\delta_{hm}+\delta_{km}\delta_{hl})$ results in

[TABLE]

which must hold for arbitrary macro strain $u^{0}_{(kh)x}$ leading to the cell problem, which reads:

Find for any $l,m\in\{1,\ldots,d\}$ the functions ${\chi}_{k}^{lm}(\bm{y})$ as the solutions of the system

[TABLE]

Since the field ${\chi}_{k}^{lm}$ enters the cell problem by $\nabla_{y}\bm{\chi}^{lm}$, it is determined up to a constant vector. For definiteness, however, uniqueness can be achieved by imposing the normalization condition

[TABLE]

It can be shown, that if (A.21) holds, the microstrain can be calculated from macrostrain by means of an elastic strain influence function $E_{kl}^{mn}(\bm{y})$

[TABLE]

with

[TABLE]

thereby realizing a separation of variables similar to (A.18).

With an eye towards the representation of the microproblem in the discrete case of FE-HMM we prefer to rewrite the cell problem (A.20) as

[TABLE]

In (A.24) $\bm{I}^{lm}$ is a vector-valued function given by

[TABLE]

The weak form of the cell problem in the format of (A.24) reads

[TABLE]

where $\mathcal{W}_{per}$ is given by

[TABLE]

Averaging the microscale elasticity law, i.e. (A.7) for $s=0$, $\sigma^{0}_{ij}=\mathbb{A}^{\epsilon}_{ijkh}\left(u^{0}_{(kh)x}+u^{1}_{(kh)y}\right)$ with an account of (A.18) leads to

[TABLE]

where $\langle{\sigma}^{0}_{ij}\rangle=1/|Y|\int_{Y}\sigma^{0}_{ij}(\bm{x},\bm{y})\,d\bm{y}$ and where the homogenized tensor $\mathbb{A}^{0}(\bm{x})=(\mathbb{A}^{0}_{ijkl})$ is identified as

[TABLE]

which equally verifies the above symmetries of elasticity.

Integration the second term in (A.11) over the unit cell domain yields

[TABLE]

With (A.30), the boundary conditions (A.14)1 and (A.15)1, integration of the $\mathcal{O}\left(1\right)$ balance equation (A.11) over the unit cell results in the the macroscopic BVP

[TABLE]

In (A.31)1 $\langle f_{i}\rangle$ is defined as the volume average of $f_{i}$, $\langle{\bm{f}(\bm{x})}\rangle=1/|Y|\int_{Y}\bm{f}(\bm{x},\bm{y})\,d\bm{y}$. In (A.31)2,3 $\langle\bullet\rangle_{\Gamma}$ defines a surface average according to

[TABLE]

where $\Gamma$ is that particular surface area ($\Gamma\in\partial\mathcal{B}$) to which the averaging is applied. Hence, $\langle\bullet\rangle_{\Gamma}$ is the surface counterpart to the volume average $\langle{\sigma}^{0}_{ij}\rangle$.

Remarks  Three remarks are in order. First, that $\bm{I}^{lm}$ is to be understood as a displacement vector which induces by its definition a unitary strain $\bm{\varepsilon}(\bm{I}^{lm})$ with components $\varepsilon_{ij}^{lm}=1/2(\partial I_{i}^{lm}/\partial y_{j}+\partial I_{j}^{lm}/\partial y_{i})$. The corresponding strain vectors in Voigt notation read

[TABLE]

where $\bm{\varepsilon}^{lm}=(\varepsilon^{lm}_{11},\varepsilon^{lm}_{22},\varepsilon^{lm}_{33},2\varepsilon^{lm}_{12},2\varepsilon^{lm}_{13},2\varepsilon^{lm}_{23})^{T}$, [46]. In analogy to the homogenization of a second order PDE for a scalar-valued field problem, [3], $(\bm{\varepsilon}^{lm})_{l,m=1}^{d}$ is the canonical basis of the symmetric, six-dimensional strain space.

Second, introducing with $\bm{U}^{lm}:=\bm{\chi}^{lm}+\bm{I}^{lm}$ the superposition of the $Y$-periodic fluctuations with the (prescribed) displacement field that induces a unitary, homogeneous strain on the cell, it follows for (A.29)

[TABLE]

Hence, component $\mathbb{A}^{0}_{ijlm}$ of the homogenized elasticity tensor equals the averaged micro stress for the solution $\bm{U}^{lm}$ of the cell problem.

Third, a note on notation. As indicated by its three indices $\chi_{k}^{lm}$ is a third order tensor; in (A.19) for example, the partial derivative of $\chi_{k}^{lm}$ with respect to $y_{h}$ yields a fourth order tensor, which is in dimensional consistency with the fourth order unity tensor in (A.19). In (A.24) however, ${\bm{\chi}}^{lm}$ as well as $\bm{I}^{lm}$ are to be understood as vector-valued functions, where the superscripts $lm$ are frozen, and only one index is active. Then the calculation of strain from a displacement field as e.g. in (A.24) is a reasonable operation. Hence, the change in notation shall reflect this characteristics.

A.2 General convergence

From homogenization theory the following convergence results and estimates are known.

The microdisplacements $\bm{u}^{\epsilon}$, solution of (1), converge for $\epsilon\rightarrow 0$ –usually in a weak sense– to $\bm{u}^{0}$, the solution of the homogenized problem (7), see [18], Sec.10. For the corresponding strong error estimate it holds

[TABLE] 2. 2.

Due to the oscillations of the fine scale solution, strong error estimates in the $H^{1}$-norm can usually not be obtained since the gradients of the oscillations are in general not $\mathcal{O}(\epsilon)$ quantities. The homogenized solution requires a correction through information of the fine scale. This can be done defining the corrector given by

[TABLE]

where $\omega^{l}$ is a $d\times d$ matrix given by $\omega^{l}=(\omega^{l}_{km})_{1\leq k,m\leq d}=\left(\chi^{lm}_{k}\right)_{1\leq k,m\leq d}$, and where $\bm{\chi}^{lm}(\bm{x},\bm{x}/\epsilon)$ is the solution of the cell problem. With this correction an error estimate in the $H^{1}$-norm can be obtained

[TABLE]

where a boundary layer term yields a $\sqrt{\epsilon}$ instead of a $\varepsilon$ convergence rate, for a thorough analysis see [18], Sec. 4.1.

Appendix B Numerical Data for Simulation Examples