K Du

L Cheng
email: long.cheng@univ-lorraine.fr

J F Barthélémy

I Sevostianov

A Giraud

A Adessina

Numerical computation of compliance contribution tensor of a concave pore embedded in a transversely isotropic matrix

Keywords: Pore inhomogeneity, Transversely isotropic matrix, Supersphere, Superspheroid, Concavity, Contribution and concentration tensors, Numerical homogenization method

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

This paper focuses on the analysis of the effect of a concave pore in a transversely-isotropic material on its overall elastic properties. For this goal we calculate compliance contribution tensors of concave pores (superspherical and superspheroidal) using adapted boundary conditions based Finite Elements Method and evaluate effect of the pore concavity. The work is motivated by multiple experimental observations on irregular character of pores shape in various natural and man-made materials. Particular attention could be payed on the SEM images of [START_REF] Grgic | The influence of CO 2 on the long-term chemo-mechanical behavior of an oolitic limestone[END_REF] showing the concave pores between the calcite grains in the oolitic limestone present (see also [START_REF] Emmanuel | Interfacial energy effects and the evolution of pore size distributions during quartz precipitation in sandstone[END_REF] for sandstones and [START_REF] Wark | Reassessment of pore shapes in microstructurally equilibrated rocks, with implications for permeability of the upper mantle[END_REF] in the case of the harzburgites), which could be described by introducing an concavity-convexity factor proposed by [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]. While the pores and inhomogeneities of irregular shape are typical for materials studied by various branches of natural science, analytical modeling of the properties of materials with microstructures formed by inhomogeneities of non-ellipsoidal shape has not been well developed. The inhomogeneities are typically assumed to be ellipsoids of identical aspect ratios and analytical micromechanical approximations of effective properties are based on the classical Eshelby solution for ellipsoidal inhomogeneities [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | Elastic inclusions and inhomogeneities[END_REF]. This approximation is largely responsible for the huge gap between methods of micromechanics and materials science applications. The reason for this lack is quite obvious: while for 2-D non-elliptical inhomogeneities the technique for evaluation of elastic fields associated with a single inhomogeneity is well developed (see [START_REF] Kachanov | Effective Moduli of Solids With Cavities of Various Shapes[END_REF]; Tsukrov andNovak (2002, 2004); [START_REF] Lanzoni | Effect of pair coalescence of circular pores on the overall elastic properties[END_REF]) only few analytical results are available for non-ellipsoidal 3-D shapes (see discussion in the book [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF]). [START_REF] Argatov | Rigid toroidal inhomogeneity in an elastic medium[END_REF] calculated stiffness contribution tensor of an absolutely rigid thin toroidal inhomogeneity, [START_REF] Kachanov | Rice's internal variables formalism and its implications for the elastic and conductive properties of cracked materials, and for the attempts to relate strength to stiffness[END_REF] obtained analytical solutions for compliance contribution tensors of a crack growing from a pore and a cracks with partial contact between the faces, [START_REF] Krasnitskii | Effect of a rigid toroidal inhomogeneity onthe elastic properties of a composite[END_REF] evaluated elastic fields associated with a rigid torus. Several results have been obtained combining numerical and analytical techniques. Trofimov et al. (2017b) used finite element calculations to analyze the effect of shape of several representative convex polyhedra on the overall elastic properties of particle-reinforced composites. Trofimov et al. (2017a); Trofimov and Sevostianov (2017) quantified effect of waviness of a helical fiber and its elastic properties on the property contribution tensors of such a fiber. In the context of inhomogeneities of concave shape, analytical approximation of compliance contribution tensor of a superspherical pore was first obtained by [START_REF] Sevostianov | On the compliance contribution tensor for a concave superspherical pore[END_REF] using numerical results of [START_REF] Sevostianov | On computation of the compliance and stiffness contribution tensors of non ellipsoidal inhomogeneities[END_REF]. Their result was corrected by [START_REF] Chen | Evaluation of the effective elastic and conductive properties of a material containing concave pores[END_REF] using higher accuracy numerical estimates. Sevostianov et al. (2016a); [START_REF] Chen | Combined effect of pores concavity and aspect ratio on the elastic properties of a porous material[END_REF] considered other types of concave pores. [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF] analyzed elastic fields associated with an inhomogeneity of superspherical pore. These results were used to calculate overall elastic properties of materials with multiple concave pores: oolitic rock [START_REF] Kalo | Effective elastic moduli of a heterogeneous oolitic rock containing 3-d irregularly shaped pores[END_REF]; [START_REF] Giraud | Micromechanical modeling of the effective elastic properties of oolitic limestone[END_REF]) and 3 -D printed S i 3 N 4 ceramics [START_REF] Lurie | Mechanical behavior of porous S i 3 N 4 ceramics manufactured with 3d printing technology[END_REF]. In all of the above mentioned works, it was assumed that the representative elementary volume is an infinite domain (similarly to the Eshelby's hypotheses), which induces in practice, especially for the numerical estimations, some expensive time consuming cost due to the volume size and the mesh refinement. This can be overcame by a recently proposed adapted boundary condition method (Adessina et al. ( 2017)) dedicated to the numerical resolution of the arbitrary shaped inhomogeneity problem. Consequently, the resulting predictions can be found to converge for a relatively small matrix domain and the process is shown to be less time consuming by holding a sufficiently accurate precision. The corrected boundary conditions in this method depend on the elastic properties of the matrix and the method, initially formulated for isotropic matrix only, is extended in the present paper to the case of a transversely isotropic matrix.

This paper is organized as follows. In Section 2, the classical Eshelby problem is reformulated for an inhomogeneity embedded in a finite transversely isotropic matrix by introducing the Green tensor based correction of boundary conditions. It is then applied in Section 3 to the numerical homogenization method for the computation of the contribution tensors. Next, we assess and validate the proposed method in Section 4 by comparing its predictions with the analytical and available numerical results to systematically justify its efficiency and accuracy with respect to the material anisotropy and the pore concavity. The whole procedure leads to some numerical estimations, as presented in Section 4, in the cases of the superspheroidal and superspherical voids planted in the transversely isotropic matrix. Particular attention should be payed to the significant combined effect of the material anisotropy and the shape of pore especially when it is concave. We finally present some concluding remarks in Section 6.

Green tensor based correction of boundary conditions

In Eshelby's footsteps [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], let us consider an infinite domain Ω comprising a matrix surrounding an inhomogeneity E of arbitrary shape. The matrix is linear elastic of homogeneous stiffness tensor denoted by C 0 and the inhomogeneity is also linear elastic but not necessarily homogeneous. Before specifying the work to a transversely isotropic behavior of the matrix embedding a porous concave domain, it is worth noticing that the main results of this section remain theoretically valid in the most general case of anisotropy of the matrix as well as arbitrary shape and internal heterogeneity of E. The infinite domain Ω is submitted to the Hashin-type boundary condition:

ξ(x) ∼ x →∞ E.x (1)
where ξ is the displacement field at the position x and E denotes the remote homogeneous strain. The above mentioned problem is described as:

(P) unbounded                div (σ(x)) = 0 (x ∈ Ω) σ(x) = C(x) : ε(x) (x ∈ Ω) ε = 1 2 gradξ + t gradξ (x ∈ Ω) ξ(x) = E.x (x ∈ ∂Ω) (2)
By introducing the polarization tensor

p = (C(x) -C 0 ) : ε(x) (3)
which is non-zero only in the inhomogeneity E, one has the displacement solution of Eq.( 2) [START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF]; [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF]; [START_REF] Barthélémy | Effective properties of ageing linear viscoelastic media with spheroidal inhomogeneities[END_REF]:

ξ(x) = E.x + x ∈E gradG 0 (x -x ) : p(x ) dΩ x (4)
where G 0 is the second-order Green tensor of the infinite matrix of elasticity C 0 . The first term in the r.h.s. of Eq.( 4) represents the remote displacement field and the second one corresponds to the disturbance caused by the inhomogeneity.

The idea then is to derive from Eq.( 4) a Taylor expansion of the displacement for remote values of x that could be eventually used at relatively large but finite distance (for instance at a mesh boundary in a Finite Element computation) where the sole higher order of the asymptotic behavior (i.e. E.x) may lack of accuracy. For this purpose it can first be noticed that the following approximation is relevant when x

