On the effective properties of foams in the framework of the couple stress theory

In the framework of the couple stress theory, we discuss the effective elastic properties of a metal open-cell foam. In this theory, we have the couple stress tensor, but the microrotations are fully described by displacements. To this end, we performed calculations for a representative volume element which give the matrices of elastic moduli relating stress and stress tensors with strain and microcurvature tensors.

The aim of this paper is to discuss the homogenization technique of foams in the framework of the couple stress theory as discussed in [40]. A real geometry of an alumina foam specimen was obtained using the tomography technique. As an open-cell foam has quite complex microstructure advanced technologies as tomography is required to get real geometry as well as extensive calculations, see, e.g., [40][41][42][43][44][45]. As a result, we obtain the complete set of anisotropic material parameters which can be used for further analysis of this foam materials.
The paper is organized as follows. In Sect. 2, we recall in brief the governing equations of the couple stress theory. In order to demonstrate the used approach, we consider a benchmark solution for an isotropic material in Sect. 3. Here, we consider deformations of a solid cube for different boundary conditions. The derivation of the effective properties is based on Hill-Mandel lemma, that is on the comparison of the mechanical energies of the non-homogeneous material and the effective media [46,47]. To this end, we propose new set of consistent kinematic boundary conditions. Finally, in Sect. 4 we apply this technique for the real tomography-based structures. Here, we obtain 21 elastic moduli for a given metal foam. As a result, we get the complete set of material parameters for the considered foam in the framework of the couple stress theory.

Constitutive equations of the couple stress theory
The kinematics of a Cosserat-type medium is based on two kinematical descriptors that are the fields of displacements and microrotations where x is the position vector of material particles. Unlike Cosserat continuum where φ is independent on u, in the couple stress theory φ depends on u as follows [29,40] where ∇ is the three-dimensional nabla operator and × stands for the cross product.
In the couple stress theory, we introduce two symmetric strain measures which are strain and microcurvature tensors, respectively. For anisotropic media with central symmetry, there is a strain energy density given by the formula where C and D are the fourth-order tensors of elastic moduli and : denotes the inner product in the space of second-order tensors. For (4), we get the stress-strain relations in the form where σ and μ are the stress and couple stress tensors, respectively. The form of C and D was discussed in [30,48,49] for various material symmetries of micropolar solids. Let us only note that this analysis should be modified as here χ is a traceless tensor.