x :

G 0 (x -x ) ∼ x →∞ G 0 (x) ∀x ∈ E (5)
Consequently, Eq.( 4) can be recast into (see also [START_REF] Sevostianov | Elastic fields generated by inhomogeneities: Far-field asymptotics, its shape dependence and relation to the effective elastic properties[END_REF]):

ξ(x) = E.x+ | E | gradG 0 (x) : P (6)
where P defines the average polarisation tensor inside the inhomogeneity which is given as:

P = σ E -C 0 : ε E (7) with ε E = 1 | E | x ∈E ε(x ) dΩ x (8)
being the average strain field of the inhomogeneity.

Correction of boundary condition for a finite domain

We focus henceforward on a finite domain D with a matrix containing an inhomogeneity E. Note once again that the following developments do not require any limitation on the material symmetry of the matrix nor on the shape or content of the inhomogeneity By taking into account Eq.( 6), one has the so called bounded problem expressed as follows:

(P) bounded                div (σ(x)) = 0 (D) σ(x) = C(x) : ε(x) (D) ε = 1 2 gradξ + t gradξ (D) ξ(x) = E.x+ | E | gradG 0 (x) : P (∂D) (9)
It can be considered as the superposition of two elementary linear elastic problems with different boundary conditions, which can be respectively described as:

(P) E bounded                div (σ(x)) = 0 (D) σ(x) = C(x) : ε(x) (D) ε = 1 2 gradξ + t gradξ (D) ξ(x) = E.x (∂D) (10)
and

(P) P bounded                div (σ(x)) = 0 (D) σ(x) = C(x) : ε(x) (D) ε = 1 2 gradξ + t gradξ (D) ξ(x) =| E | gradG 0 (x) : P (∂D) (11)
By separately solving the subproblems (P) E bounded and (P) P bounded , one has:

ε E = A E : E, σ E = B E : E (12)
ε P = A p : P, σ P = B p : P

where A E , B E , A p and B p are the strain and stress concentration tensors in the (P) E bounded and (P) P bounded problems, respectively. Due to their linearity, the solution of the initial problem (P) bounded (i.e. Eq.( 9)) can be obtained by the following superposition:

   ε E = A E : E + A p : P σ E = B E : E + B p : P (14)
Next, by inserting Eq.( 14) into (7), one obtains:

P = D : E (15)
with

D = I -B p + C 0 : A p -1 : B E -C 0 : A E (16)
Consequently, Eq.( 14) can be rewritten as:

   ε E = A E 0 : E, A E 0 = A E + A p : D σ E = B E 0 : E, B E 0 = B E + B p : D (17)
where A E 0 and B E 0 are respectively the average strain and stress concentration tensors of the bounded problem (P) bounded (see also Eq.( 9)).

Determination of elastic compliance and stiffness contribution tensors

As recalled in [START_REF] Kachanov | Effective Moduli of Solids With Cavities of Various Shapes[END_REF]; [START_REF] Sevostianov | Compliance tensor of ellipsoidal inclusion[END_REF]; [START_REF] Sevostianov | Generalization of Maxwell homogenization scheme for elastic material containing inhomogeneities of diverse shape[END_REF], the compliance and stiffness contribution tensors denoted respectively by H E 0 and N E 0 allow to calculate the extra strain and stress induced by the presence of the inhomogeneity in a dilute situation such that

∆ε = f H E 0 : Σ, ∆σ = f N E 0 : E with f = | E | | Ω | ( 18
)
where Σ is the remotely applied stress and as aforementioned, E is the remotely applied strain. Moreover, by applying the consistency laws ensuring that Σ and E are also the average stress and strain within a representative elementary volume, it is useful to notice that the extra stress can be expressed by means of the averages of stress and strain within the inhomogeneity phase in the second term of the r.h.s. of the following decomposition

Σ = C 0 : E + f σ E -C 0 : ε E (19)
Eq.( 19) casts a new light on the definition of the average polarisation tensor in ( 16) and consequently on the stiffness contribution tensor N E 0 which is approximated here by no other than D introduced in Eq.( 15) and related to the partial concentration tensors in Eq.( 16)

N E 0 = D = I -B p + C 0 : A p -1 : B E -C 0 : A E (20)
In the general case of non ellipsoidal shapes, contribution and concentration tensors related to an inhomogeneity need to be calculated numerically as it is presented in this paper (see also in [START_REF] Eroshkin | On micromechanical modeling of particulate composites with inclusions of various shapes[END_REF] details on such calculations). Here whereas the average concentration tensors A E 0 and B E 0 are estimated by Eq.( 17), the contribution ones H E 0 and N E 0 can be interrelated as:

H E 0 = -S 0 : N E 0 : S 0 , N E 0 = -C 0 : H E 0 : C 0 (21)

Case of a transversely isotropic matrix

The practical implementation of the reasoning presented in the previous paragraphs relies on the calculation of the Green tensor and its gradient in the problem of Eq.( 11). In the sequel a transversely isotropic matrix is particularly considered. It is worthy to point out that the three dimensional (3D) solution of G 0 as well as its gradient gradG 0 in the case of transversely isotropy could be analytically calculated by using the results published in literatures [START_REF] Elliott | Three-dimensional stress distributions in hexagonal aeolotropic crystals[END_REF]; [START_REF] Yoo | Elastic interaction of small defects and defect clusters in hexagonal crystals[END_REF]; [START_REF] Pan | Point force solution for an infinite transversely isotropic solid[END_REF]; [START_REF] Mura | Micromechanics of defects in solids[END_REF]; [START_REF] Karapetian | Green's functions for the isotropic or transversely isotropic space containing a circular crack[END_REF]; Pouya (2007a[START_REF] Pouya | Ellipsoidal anisotropy in linear elasticity: Approximation models and analytical solutions[END_REF]; [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF]. However, some of them might present accuracy problems. The misprint and an apparent degenerate case, as quoted by many authors, in Pan and Chou solution [START_REF] Pan | Point force solution for an infinite transversely isotropic solid[END_REF] were corrected by Pouya (2007a[START_REF] Pouya | Ellipsoidal anisotropy in linear elasticity: Approximation models and analytical solutions[END_REF]. It can also be observed that the gradient gradG 0 obtained by [START_REF] Pan | Point force solution for an infinite transversely isotropic solid[END_REF]; Pouya (2007a) are singular on the symmetry axis e 3 . In this paper, the Green functions G 0 i j as well as their gradients G 0 i j,k , are calculated on the one hand by using the solution given in Pouya (2007a) that is also brievely recalled in Appendix B.3, and on the other hand by applying a Fourier transform based integration in order to overcome the singularity problem on the symmetry axis e 3 . The latter is described in Appendix B.2 and validated by comparison with Pouya's solution Pouya (2007a) 1 . Note that the Gaussian integration rule is adopted in the proposed Fourier transform integral based method that is shown to be very efficient and robust in the corresponding numerical implementations and applications. It must be emphasized that the Fourier transform based solution, for the Green tensor and its gradient is valid in the general anisotropic case including all others classes of symmetry.

Case of an homogeneous inhomogeneity

The previous developments, leading to the construction of concentration tensors in Eq.( 17) and contribution tensors in Eq.( 20) and ( 21) from the partial tensors in Eq.( 12) and ( 13), apply to an arbitrary inhomogeneity in terms of shape or content that is possible made of heterogeneous material. However, it may be interesting for practical implementation to examine how these tensors write in the case of an inhomogeneity of uniform stiffness tensor C E . First it is clear that the partial concentration tensors are related by B E = C E : A E and B p = C E : A p . It follows that Eq.( 20) becomes:

N E 0 = D = (C E -C 0 ) -1 -A p -1 : A E (22)
The relationships between the compliance and stiffness contribution tensors as expressed in Eq.( 21) still hold as well as the concentration tensors A E 0 and B E 0 in Eq.( 17) with D = N E 0 given by Eq.( 22). However they can alternatively be written here:

A E 0 = (C E -C 0 ) -1 : N E 0 , B E 0 = C E : A E 0 (23)

Case of an ellipsoidal homogeneous inhomogeneity