Benchmark solution
For the determination of effective material properties using the homogenization, the volume of the representative volume element plays a very important role. Another important issue is the type of the loading which can be considering in the form of constant traction, kinematic assumptions (prescribed displacements) or a periodic boundary conditions (PBCs). With the increased size of representative volume element, these three types should lead to the same solution. However, a large size of RVE requires solving large boundary-value problems which is both expensive and time-consuming. In most of the problems, the periodic boundary conditions provide the best convergence in the determination of the effective material data (unless the analyzed structure is periodical).
In the analysis of the high-porosity foams, constant traction may generate locally very large deformations of the foam skeleton which results in plasticity and/or buckling phenomena. Because of this strain localization, this type of BCs is not considered here. Thus, in this research prescribed displacements (kinematic load) are assumed. The general approach to determination of the effective properties of porous media in the framework of the couple stress theory has been recently proposed by Goda et al. [40]. In [40], twenty-one numerical tests were performed in order to determine twenty-one couple-stress material elastic parameters. In what follows, we use the approach of [40] with some modifications related to boundary conditions. In order to discuss the proposed approach in more details, we consider a benchmark test for an isotropic homogeneous material and performed all numerical tests for it.
Following Hill-Mandel lemma, here we compare stored energies of a non-homogeneous material with stored energies of heterogeneous material with effective elastic moduli. For example, assuming that displacements are linear functions of coordinates such as where ε is a constant tensor and dot means the scalar product, or in Cartesian coordinates we get that χ = 0 and strain energy U of RVE takes the form where V is the volume of RVE. So choosing one or two components of ε i j equal to 1, we can easily calculate the corresponding elastic moduli. To this end, we have to applied proper boundary conditions on the faces of RVE.
With C in hands, we can determine D in a similar way. Again, assuming χ i j constant we can consider the following field of rotations φ φ = χ · x, (9) or, in the coordinate form, we get the average energy in the form where ... denotes the averaged (mean) value. Again, assuming simple form of χ i j we can obtain from calculated energy the value of corresponding D i jkl . Let us note that here we cannot use both ε i j and χ i j constant as in the case of the couple stress theory constant χ corresponds to rotations given as linear functions of coordinates. As the rotations relate to the displacements through (2), in this case we have a quadratic dependence for u and linear one for ε.
As an example, we consider steel with the following properties: Young modulus E = 2e5 MPa and Poisson's ratio ν = 0.33. In Voigt's notation, C takes the form of 6 by 6 matrix are the bulk and shear moduli. So we have K + 4/3G = 2.9633e5 MPa, K − 2/3G = 1.4595e5 MPa, and G = 7.5188e4 MPa. Comparison of given material parameters and obtained ones in numerical computations can prove the reliability of the proposed strain energy approach. In the numerical simulations, the 1 × 1 × 1 mm cube is considered, see Fig. 1. So the volume of the RVE is V = 1 mm 3 . The Cartesian coordinate system is located in the RVE centroid.
Here, we used an idea to compare stored energies of a non-homogeneous material with stored energies of heterogeneous material with effective elastic moduli.
1. For computation of C 11 , we assume an uniaxial tension with a uniform strain ε x x = 1. To this end, we assume the following boundary conditions u x = x on n x faces, u y = 0 on n y faces, u z = 0 on n z faces.
The corresponding shape is given in Fig. 2. The strain energy from the FEM analysis and the corresponding modulus are U FEM = 148165 mJ, The latter value coincides with the value of K + 4/3G. 2. For computation of C 22 , we assume an uniaxial tension with a uniform strain ε yy = 1. Here, we assume the following boundary conditions u x = 0 on n x faces, u y = y on n y faces, u z = 0 on n z faces.
The calculated strain energy and the corresponding modulus are 3. Computation of C 33 -uniaxial tension with a uniform strain ε zz = 1. Assumed boundary conditions (BCs) are u x = 0 on n x faces, u y = 0 on n y faces, u z = z on n z faces.
Here, we get again U FEM = 148165 mJ, In the last two cases, moduli coincide again with the value of K + 4/3G. 4. Computation of C 12 -biaxial tension with a uniform strains ε x x = ε yy = 1 (Fig. 3). BCs are u x = x on n x faces, u y = y on n y faces, u z = 0 on n z faces.
This results in which coincides with the value of K − 2/3G. 5. Computation of C 23 -biaxial tension with a uniform strains ε yy = ε zz = 1. BCs are u x = 0 on n x faces, u y = y on n y faces, u z = z on n z faces.

This results in
This results in It is seen that C 12 = C 23 = C 13 = K − 2/3G. 7. Computation of C 44 -shear deformation with ε xy = 1 (Fig. 4). BCs are u x = y/2 on n x and n z faces, u y = x/2 on n y and n x faces, u z = 0 on n z faces.