The ellipsoidal homogeneous inhomogeneity is of particular interest in the present since analytical expressions of contribution and concentration tensors are available and can then further be compared to the numerical ones to validate the methodology in Section 4.1.In the particular case of an ellipsoidal inhomogeneity E embedded in an infinite matrix 0 of stiffness C 0 and compliance S 0 tensors, compliance H E 0 and stiffness N E 0 contribution tensors write (see [START_REF] Kachanov | Explicit cross-property correlations for porous materials with anisotropic microstructures[END_REF]; [START_REF] Kachanov | Micromechanics of materials, with applications[END_REF] for details) :

H E 0 = (S E -S 0 ) -1 + Q E 0 -1 , N E 0 = (C E -C 0 ) -1 + P E 0 -1 (24)
where P E 0 and Q E 0 denote the fourth order Hill's tensors [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] of the inhomogeneity. Strain concentration tensor of the ellipsoidal inhomogeneity writes

A E 0 = I + P E 0 : (C E -C 0 )
-1

(25)

Numerical framework of the compliance contribution tensor estimate

This section deals with the numerical procedure for the computation of the compliance contribution tensors H E 0 . For the sake of keeping this work focused and concise, we especially consider the inhomogeneity in the case of pore (denoted also by E) that is surrounded by a bounded domain D with a transversely isotropic matrix D M = D\E. Note again that the matrix anisotropy is defined around the axis e 3 in Cartesian system. Different shapes of the pore, in particular the concave ones, will be considered in the present work whose 3D geometries are realized by adopting a user-defined Matlab script. Moreover, as will be detailed in the following part of this section, the bounded domain will be represented either by a spherical model or a cubic one. They are meshed by utilizing the Netgen software [START_REF] Schöberl | Netgen an advancing front 2d /3d-mesh generator based on abstract rules[END_REF] with quadratic 3D elements (C3D10) that are compatible with the Finite Elements computations via Abaqus/Standard software [START_REF] Smith | ABAQUS/Standard User's Manual[END_REF]. Furthermore, the boundary conditions on the external surface ∂D are given by Eq.( 9) that could be decoupled in Eqs.( 10) and ( 11). As a consequence, the numerical procedure will be carried out simultaneously in two parts and through 8 different numerical computations. More specifically,

• for (P) E bounded problem, different boundary strain fields are respectively defined for two simple tension loadings in the directions of e 1 (i.e. E = E 0 e 1 ⊗ e 1 ) and e 3 (i.e. E = E 0 e 3 ⊗ e 3 ), and two simple shear loadings in the planes of e 1e 2 (i.e. E = E 0 (e 1 ⊗ e 2 + e 2 ⊗ e 1 )) and e 1e 3 (i.e. E = E 0 (e 1 ⊗ e 3 + e 3 ⊗ e 1 )) with an arbitrary small constant amplitude E 0 = 10 -5 . The displacement field subjected on the external boundary is calculated by

ξ E i = E i j x j (26)
• for (P) p bounded problem, four analogical computations are realized: two simple tractions with polarization tensor P = P 0 e 1 ⊗ e 1 and P = P 0 e 3 ⊗ e 3 , and two simple shear loadings with P = P 0 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) and P = P 0 (e 1 ⊗ e 3 + e 3 ⊗ e 1 ). A normalized value of P 0 is taken as 1GPa in the numerical computations. The displacements subjected on the external boundary are given by

ξ P i =| E | ∂G i j ∂x k P k j (27)
It is important to emphasize that the proposed method can be applied to any form of the bounded media with transversely isotropic matrix2 . We display in Fig. 1(a) the spherical model comprising the ellipsoidal pore represented by its 1/8 geometry, which is described as

x 1 a 2 + x 2 a 2 + x 3 c 2 = 1 s.t. x 1 ≥ 0, x 2 ≥ 0, x 3 ≥ 0 (28)
with a and c being the major and minor radii, respectively. It can be readily obtained that the aspect ratio γ = c/a. Moreover, the superspherical and axisymmetrical superspheroidal shapes are considered as two promising candidates to comply with the benchmarking representation of the pore concavity, whose geometries are respectively expressed as:

x 1 a 2p + x 2 a 2p + x 3 a 2p = 1 (29) and x 2 1 + x 2 2 a 2 p + x 3 a 2p = 1 ( 30
)
where p is the concavity parameter, a and c are the radii of the principal axes. Note that both the superspherical and axisymmetrical superspheroidal pores are concave when p < 0.5 and convex if p > 0.5. Note that for the sake of efficient mesh generation and convenient numerical homogenization that will be introduced in next sections, we adopt the whole cubic model in the cases of concave pores (see Figs.

Numerical homogenization method

As the considered media being defined in a finite domain (i.e. porous media), the numerical homogenization method is adopted in order to estimate the compliance contribution tensor H. Following Eq.( 17), the strain concentration tensor A E 0 will be first computed. In the framework of homogenization, one has:

ε E D = E = (1 -f ) ε E D M + f ε E E ε P D = (1 -f ) ε P D M + f ε P E (31)
where f =| E | / | D | is the porosity, ε E D and ε P D denote the subjected macroscopic strain respectively for the (P) E bounded and (P) P bounded problems, ε E E , ε P E , ε E D M and ε P D M are the corresponding average strains in the porous phase and those in the matrix phase. Following the Gauss theorem, ε P D M can be computed as an integral over the external boundary ∂D:

ε i j P D = 1 2 | D | M m=1 ξ P i n j + n i ξ P j (m) S (m) ( 32
)
where M is the total number of the surface elements and S (m) is the area of the m-th one. ε i j E D M and ε i j P D M can be obtained, respectively for the (P) E bounded and (P) P bounded problems, by averaging the matrix strains of the elements as:

ε i j D M = 1 | D M | N n=1 ε i j V n (33)
with N being the total number of the volume elements and V (m) giving the volume of the n-th one.

Having in hand the above computed quantities, the average strain field in the porous phase can be obtained as:

ε E E = E -(1 -f ) ε E D M f ε P E = ε P D -(1 -f ) ε P D M f (34)

Strain concentration tensor and compliance contribution tensor

In this section, the strain concentration tensor A E 0 and the compliance contribution one H will be calculated based on the numerical computation as described in Section 3.1. To this end, we aim first at computing the A E and A p that both have 7 independent non-zero components. More specifically, for each of them, the components denoted as A 1111 , A 3333 , A 1122 , A 1133 and A 3311 are calculated from 2 tension and traction loadings and those of A 1212 and A 1313 can be obtained from 2 shear ones. Note that the particular consideration of the pore inhomogeneity leads to the fact that the stress concentration tensors B E and B p both vanished. Consequently, the stiffness contribution tensor N E 0 (see also Eq.( 22)) and the compliance contribution tensor H E 0 can be simplified as:

N E 0 = -S 0 + A p -1 : A E , H E 0 = S 0 : S 0 + A p -1 : A E : S 0 (35)

Assessment and validation of the proposed numerical procedure

In this section, the estimations of compliance contribution tensor H obtained from the proposed numerical procedure will be assessed and validated by comparison with the avalaible analytical and numerical results in literatures. In this light, we systematically consider the spheroidal pores embedded in a transversely isotropic matrix as well as the superspherical ones surrounded by an isotropic matrix to respectively justify its accuracy on the anisotropy of the matrix material and that on the concavity of the pores. More specifically, in the spheroid case, comparison will be made with respect to the analytical results proposed by [START_REF] Mura | Micromechanics of defects in solids[END_REF]. Whereas in the concave superspherical cases, due to the lack of analytical results, we compare the numerical predictions with the available FEM simulations that were recently obtained in [START_REF] Chen | Accuracy of the replacement relations for materials with non-ellipsoidal inhomogeneities[END_REF]; Sevostianov et al. (2016b) (see also [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]). The transverse isotropy of the elastic matrix material is described by the parameters as shown in Table1. We consider a spherically bounded model comprising an ellipsoidal pore, due to its symmetry, whose 1/8 geometry and mesh are illustrated in Fig. 2. The corrected boundary conditions ( 9) and the uncorrected ones (10) will be respectively subjected to the same mesh to study the efficiency and the accuracy of the proposed numerical procedure. Note that in order to evaluate the efficiency of the numerical method with corrected boundary conditions, the computations will be carried out by adopting different spherical models with a fixed aspect ratio γ by varying the the scale ratio a/L between the major radius of the ellipsoid a and the the radius of the spherical model L. Each independent components of the compliance contribution tensor H i jkl will be numerically computed and compared with the analytical solution [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF].

Properties of material

E 1 (GPa) E 3 (GPa) ν 12 ν 31 G 13 (
Figure 2: Geometry and mesh of 1/8 spherical inclusion in the center of spherical matrix Fig. 3 illustrates the numerical estimation for each independent component of H tensor in the case of the aspect ratio γ = 1/5. A striking observation from this figure is that in the case of transversely isotropic host matrix, the numerical predictions obtained by switching on the correction of the boundary conditions converge more quickly than those obtained from the classical modeling. The convergence between them can observed when a/L 8 and that with corrections of boundary conditions is shown to be more efficient. Additionarilly, we show in Fig. 4 the relative errors of the numerical estimations with respect to the analytical solution [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] (recalled in appendix Appendix B). First of all, it can be readily observed that when a/L is approximately in the range of [4,8], the relative errors obtained from the corrected boundary conditions are shown to be very small (as around the value of 10 -4 ), whereas those obtained from the classical modeling is relatively significant. It is logical and evident that the corrections of the boundary conditions allows to accelerate the numerical convergence without degrading the computation accuracy. It may be noticed that non monotonous evolutions of relative errors as function of ratio a/L may be observed for small values of a/L and components (H 3333 , H 1122 , H 1133 ), in the case of corrected boundary condition. In the paper [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF], this aspect of the curve is assumed to be probably due to the error compensations, and is interpreted as an compensation between the effect of the mesh refinement and that of the boundary condition correction. By contrast, when a/L ∈ [8, 18] , the relative errors asymptotically converge and attaining a value around 10 -3 . This can be interpreted as the fact that when the scale ratio a/L is sufficiently big, the representative bounded model tends to be an infinite one such that the correction of the boundary conditions is hence neither efficient nor useful in the numerical modeling. For completeness, we provide in Appendix D the computation results of the strain contribution tensor A in the case of γ = 1/5 as well as some supplementary ones of the compliance contribution tensor H for some other ellipsoidal pores such as γ = 1 (i.e. spherical pore), γ = 1/2 and γ = 1/10. By considering all of the above mentioned computation results, the scale ratio a/L = 8 is supposed to be fixed in the next part of this work. This is of course an approximation but will be shown as sufficiently accurate in the following numerical estimations even in the case of cubically bounded matrix comprising a concave pore.

Superspherical pore in the cubically bounded model

This section deals with the assessment and validation of the numerical procedure by paying particular attention to the concavity of the pore shape. In this light, we switch off the matrix anisotropy 3 and consider that the pore is in a superspherical form (see also Eq.( 29)). As aforementioned, the bounded domain is represented by a cubic model with the scale ratio a/L = 8 for the convenient mesh generation and numerical homogenization. Different geometries and meshes are realized by varying the value of the concavity p. We show in Fig. 5 an example in the case of p = 0.4 and its mesh that is sufficiently refined in the transition zone between the matrix and porous phases. Moreover, the bounded cubic model is subjected only to the corrected boundary conditions. The predictions are compared with the FEM results that have been published in [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF]. In Table .2, we display several numerical estimations computed with different values of the concavity p as well as the comparisons with the available results obtained by [START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF] via the classical numerical modeling without corrections of boundary conditions. A very good agreement can be found and the relative errors are shown to be minor especially by paying attention to the cases of p ≤ 1. In turn, it justifies that the accuracy of the proposed numerical method is unaffected by the pore concavity.

Numerical estimation in the case of concave pores

In this section, we propose to carry out the study of the compliance contribution tensor H in the case of transversely isotropic matrix comprising concave pore to understand in more detail their combined effect. The superspheroidal and superspherical shapes of the pores will be respectively considered by varying the concavity parameter p in a relatively large interval such as [0.2, 5]. Again, we restrict the study, particularly in the superspheroidal case (i.e. shape verifying the symmetry of revolution), to the assumption that the directions of the symmetry between the matrix anisotropy and that of the pore are aligned on the same direction. In Figs. 6 and5 , we show the cubic geometries that comprising the superspheroidal and superspherical pores as well as the corresponding surface meshes on internal boundaries. Moreover, for the sake of prediction accuracy, the numerical computation will be carried out based on the sufficiently refined meshes, for which the number of nodes and elements is detailed in Appendix C. Note that we keep the material parameters fixed as previously introduced in Table .1.

Superspheroidal pores

We provide in Table .4 the numerical estimations of each independent non-zero component of H i jkl for the superspheroidal pores with different values of concavity parameter p. Since the transverse isotropy of the matrix material and that of the superspheroidal pore are both around the e 3 , we can theoretically conclude that the overall response such as the compliance contribution tensor must satisfy the equality (H 1111 -H 1122 )/2 = H 1212 . This property is well verified from the computation results that indirectly demonstrate their accuracy. Fig. 7 illustrates the H i jkl components as a function of the concavity parameter p. It can be observed that the evolutions, in particular those of H 3333 and H 1313 , are more important when p < 0.5. In other words, the concave form of the porous heterogeneity has a significant effect on the compliance contribution tensor H. Whereas it becomes qualitatively stabled when p > 0.5 (i.e. convex forms).

Superspherical pores

We investigate the numerical estimation of H i jkl for the superspherical pores surrounded by a cubically bounded matrix whose geometry and mesh is illustrated in Fig. 5(b). Detailed results for different values of concavity p ∈ [0.2, 5] are summarized in Table .4. First of all, it should be emphasized that the supersphere, as described in Eq.( 29), is obviously not transverse isotropic. As a consequence, the aforementioned equality in the case of superspheroidal pores is no longer satisfied. Fig. 8 shows the evolution of the H i jkl components on the concavity parameter p. As same as the superspheroidal cases, it can also be finally concluded that the concavity of the superspherical pore significantly affect its compliance contribution tensor H. For completeness, the numerical evaluation of the strain concentration tensor A are reported in Appendix E.

Concluding remarks

In the present work, we have numerically evaluated the compliance contribution tensor of the concave pore inhomogeneity embedded in a transversely isotropic matrix. This has been realized by use of an original developed numerical homogenization method complying with the adapted boundary conditions based method recently formulated by [START_REF] Adessina | Effective elastic properties of materials with inclusions of complex structure[END_REF]. The proposed numerical procedure was carried out for an arbitrarily bounded representative elementary volume and is shown to be efficient and accurate in the numerical modeling. By paying particular attentions to the pore concavity and the matrix anisotropy, a major contribution of this work is found as the sufficiently exact computation results and analysis that illustrate the significant effect of the pore concavity on the elastic properties of the matrix.

Specifically, the adapted boundary conditions based method was extended in the context of the matrix anisotropy thanks to the Green functions and their gradients applied in the correction of the boundary conditions, which have been rigorously reformulated via the Fourier transform based integral method in particular by solving the singularity problem on the axis of the transverse isotropy. Moreover, the numerical homogenization method has been utilized in the proposed numerical procedure that has been firstly assessed and validated by comparing its predictions with the analytical and existing numerical results in particular cases. It is then used to the investigate the effect of the concave pore on the transversely isotropic matrix from the quantitative estimates of the compliance contribution tensor, which was found to be of critical importance especially in the case of the concavity being less than 0.5.

Last but not least, the proposed numerical method is able to deal with any general anisotropy of the matrix material but such a study has not been attempted here for the sake of keeping the work focused and concise. In the perspective point of view, effective properties such as those predicted from semi-analytical homogenization models could be developed based on the obtained numerical computations. Alternatively, new formulations of the "adapted boundary conditions" as well as their applications in the context of non-linear homogenization may also constitute a challenging extension.

Appendix B.2. Exact results on symmetry axis in transversely isotropic case obtained by using Fourier transform solution In the particular case of transversely isotropic matrix stiffness tensor with symmetry axis x 3 , see relations (A.7-A.8-A.9)

C 0 = 6 i=1 c i E i (B.8)
One considers a position vector located on symmetry axis

x = x 3 e 3 , x = r = |x 3 |, v = x x = e 3 , u ψ = cos ψe 1 + sin ψe 2 (B.9)
Inverse of acoustic tensor writes

K -1 =      c 2 (1-cos(2ψ))+2c 5 c 5 (c 2 +c 5 ) -c 2 sin(2ψ) c 5 (c 2 +c 5 ) 0 -c 2 sin(2ψ) c 5 (c 2 +c 5 ) c 2 (1+cos(2ψ))+2c 5 c 5 (c 2 +c 5 ) 0 0 0 2 c 6      (B.10)
and corresponding Green function

G 0 i j (x = x 3 e 3 ) = 1 8 π 2 x 2 π 0 K -1 i j (u ψ ) dψ (B.11)
Non zero components write

G 0 11 (x = x 3 e 3 ) = G 0 22 (x = x 3 e 3 ) = 1 4 π x c 2 + 2c 5 c 5 (c 2 + c 5 ) , G 0 33 (x = x 3 e 3 ) = 1 2 π x c 6 , x = |x 3 | (B.12)
and it coincides with relations (B.33-B.34).

K -1 . v.C 0 .u ψ + u ψ .C 0 .v .K -1 = 2 √ 2c 3 + c 6 (c 2 + c 5 ) c 6      0 0 cos(ψ) 0 0 sin(ψ) cos(ψ) sin(ψ) 0      (B.13) G 0 i j,k (x = x 3 e 3 ) = 1 8 π 2 x 2 2 π 0 K -1 . v.C 0 .u ψ + u ψ .C 0 .v .K -1 i j u ψ k -K -1 i j v k dψ, x = |x 3 | (B.14)
Non zero components write

G 0 11,3 (x = x 3 e 3 ) = G 0 22,3 (x = x 3 e 3 ) = - 1 4 π x 2 c 2 + 2c 5 c 5 (c 2 + c 5 ) , G 0 33,3 (x = x 3 e 3 ) = - 1 2 π x 2 c 6 (B.15) G 0 32,2 (x = x 3 e 3 ) = G 0 31,1 (x = x 3 e 3 ) = G 0 13,1 (x = x 3 e 3 ) = G 0 23,2 (x = x 3 e 3 ) = 1 4 π x 2 √ 2c 3 + c 6 (c 2 + c 5 ) c 6 (B.16)
or, in terms of usual components

G 0 11 (x = x 3 e 3 ) = 1 8 π x 3C 0 1111 -C 0 1122 C 0 1111 C 0 1111 -C 0 1122 , G 0 33 (x = x 3 e 3 ) = 1 4 π x C 0 2323 (B.17) G 0 11,3 (x = x 3 e 3 ) = - 1 8 π x 2 3C 0 1111 -C 0 1122 C 0 1111 C 0 1111 -C 0 1122 , G 0 33,3 (x = x 3 e 3 ) = - 1 4 π x 2 C 0 2323 (B.18) G 0 23,2 (x = x 3 e 3 ) = 1 8 π x 2 C 0 1133 + C 0 2323 C 0 1111 C 0 2323 (B.19)
Appendix B.3. Exact 3D elastic Green function in the transversely isotropic case We briefly recall exact Green function given in Pouya (2007b). We only consider non degenerate case c-c 3 √ 2 -c 6 0 with n = e 3 (n denotes unit vector on symmetry axis of transversely isotropic material). See Pouya (2007b) for discussion on more general and non degenerate cases.

c = c 1 (c 2 + c 5 ) 2 (B.20) υ 1 =   c -c 3 √ 2 c + c 3 √ 2 + c 6 2c 1 c 6   1/2 +   c + c 3 √ 2 c -c 3 √ 2 -c 6 2c 1 c 6   1/2 (B.21) υ 1 =   c -c 3 √ 2 c + c 3 √ 2 + c 6 2c 1 c 6   1/2 -   c + c 3 √ 2 c -c 3 √ 2 -c 6 2c 1 c 6   1/2 (B.22) υ 3 = c 5 c 6 1/2 , υ 4 = c 2 + c 5 √ 2c 3 + 2c 6 1/2 (B.23)
In the particular case investigated

ζ = x.n = x 3 , ρ = x.x -ζ 2 = x 2 1 + x 2 2 (B.24)
One defines 6 functions R α and R * α (α = 1, 2, 3)

R α = ρ 2 + υ 2 α ζ 2 , R * α = R α + υ α ζ (B.25)
In the non degenerate case c -c 3 √ 2c 6 0

w 1 = 2 α=1 A α R * α , w 2 = 2 α=1 - A α R α R * 2 α , w 3 = 2 α=1 - υ 2 α A α ρ 2 R α , w 4 = 2 α=1 c 2 + c 5 -c 6 υ 2 α √ 2c 3 + c 6 A α R α (B.26) A α = (-1) α c 6 - √ 2c 3 υ 2 α 2c 1 υ 2 1 -υ 2 2 υ α , A α = (-1) α √ 2c 3 + c 6 2c 1 υ 2 1 -υ 2 2 υ α , α = 1, 2 (B.27)
As indicated in Pouya (2007b), the singularity on x 3 axis (ρ = 0) may be removed by using the following expression for w 3

w 3 = - √ 2c 3 + c 6 2c 1 R 1 R 2 (υ 2 R 1 + υ 1 R 2 ) (B.28)
Green function's for the infinite transversely isotropic medium writes

G 0 i j (x) = 1 2 π c 6 g 1 δ i j + g 2 x i x j + g 3 x i δ 3 j + x j δ 3i + g 4 δ 3i δ 3 j (B.29)
with Calculation of gradient G 0 i j,k (x) is straightforward by using formal calculation tools.

g 1 = w 1 + ζ R 3 R * 3 , g 2 = w 2 + 1 υ 3 R 3 R * 2 3 , g 3 = (w 3 -g 2 ) ζ, g 4 = w 4 -g 1 -g 2 ζ 2 -2g 3 ζ (B.30) Table B.5: G 0 i jk : numerical results at (x 1 = -1, x 2 = 0.8, x 3 = 1.5) k = 1 k = 2 k = 3
G 0 i j,k (x) = 1 2 π c 6 g 1,k δ i j + g 2,k x i x j + g 3,k x i δ 3 j + x j δ 3i + g 4,k δ 3i δ 3 j + g 2 δ ik x j + δ jk x i + g 3 δ ik δ 3 j + δ jk δ 3i (B.31)
Considering the data given in table (1) and position vector (x 1 = -1, x 2 = 0.8, x 3 = 1.5) one obtains 39800.47 -903.6818 -7004.161 -903.6818 39393.815 5603.329 -7004.161 5603.329S 134599.77

10 7 G 0 =     
     (B.32)
We recall that Green function is not singular on the axis x 3 , except at the origin x 3 = 0. In the case ρ → 0 and x 3 0, it writes ρ=0) 11 δ 1i δ 1 j + δ 2i δ 2 j + g (ρ=0) 33 δ 3i δ 3 j (B.33)

G 0 i j x 3 e 3 = 1 2 π c 6 x g (
g (ρ=0) 11 = 1 2 A 1 υ 1 + A 2 υ 2 + 1 υ 2 3 , g (ρ=0) 33 = A 1 υ 1 c 2 + c 5 -c 6 υ 2 1 √ 2c 3 + c 6 + A 2 υ 2 c 2 + c 5 -c 6 υ 2 2 √ 2c 3 + c 6 (B.34)
Appendix B.4. Strain Hill polarization tensor of a spheroidal inclusion aligned in a transversely isotropic host matrix We only recall in this section the solution of strain Hill polarization tensor of a spheroidal inclusion aligned in a transversely isotropic host matrix, in the non degenerate case (c -c 3 √ 2c 6 0) with n = e 3 . We use constants c, υ i defined in section (Appendix B.3) and we present Withers solution [START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] as in [START_REF] Parnell | The Eshelby, Hill, Moment and Concentration Tensors for Ellipsoidal Inhomogeneities in the Newtonian Potential Problem and Linear Elastostatics[END_REF]. See also [START_REF] Laws | A note on penny-shaped cracks in transversely isotropic materials[END_REF]; [START_REF] Sevostianov | Effective elastic properties of matrix composites with transversely-isotropic phases[END_REF] for exact solutions derived from Fourier transfom integral. The function S(x) (with lim x→1 S(x) = 1 3 ) characterizes the influence of the shape of the spheroid)

S(x) = 1 1 -x 2 - x 1 -x 2    1 √ 1-x 2 arccos (x) , x < 1 1 √ x 2 -1 arcosh (x) , x > 1 (B.35)
Functions I 1 , I 3 (a misprint in relation 5. 124 Parnell (2016) has been fixed) write:

I 3 (υ i ) = 4π υ i S(υ i γ), I 1 (υ i ) = 2π υ i - I 3 (υ i ) 2 = 2π υ i (1 -S(υ i γ)) (B.36)
D = 1 2 πc 6 υ 3 , L i = (c 2 + c 5 ) /υ 2 i -c 6 √ 2c 3 + c 6 , M i = (-1) i c 6 -2c 1 υ 2 i 8 π c 1 c 6 υ 2 1 -υ 2 2 υ 2 i (B.37)
Components p i of strain Hill tensor P E 0 in TI tensor basis write Appendix C. Information of meshes for the FEM computations in the case of cubic model containing concave pore

p 1 = -2 2 i=1 L 2 i M i υ 5 i I 3 (υ i ), p 2 = 2 2 i=1 M i υ i I 1 (υ i ), p 3 = √ 2 2 i=1 L i M i υ 3 i I 3 (υ i ), p 4 = p 3 (B.38) p 5 = 2 i=1 M i υ i I 1 (υ i ) + D I 1 (υ 3 ) 2 , p 6 = 1 2 2 i=1 (1 + L i )M i υ 3 i (I 3 (υ i ) -2L i I 1 (υ i )) + Dυ 2 3 I 3 (υ 3 ) 4 (B.
We provide in Tables C.8 and C.9 the mesh information during the FEM computation for cubically bounded representative elementary volume containing respectively the superspheroidal and the superspherical pore. The displayed numbers of nodes and those of 3D quadratic elements shows that each mesh is well refined for the corresponding computations to obtain a precision of computation as accurate as possible. (i.e. spherical pore), γ = 1/2 and 1/10.

Appendix D.1. Strain concentration tensors in the spherical case and spheroidal one with γ = 1/5 Fig.D.9 displays the numerical estimations of each independent component of the strain concentration tensor A for the ellipsoidal pore with the aspect ratio γ = 1/5. As can be found in the computation results of H, similar evolution trends in function of the scale ratio a/L, can be observed. It should be note here that, unlike the compliance contribution tensor H, the strain concentration tensor A does not present the major symmetry between the components A 1133 and A 3311 . As a consequence, it has 7 non-zero independent components.

Additionally, we show in Fig. D.10 the relative errors of their numerical estimations with respect to the analytcial solutions proposed by [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] as well as those obtained from the classical numerical computations without the correction of the boundary conditions. A very good precision can be found for the proposed numerical procedure in the case of smaller scale ratios a/L ∈ [4, 8], whereas the correction of the boundary conditions is shown to be useless when a/L is sufficiently important. In Figs.D.11 -D.16, we show some supplementary results of the compliance contribution tensor H for the ellipsoidal pores respectively with the aspect ratio γ = 1, 1/2 and 1/10. By simultaneously considering those of γ = 1/5 (see Section 4.1), we focus on the value of a/L at the starting point of the convergence between the related evolution respectively obtained from the proposed numerical procedure and the classical one. Finally, we fix a/L = 8 for the numerical computations in the cases of the concave pores. .9275 -0.2891 -0.2936 15.9062 47.1787 2.2130 5.3414 0.25 2.0047 -0.2903 -0.1911 6.3536 20.3159 2.2940 2.6196 0.3 2.1128 -0.2966 -0.1204 3.3695 11.7872

Figure 1 :

 1 Figure 1: Geometries of different representative models: (a) the spherical model with ellipsoidal pore; (b) the cubic model comprising the superspherical pore; (c):the cubic model embedded with superspheroidal pore.

Figure 3 :

 3 Figure 3: Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratio a/L ∈ [4, 18].

Figure 4 :

 4 Figure 4: Relative errors of numerical computations of H i jkl with respect to the analytical results[START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratios a/L ∈ [4, 18].

Figure 5 :

 5 Figure 5: Geometry and mesh for a cubic model comprising a superspherical pore with p = 0.4.

Figure 6 :

 6 Figure 6: Geometry and mesh of cubic model comprising a superspheroidal pore.

Figure 7 :

 7 Figure 7: Evolution of components H i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5]

Figure 8 :

 8 Figure 8: Evolution of components H i jkl for the superspherical pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5]

Figure D. 9 :

 9 Figure D.9: Numerical estimations of A i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/5 and different scale ratio a/L ∈ [4, 18].

Figure D. 15 :

 15 Figure D.15: Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/10 and different scale ratio a/L ∈ [4, 18].

Table 2 :

 2 Numerical estimations for the independent components of H tensor in the case of isotropic host matrix embedded with a superspherical pore and the comparison with the results obtained in[START_REF] Trofimov | Overall elastic properties of a material containing inhomogeneities of concave shape[END_REF].

Table 3 :

 3

	p	H FEM 1111	H FEM 1122	H FEM 1133	H FEM 3333	H FEM 1212	H FEM 1313	H FEM 1111 -H FEM 1122 2	Error 1212	1
	0.2 0.1004 -0.0192 -0.0520 3.6670 0.0597 0.8424	0.0598	0.16%
	0.25 0.1026 -0.0212 -0.0442 1.5948 0.0619 0.4132	0.0619	0.04%
	0.3 0.1068 -0.0232 -0.0395 0.9354 0.0650 0.2928	0.0650	0.01%
	0.35 0.1102 -0.0248 -0.0359 0.6445 0.0675 0.2495	0.0675	0.00%
	0.4 0.1128 -0.0260 -0.0331 0.4929 0.0694 0.2324	0.0694	0.00%
	0.45 0.1147 -0.0270 -0.0310 0.4047 0.0708 0.2248	0.0708	0.00%
	0.5 0.1159 -0.0276 -0.0293 0.3492 0.0718 0.2211	0.0718	0.00%
	0.6 0.1175 -0.0284 -0.0271 0.2900 0.0730 0.2187	0.0730	0.01%
	0.7 0.1182 -0.0288 -0.0258 0.2604 0.0735 0.2188	0.0735	0.01%
	0.8 0.1186 -0.0291 -0.0249 0.2436 0.0738 0.2197	0.0738	0.01%
	0.9 0.1187 -0.0292 -0.0243 0.2333 0.0739 0.2211	0.0739	0.01%
	1	0.1187 -0.0293 -0.0239 0.2265 0.0740 0.2226	0.0740	0.01%
	1.5 0.1184 -0.0293 -0.0225 0.2129 0.0739 0.2292	0.0739	0.01%
	2	0.1182 -0.0294 -0.0219 0.2098 0.0738 0.2343	0.0738	0.02%
	2.5 0.1181 -0.0294 -0.0214 0.2091 0.0737 0.2381	0.0738	0.01%
	3	0.1180 -0.0295 -0.0211 0.2092 0.0738 0.2411	0.0738	0.02%
	4	0.1180 -0.0296 -0.0207 0.2099 0.0738 0.2453	0.0738	0.02%
	5	0.1181 -0.0297 -0.0204 0.2108 0.0739 0.2481	0.0739	0.02%
	1 Relative error of	H FEM 1111 -H FEM 1122 2	with respect to H FEM 1212

Numerical estimation of H i jkl for the superspheroidal pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5].