This results in
which confirms G value as it should be. It is worth to note that both prescribed displacements must be applied to all n x and n y faces in order to extort the pure shear state. If u x is applied only on n y faces and u y is applied only of n x faces as in [40] the distorted shape of RVE seems to be correct, but locally near the corners the tension and compression stresses arise and the stress state is far from being the pure shear.   Now we calculate D i j . Here, we also use Voigt's notation as described in [40].  Note that here we have the same value as for D 11 . 12. Computation of D 33 -uniaxial torsional rotations with χ zz = 1. Here Here, we have again the same value as for D 11 and for D 22 . 13. Computation of D 12 -biaxial torsional rotations with χ x x = 1 and χ yy = −1 (Fig. 6). Here, It is worth to note that one of the microrotation gradients (or equivalent torsional moments) should be negative. Otherwise, on the common edges between n x and n y faces the prescribed displacements act in the opposite directions which relates to a certain kinematic inconsistency. 14. Computation of D 23 -biaxial torsional rotations with χ yy = 1 and χ zz = −1 . Here, BCs are u x = yz on n z faces, u y = −zx on n x faces, u z = x y on n x faces, u y = −xz on n z faces. 16. Computation of D 44 -uniform curvature χ xy = 1 (Fig. 7). Here, BCs are u x = xz on n x faces, u y = 0 on n y faces (optionally), Optionally we can also apply u y = 0 on n y faces (Fig. 7). This assumption will increase the magnitude of D 44 BCs are u x = 0 on n x faces (optionally), u y = −yz on n y faces, u z = y 2 2 on n z faces.
18. Computation of D 66 -uniform curvature χ yz = 1. Here, BCs are on n x faces, u y = x y on n y faces, u z = 0 on n z faces (optionally).
19. Computation of D 77 -uniform curvature χ zy = 1. Here, BCs are on n x faces, u y = 0 on n y faces (optionally), u z = −zx on n z faces.
20. Computation of D 88 -uniform curvature χ xz = 1. Here, BCs are u x = −x y on n x faces, u y = x 2 2 on n y faces, u z = 0 on n z faces (optionally).

21.
Computation of D 99 -uniform curvature χ zx = 1. Here, BCs are u x = 0 on n x faces (optionally), u y = − z 2 2 on n y faces, u z = yz on n z faces.
Comparing used boundary conditions with used in [40], one can see that we propose changes in BCs for C 44

Effective properties for a metal foam
Following the benchmark test scheme, we perform here the investigations of mechanical properties of the metallic foam. We consider the representative volume element of dimensions 2.5075 × 2.5075 × 2.5075 mm as shown in Fig. 8. The geometrical model of the foam is obtained by microcomputer tomography with the ScanIP+Fe software [45]. The porosity of the foam (the volume of pores related to the RVE volume) is 93.66%. The isotropic material data of the foam are as follows: Young's modulus is 110 GPa, Poisson's ratio is 0.34, and the mass density is 8960 kg/m3.
The finite element method analysis is made in ABAQUS commercial software. The dense, quite regular mesh consists of tetrahedral elements, see Fig. 9. The size of the problem is large, but reasonable from the point of view of the computation time: the number of nodes is 965,126, the number of elements is 4,077,000, and the total number of unknowns is 2,895,606. The computation time of each problem executed on an average PC equipped with 32 GB memory was about 15 min, which is surprisingly short. It was caused by the effective optimization made by the sparse solver for the analysis of highly porous media. If the solid structure without pores of similar size (about 3 millions of DOFs) is investigated, the execution time would be at least one order longer. Twenty-one numerical cases necessary for the determination of the couple-stress theory material data are the same as in the benchmark problem. In the computation of D i j parameters, the boundary conditions marked in text as optional are not prescribed. It is caused by the fact that the deformed foam does not preserve the flatness of the not loaded surfaces. Below, we present the foam deformations for all considered kinematic loading cases as well as computed strain energy and the value of the corresponding material parameter. The Cartesian coordinate system is displayed only for visualization of the spatial orientation of the RVE. (The global coordinate system is located in the RVE centroid.) The performed calculations are described in the sequel.
1. Computation of C 11 -uniaxial tension with a uniform strain ε x x = 1. The corresponding shape of the RVE is given in Fig. 10. This leads to the following values: U FEM = 3219.78 mJ, C 11 = 2U FEM /V = 408.45 MPa.

Conclusions
Using the linear homogenization technique similar to developed in [40], we have obtained full matrices of elastic moduli which relate the stress and couple stress tensors with the tensors of strain and microcurvature. Here, we have considered a metal foam of porosity 93.66%. We have shown that this foam can be treated as material couple stress. In a certain sense, the presented results solved a problem of reconstruction of constitutive equations for a given microstructured material, which was treated as one of principal problem of continuum mechanics [14,15]. In fact, in addition to other material data these matrices give a source for further modeling of the considered foam and other porous materials in the framework of the couple stress theory, whereas the developed technique can be applied for other materials with complex microstructure. In particular, such beam-lattice structures could be interesting also for non-mechanical applications considering electroelastic coupling, see, e.g., [50,51].