Table 4 :

 4 Numerical estimation of H i jkl for the superspheroidal pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5].

	p	H FEM 1111	H FEM 3333	H FEM 1122	H FEM 1133	H FEM 1212	H FEM 1313
	0.2 0.5401 0.9682 -0.1386 -0.1093 0.2043 0.4972
	0.25 0.3391 0.6106 -0.0895 -0.0676 0.1288 0.3327
	0.3 0.2521 0.4567 -0.0682 -0.0496 0.0984 0.2715
	0.35 0.2052 0.3749 -0.0567 -0.0400 0.0844 0.2454
	0.4 0.1770 0.3269 -0.0497 -0.0343 0.0779 0.2330
	0.45 0.1589 0.2967 -0.0450 -0.0308 0.0749 0.2265
	0.5 0.1460 0.2747 -0.0413 -0.0284 0.0735 0.2229
	0.6 0.1337 0.2542 -0.0371 -0.0263 0.0727 0.2199
	0.7 0.1269 0.2419 -0.0342 -0.0252 0.0728 0.2194
	0.8 0.1229 0.2345 -0.0321 -0.0246 0.0731 0.2200
	0.9 0.1204 0.2297 -0.0305 -0.0242 0.0735 0.2212
	1	0.1188 0.2266 -0.0293 -0.0239 0.0740 0.2226
	1.5 0.1155 0.2202 -0.0257 -0.0228 0.0764 0.2299
	2	0.1149 0.2191 -0.0239 -0.0222 0.0783 0.2359
	2.5 0.1150 0.2192 -0.0228 -0.0218 0.0799 0.2406
	3	0.1152 0.2196 -0.0221 -0.0214 0.0811 0.2442
	4	0.1158 0.2207 -0.0211 -0.0209 0.0829 0.2495
	5	0.1162 0.2215 -0.0205 -0.0206 0.0841 0.2530

Table B .

 B 6: G 0 i j and G 0 i jk : numerical results at (x 1 = 0, x 2 = 0, x 3 = 1.)

		k = 1	k = 2	k = 3
	10 7 g kk (no summation) 60722.101 60722.101 502034.39
	10 7 g 11k	0.	0.	-60722.101
	10 7 g 22k	0.	0.	-60722.101
	10 7 g 33k	0.	0.	-502034.39
	10 7 g 23k	0.	70648.848	0.
	10 7 g 31k	70648.848	0.	0.
	10 7 g 12k	0.	0.	0.

Table B .

 B 7: Reference transversely isotropic elastic parameters: C 0 i components in TI tensor basis and components of Hill tensor related to an oblate spheroid of aspect ratio γ = 0.5

		c 0 1111 (GPa)	C 0 3333 (GPa)	C 0 1122 (GPa)	C 0 1133 (GPa)	C 0 2323 (GPa)
		22.3639	12.9994	3.8275		4.7092	1.5851
		c 1 (GPa)	c 2 (GPa)	c 3 = c 4 (GPa)	c 5 (GPa)	c 6 (GPa)
		12.9994	26.1914	6.65983		18.5363	3.1702
		p 1 (GPa -1 )	p 2 (GPa -1 )	p 3 = p 4 (GPa -1 )	p 5 (GPa -1 )	p 6 (GPa -1 )
		0.0656356	0.0199043	-0.0100943	0.0292588	0.106492
		p 1111 (GPa -1 ) p 3333 (GPa -1 )	p 1122 (GPa -1 )	p 1133 (GPa -1 ) p 2323 (GPa -1 )
		0.0245816	0.0656356	-0.00467723	-0.00713773	0.0532462
	p	0.2	0.25	0.3	0.35	0.4	0.45	0.5	0.6	0.7
	Num. N. 1 5580080 5062346 4246894 3412606 3777016 1965376 1272710 2049212 2125594
	Num. E. 2 3949744 3243744 2958504 2289744 2575736 1377480 843280 1361640 1419840
	p	0.8	0.9	1	1.5	2	2.5	3	4	5
	Num. N. 1047473 1716790 798157 2284278 2724190 2072024 2233528 2422684 2668330
	Num. E. 1177376 1146672 545358 1531896 1823384 1395296 1501464 1629936 1801512
	1 Number of nodes
	2 Number of elements
		Table C.8: Number of nodes and elements in the meshes of the cubic models comprising different superspheroid pores
	p	0.2	0.25	0.3	0.35	0.4	0.45	0.5	0.6	0.7
	Num. N. 4829072 2859018 3539342 3025846 3087712 2166076 1029954 1893940 2447226
	Num. E. 3260176 1924360 2473320 2020488 2049608 1437104 762880 1309424 1716608
	p	0.8	0.9	1	1.5	2	2.5	3	4	5
	Num. N 2487370 2288388 682512 1079731 2769712 2789860 3094708 1281679 749678
	Num. E. 1716568 1579872 481944	739903 1866728 1883248 2091464 889557	505988

Table C .

 C 9: Number of nodes and elements in the meshes of the cubic models comprising different superspherical pores

  Relative errors of numerical computations of H i jkl with respect to the analytical results[START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the spherical pore γ = 1 embedded in a transversely isotropic matrix with different scale ratios a/L ∈ [4, 18]. Numerical estimations of H i jkl for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/2 and different scale ratio a/L ∈ [4, 18]. Relative errors of numerical computations of H i jkl with respect to the analytical results[START_REF] Withers | The determination of the elastic field of an ellipsoidal inclusion in a transversely isotropic medium, and its relevance to composite materials[END_REF] for the ellipsoidal pore embedded in a transversely isotropic matrix with fixed aspect ratio γ = 1/2 and different scale ratios a/L ∈ [4, 18].

		0.12 10 0 0.105 10 0 0.0634		0.23 10 0 0.45 10 0 1.64
		0.1
		0.118 0.06335		10 -1 0.4 1.62
	Relative error Relative error	0.114 0.116 10 -4 10 -2 0.08 0.085 0.09 0.095 10 -4 10 -2 0.0632 0.06325 0.0633	Relative error Relative error	0.22 0.225 10 -3 10 -2 0.3 0.35 10 -4 10 -2 1.58 1.6
	• Spherical pore: γ = 1 4 6 8 0.112 (a) H 1111 10 0.075 10 -6 0.07 0.06315 4 6 8 10 -0.03 -0.029 -0.028 -0.027 -0.026 -0.025 (c) H 4 6 8 10 0.069 0.07 0.071 0.072 0.073 0.074 0.075 (e) H 1212 Figure D.11: Numerical estimations of H i jkl for the spherical pore γ = 1 embedded in a transversely isotropic matrix with different scale ratio 12 14 16 18 4 6 8 10 12 14 16 0.215 10 -4 0.25 1.56 18 (b) H 3333 12 14 16 18 4 6 8 10 12 14 16 18 -0.024 -0.0235 -0.023 -0.0225 (d) H 1133 12 14 16 18 4 6 8 10 12 14 16 18 0.217 0.218 0.219 0.22 0.221 0.222 0.223 0.224 (f) H 1313 a/L ∈ [4, 18]. 4 6 8 10 12 14 16 18 (a) H 1212 4 6 8 10 12 14 16 18 (b) H 1313 Figure D.12: 2 4 6 8 10 12 14 16 18 (a) H 1111 2 4 6 8 10 12 14 16 18 (b) H 3333 2 4 6 8 10 12 14 16 18 -0.025 -0.02 -0.015 -0.01 (c) H 1122 2 4 6 8 10 12 14 16 18 -0.034 -0.032 -0.03 -0.028 -0.026 -0.024 -0.022 -0.02 (d) H 1133 2 4 6 8 10 12 14 16 18 0.04 0.045 0.05 0.055 0.06 0.065 (e) H 1212 2 4 6 8 10 12 14 16 18 0.2 0.21 0.22 0.23 0.24 0.25 (f) H 1313 4 6 8 10 12 14 16 18 (a) 2 4 6 8 10 12 14 16 18 (b) Figure D.14: 4 6 8 10 12 14 16 18 (a) H 1111 4 6 8 10 12 14 16 18 (b) H 3333 4 6 8 10 12 14 16 18 -6.96 -6.95 -6.94 -6.93 -6.92 -6.91 -6.9 10 -3 (c) H 1122 4 6 8 10 12 14 16 18 -0.0365 -0.036 -0.0355 -0.035 (d) H 1133 4 6 8 10 12 14 16 18 0.03505 0.0351 0.03515 0.0352 4 6 8 10 12 14 16 18 0.493 0.4935 0.494 0.4945 0.495 0.4955 0.496 0.4965 Figure D.13: 2 (e) H 1212 (f) H 1313

Table E .

 E 10: Evolution of components A i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the Logarithm value of the concavity parameter log(p) such that p ∈ [0.2, 5] Figure E.17: Evolution of components A i jkl for the superspheroidal pore embedded in a transversely isotropic matrix on the concavity parameter p such that p ∈ [0.2, 5]

	p	50 A FEM 1111	A FEM 1122	A FEM 1133	A FEM 3311	A FEM 3333	A FEM 1212	A FEM 1313
	0.2 11.0342 -1.5470 0.4701 1.6971 11.5571 7.5756 3.1521
	40 0.25 6.9234 -1.0207 0.2974 1.1060 7.3017 4.7760 2.1095
	0.3	30 5.1440 -0.7928 0.2216 0.8517 5.4701 3.6464 1.7214 2.4092 1.8562
	0.35 2.2013 -0.3019 -0.0646 2.0944 0.35 4.1841 -0.6706 0.1799 0.7189 4.4972 3.1302 1.5559 8.0402 2.5032 1.5821
	0.4 2.2674 -0.3064 -0.0218 1.4539 20 0.4 3.6067 -0.5955 0.1537 0.6410 3.9260 2.8884 1.4771 6.0956 2.5737 1.4733
	0.45 2.3159 -0.3101 0.0106 0.45 3.2362 -0.5423 0.1364 0.5910 3.5668 2.7761 1.4360 1.0948 4.9690 2.6259 1.4255 0.5 2.3494 -0.3117 0.0355 0.8777 4.2632 2.6610 1.4019 10 0.5 2.9741 -0.4975 0.1243 0.5499 3.3036 2.7236 1.4132
	0.6 2.3910 -0.3139 0.0670 0 0.6 2.7249 -0.4420 0.1130 0.5078 3.0561 2.6957 1.3944 0.6557 3.5146 2.7046 1.3865
	0.7 2.4121 -0.3140 0.0856 -2 -1 0.7 2.5883 -0.3976 0.1086 0.4781 2.9065 2.6973 1.3912 0.5505 3.1416 2.7258 1.3870 0 1 2
	0.8 2.4228 -0.3134 0.0974 0.8 2.5102 -0.3629 0.1078 0.4595 2.8163 2.7094 1.3951 0.4942 2.9321 2.7358 1.3933
	0.9 2.4281 -0.3126 0.1054 0.9 2.4623 -0.3351 0.1089 0.4482 2.7586 2.7256 1.4024 0.4616 2.8035 2.7403 1.4018
	1 1	2.4312 -0.3122 0.1112 2.4314 -0.3123 0.1111 0.4417 2.7202 2.7435 1.4113 0.4417 2.7200 2.7431 1.4112
	1.5 2.4300 -0.3090 0.1265 • Superspherical pores 1.5 2.3776 -0.2394 0.1262 0.4389 2.6479 2.8318 1.4578 0.4123 2.5554 2.7386 1.4534
	2 2	2.4281 -0.3081 0.1340 2.3733 -0.1983 0.1402 0.4500 2.6375 2.9019 1.4953 0.4152 2.5211 2.7356 1.4858
	2.5 2.4267 -0.3074 0.1388 2.5 2.3816 -0.1723 0.1513 0.4623 2.6441 2.9613 1.5252 0.4234 2.5159 2.7338 1.5098
	3 3	2.4272 -0.3075 0.1425 2.3911 -0.1535 0.1602 0.4734 2.6531 3.0068 1.5483 0.4321 2.5203 2.7342 1.5285
	4 4	2.4290 -0.3075 0.1476 2.4095 -0.1283 0.1735 0.4912 2.6718 3.0745 1.5817 0.4470 2.5342 2.7361 1.5551
	5 5	2.4314 -0.3078 0.1512 2.4230 -0.1114 0.1829 0.5042 2.6860 3.1190 1.6040 0.4584 2.5478 2.7385 1.5731

Table E .

 E 11: Numerical estimation of A i jkl for the superspherical pore embedded in a transversely isotropic corrected model with different values of concavity p ∈ [0.2, 5].

This comparison is made except for the gradient G 0 i j,k on the symmetry axis due to the singularity of the reference solution derived fromPouya (2007a).

Note also that the corrected boundary conditions can also be applied to any type of anisotropy. It probably requires supplementary simulations and is not further pursued here.

In the case of isotropic elasticity of the matrix material, the Young modulus and coefficient of Poisson are respectively supposed as E = 1GPa, ν = 0.3. Note that the computation result of the contribution tensr does not depend on their values.

Appendix A. Background

Appendix A.1. Notations Barred letters A, C, D, Q refer to fourth order tensors, bold letters ε, σ, i refer to second order tensors, underlined letters x, x refer to first order tensors. Einstein's summation convention over repeated indices is used unless otherwise indicated. ⊗, : and :: respectively represent tensor product, (dot product), contracted products on two and four indices. i, I, J and K = I -J respectively represent the second-rank identity tensor, the fourth-rank symmetric identity tensor, and fourth-rank spherical and deviatoric isotropic projectors (δ i j denotes Kronecker delta symbol,

Appendix A.2. Fourth-order transversely isotropic tensors in Walpole tensor basis We use in this paper the following fourth order tensor basis for transversely isotropic TI tensors (see [START_REF] Kunin | Elastic Media with Microstructure[END_REF]; [START_REF] Walpole | Fourth-rank tensors in the thirty-two crystal classes: multiplication tables[END_REF]) and post of Sébastien Brisard on github, http://sbrisard.github.io/, entitled Decomposition of transverse isotropic, fourth-rank tensors. By denoting n the unit vector of symmetry axis of the material, let us introduce the second-order tensors

In the particular case of n = e 3 , (A.5) writes

One introduces fourth-order tensors

It may be shown that any transversely isotropic fourth-order tensor H can be decomposed as

Considering symmetry axis equal to n = e 3 , corresponding relations with usual components write

Appendix B.1. Green tensor of the infinite medium and its gradient in the general anisotropic case Appendix B.1.1. Green tensor The general expression of the Green tensor derived from a reasoning based on the Fourier transform [START_REF] Mura | Micromechanics of defects in solids[END_REF] or plane-wave expansion [START_REF] Willis | Micromechanics determination of electroelastic properties of piezoelectric materials containing voids[END_REF] writes

where δ is the scalar Dirac distribution and K(ξ) is the acoustic tensor defined from the fourth order elastic tensor C by

2)

The expression (B.1) can be further simplified by using the spherical coordinates associated to the pole x 0 in which the integration variable ξ writes with the variable change z = cos θ ξ = cos θ e + sin θ u ϕ = z e + 1z 2 u ϕ (B.3)

where e = x/ x and u ϕ is the unit vector parametrized by ϕ describing the unit circle centered on 0 and orthogonal to x.

Recalling that δ(λ z) = δ(z)/λ for all λ > 0, it finally comes that Bonnet ( 2009)

. Gradient of the Green tensor

Taking the gradient of (B.1) provides the following third order tensor grad

Observing now that δ (λ z) = δ (z)/λ 2 for all λ > 0, using the definition of the derivative in the sense of distributions [START_REF] Schwartz | Théorie des distributions[END_REF][START_REF] Gel'fand | Generalized Functions[END_REF]) and still adopting the spherical coordinates associated to the pole x 0 yield the following expression

and finally using (B.2), (B.3) and the derivative of the inverse (

It may be noticed that B.7 corresponds to relation 5.58 given in [START_REF] Mura | Micromechanics of defects in solids[END_REF], p. 33. Fourier transform based derivation and change of variable on the unit sphere used in this paper are similar to the ones detailed in [START_REF] Mura | Micromechanics of defects in solids[END_REF].

Appendix D. Complementary results concerning the strain concentration tensor in the case of γ = 1/5 and the compliance contribution tensor for the ellipsoidal pores with γ = 1, 1/2, and 1/10

We aim at display some complementary results respectively in Appendix D.1 for the strain concentration tensor A in the case of γ = 1/5 and in Appendix D.2 concerning some other shapes of the ellipsoidal pore, such as γ = 1 Appendix E. Complementary results of the strain concentration tensor A for the superspheroidal and superspherical pores

In this section, we report respectively in E.10 and E.11 the numerical computation of the strain concentration tensor of the superspheroidal and superspherical pores as described in Section 5. The corresponding evolution of their components in function of the concavity parameter p ∈ [0.2, 5] are also illustrated in Figs.

• Superspheroidal pores