Mingchao Liu

Jian Wu

Yixiang Gan
email: yixiang.gan@sydney.edu.au

Dorian Ah Hanaor

C Q Chen
email: chencq@tsinghua.edu.cn

Multiscale modeling of effective elastic properties of fluid-filled porous materials

Keywords: Porous materials, Effective elastic properties, Micromechanical model, Double-porosity, Fluid diffusion

Highlights

• A multiscale micromechanical homogenization framework is proposed for fluid-filled porous materials.

• The overall elastic properties of the fluid-filled porous materials are determined in the framework.

• The effects of pore distribution and fluid pressure on the effective elastic properties are quantified.

• By considering the fluid diffusion from macro-to micro-scale pores, evolution of the effective properties is predicted

Introduction

As a type of heterogeneous materials, porous materials are usually composed of two components, a solid skeleton and pores.

Pores may be either empty or filled with fluid (i.e., gas or liquid) [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF][START_REF] Gibson | Cellular Materials in Nature and Medicine[END_REF]. The latter type is classified as fluid-filled porous materials, which are widely available in not only natural geophysical and biological systems, e.g., rock, plant, and bone, but also synthetic structures, such as hydrogels, tofu, and fluidic origami (Kim and Guyer, 2014).

Predicting the effective elastic responses of fluid-filled porous materials is a challenging problem, which has attracted considerable attention in recent years, driven by its importance in diverse contexts such as estimating the stored amount of carbon dioxide (CO2) during the injection in underground saline aquifer [START_REF] Streit | Estimating fault stability and sustainable fluid pressures for underground storage of CO2 in porous rock[END_REF][START_REF] Leung | An overview of current status of carbon dioxide capture and storage technologies[END_REF][START_REF] Szulczewski | Theoretical analysis of how pressure buildup and CO2 migration can both constrain storage capacity in deep saline aquifers[END_REF], predicting the time-dependent response of hydraulic actuators [START_REF] Guiducci | Pressurized honeycombs as soft-actuators: a theoretical study[END_REF][START_REF] Yuk | Hydraulic hydrogel actuators and robots optically and sonically camouflaged in water[END_REF], and designing fluid-filled acoustic metamaterials [START_REF] Spadoni | Closed-cell crystalline foams: self-assembling, resonant metamaterials[END_REF][START_REF] Dorodnitsyn | Twodimensional fluid-filled closed-cell cellular solid as an acoustic metamaterial with negative index[END_REF].

The particularity of this kind of materials is the strong coupling between the fluid pressure in the pores and the elastic deformation of the solid skeleton. This coupling induces a complex problem, as the overall elastic responses of fluid-filled porous materials are governed by the fluid pressure in the pores [START_REF] Shafiro | Materials with fluidfilled pores of various shapes: Effective elastic properties and fluid pressure polarization[END_REF][START_REF] Warner | The elasticity and failure of fluid-filled cellular solids: theory and experiment[END_REF][START_REF] Ayyagari | On the nature of pressure dependence in foams[END_REF]. To describe the coupling of pore pressure, solid deformation, and macroscale stresses, in fluid saturated porous media, linear poroelasticity was established by [START_REF] Biot | General theory of three-dimensional consolidation[END_REF] through a phenomenological approach. The thermodynamic background of this theory has since been confirmed by the mixture theory [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF] and the homogenization method [START_REF] Thompson | A reformation of the equations of anisotropic poroelasticity[END_REF], as well as by the micromechanical approach [START_REF] Cheng | Material coefficients of anisotropic poroelasticity[END_REF].

More recently, Biot's constitutive relations for a linearly poroelastic material were also extended to cover the cases on finite deformation [START_REF] Brown | Effective equations for fluid-structure interaction with applications to poroelasticity[END_REF]. Furthermore, if the solid skeleton is not a 4 continuum, i.e., there are some micro-scale pores in the solid matrix with a significantly smaller size than the macro-scale pores, the fluid in macro-scale pores will diffuse into the matrix (i.e., the small pores) driven by the pressure gradient between the two pore scales, resulting in the evolution of effective properties [START_REF] Berryman | The elastic coefficients of double-porosity models for fluid transport in jointed rock[END_REF][START_REF] Rohan | Multiscale modeling of a fluid saturated medium with double porosity: Relevance to the compact bone[END_REF][START_REF] Song | Shear properties of heterogeneous fluid-filled porous media with spherical inclusions[END_REF]. The evolving effective properties of fluid-filled porous materials are of importance to many practical applications and deserve further investigation towards a systematic understanding.

There has been a series of efforts made to understand the overall elastic responses of porous materials. As a simple case, the elastic constants and their porosity dependence in dry porous materials, i.e., pores are empty or filled with air, have attracted much attention [START_REF] Mackenzie | The elastic constants of a solid containing spherical holes[END_REF][START_REF] Walsh | Effect of porosity on compressibility of glass[END_REF][START_REF] Wang | Young's modulus of porous materials[END_REF][START_REF] Day | The elastic moduli of a sheet containing circular holes[END_REF][START_REF] Hu | Effective elastic properties of 2-D solids with circular holes: numerical simulations[END_REF][START_REF] Pabst | Elasticity of porous ceramics-A critical study of modulus -porosity relations[END_REF][START_REF] Li | Effective elastic properties of randomly distributed void models for porous materials[END_REF][START_REF] Chen | Microstructural characteristics and elastic modulus of porous solids[END_REF]Chen et al., 2017).

However, the role of the spatial distribution of pores has rarely been considered. [START_REF] Goussev | Void-containing materials with tailored Poisson's ratio[END_REF] showed that the transverse Poisson's ratio of porous materials depends strongly upon the spatial arrangement of the pores. More recently, Liu et al. (2016 (a) and (c)) also identified the effect of pore distribution on the elastic responses of ordered porous materials to inner pressure.

Naturally occurring porous materials, such as rocks and plants, generally contain fluidfilled pores, and effects of pressure on their overall elastic properties cannot be ignored [START_REF] Brown | On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid[END_REF][START_REF] Christensen | The influence of pore pressure and confining pressure on dynamic elastic properties of Berea sandstone[END_REF][START_REF] Georget | Modelling of carrot tissue as a fluid-filled foam[END_REF][START_REF] Gibson | Cellular Materials in Nature and Medicine[END_REF]. More recently, synthetic porous materials have been designed by incorporating pressurized fluid in pores for use as tunable systems [START_REF] Guiducci | Pressurized honeycombs as soft-actuators: a theoretical study[END_REF][START_REF] Lv | A multiscale corotational method for geometrically nonlinear shape morphing of 2D fluid actuated cellular structures[END_REF][START_REF] Yuk | Hydraulic hydrogel actuators and robots optically and sonically camouflaged in water[END_REF]. Although intrinsically phenomenological, substantial work on the theoretical modeling of effective responses in fluid-filled porous materials can be found in the literature [START_REF] Shafiro | Materials with fluidfilled pores of various shapes: Effective elastic properties and fluid pressure polarization[END_REF][START_REF] Warner | The elasticity and failure of fluid-filled cellular solids: theory and experiment[END_REF][START_REF] Kitazono | Application of mean-field approximation to elasticplastic behavior for closed-cell metal foams[END_REF][START_REF] Vincent | Porous materials with two populations of voids under internal pressure: I. Instantaneous constitutive relations[END_REF][START_REF] Ma | Elastoplastic mechanics of porous materials with varied inner pressures[END_REF][START_REF] Guo | Analysis of vapor pressure and void volume fraction evolution in porous polymers: a micromechanics approach[END_REF][START_REF] Su | A pressure-dependent phenomenological constitutive model for transversely isotropic foams[END_REF]. Different methods for the prediction It must be noted that, in the various previous studies, the solid skeleton was assumed as a continuous phase with its interaction with pore fluid occurring only at the pore wall, assumed to be impermeable. However, many porous materials, such as rocks and bones, exhibit multiscale pore structures [START_REF] Cowin | Bone Mechanics Handbook[END_REF][START_REF] Tsakiroglou | A new approach for the characterization of the pore structure of dual porosity rocks[END_REF][START_REF] Borgomano | Dispersion and attenuation measurements of the elastic moduli of a dual-porosity limestone[END_REF]. In such materials, the structure can be considered as two or more interacting pore systems, which collectively have a strong influence on fluid transfer and effective elastic properties. The transport of fluid in such multiscale porous materials, especially the double-scale cases has been the subject of numerous studies [START_REF] Moutsopoulos | Hydraulic behavior and contaminant transport in multiple porosity media[END_REF][START_REF] Ba | Mesoscopic fluid flow simulation in double-porosity rocks[END_REF][START_REF] Choo | Hydromechanical modeling of unsaturated flow in double porosity media[END_REF].

Homogenization technologies have been

proposed to explore the elastic behaviors of multi-(or double) scale porous materials (Auriault andBoutin, 1992, 1993;[START_REF] Boutin | On models of double porosity poroelastic media[END_REF][START_REF] Rohan | Double porosity in fluid-saturated elastic media: deriving effective parameters by hierarchical homogenization of static problem[END_REF]. However, to date, few attempts have been made to incorporate the diffusion effect into 3

1 r 1 R  0 r 0 R  x y (a) (b) (c)

Effective elastic properties of fluidfilled porous materials

Fluid pressure acting on pore walls can have a significant effect on the homogenized mechanical behavior of porous materials. In this section, we examine the overall elastic responses of fluid-filled porous materials based on the single-scale micromechanical model as described in Section 2.

Elastic properties of dry porous materials

We first consider a dry porous material, in which the pores are empty or filled with air.

 o  p  o  p     (a) (b) (c) o p o  i r i R
Cite as: Liu, M., Wu, J., Gan, Y., Hanaor, D. A., & Chen, C. Q. (2019). Multiscale modeling of the effective elastic properties of fluid-filled porous materials. International Journal of Solids and Structures,162,[36][37][38][39][40][41][42][43][44] 8

In order to calculate the bulk modulus of porous materials through the single-pore model, the inner and outer boundaries have By superposing the two components, the outer boundary condition of the unit cylinder in Fig. 2(a) can be obtained as

  1 o o i p   .
In the case of plane stress, the radial displacement of the thick-wall cylindrical unit cell can be calculated according to the classical theory of elasticity [START_REF] Timoshenko | Theory of Elasticity[END_REF])

      2 2 2 2 2 2 2 1 1 1 1 i o i i o r i i i i i i i R p R r p u r v r v E R r R r r             , (1)
where

    2 2 =2 r i i r i ii u R R u R RR      . ( 2
)
9 Substituting Eq. ( 1) into Eq. ( 2), the effective bulk modulus of dry porous materials can be obtained by the definition

= o K  as         1 2 1 1 1 i i i i i i i E K vv            . (3)
It can be found that the effective bulk modulus of dry porous materials (i.e., It is noted that the geometrical arrangement of pores can have a significant effect on the overall elastic deformation. In particular, the transverse Poisson's ratio of a system having triangularly packed pores increases with porosity while that of a square packed array decreases [START_REF] Goussev | Void-containing materials with tailored Poisson's ratio[END_REF]. This porosity dependence of Poisson's ratio cannot be captured by existing theoretical models. To fill this gap, we present a semiempirical equation by adopting the geometric factor α, introduced by Liu et al. (2016 (a) and (c)) to describe the effect of pore distribution on the pore-load modulus, to relate the effective Poisson's ratio to the porosity of the porous material as

          2 3 2 1 1 1 1 1 Exp 2 2 2 i i i i i i i i v v v                                        . ( 4
)
For an isotropic elastic material, there are only two independent elastic parameters. For a planar problem, the effective Young's modulus can be obtained through the relation 3) and ( 4) as

of 1 = 2 i i i v E K  combining with Eqs. (
      3 3 2 1 1 Exp 1 2 2 2 i i i i i i i EE                                  . ( 5
)
By means of Eqs.

(3)-( 5), the elastic properties of a porous material can be related to the porosity and the arrangement of pores.

It should be mentioned that it should also be noted that the stiffness of a 10 porous material having a square lattice can be considered as isotropic only under hypothetical conditions in which porosity is sufficiently low and the so-called macropores are distributed quite sparsely. However, from a more accurate perspective, the stiffness of materials with a square lattice pore arrangement is orthotropic, and two stiff and two compliant directions can be found [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF]. Due to the symmetrical pore distribution, the relation of

1 = 2 i i i v E K  is also

Effect of fluid pressure

When pores are filled with pressurized fluid, the deformation of a solid skeleton can be significantly affected by the fluid pressure.

Accordingly, the overall elastic properties of these materials depend upon fluid pressure [START_REF] Warner | The elasticity and failure of fluid-filled cellular solids: theory and experiment[END_REF][START_REF] Gibson | Cellular Materials in Nature and Medicine[END_REF].

Here

      10 1 x f i x L P M L  , ( 7
)       20 +1 x f i e i x L P M E L   . ( 8
)
Here,

  0 x L is the x-direction
         1 1 1 1 1 ii i i i i i i E M vv             . ( 9
)
By substituting Eqs. ( 7) and (8) into Eq. ( 6), the effective Young's modulus can be with an explicit form as

1 i i f i E E P M  , ( 10
)
where The effective bulk modulus can be obtained according to the similar analysis through the micromechanical model, and the corresponding expression is

1 i i f i K K P M  , ( 11
)
where i K is the effective bulk modulus of the dry porous material given by Eq. ( 3).

The linear dependence of the effective bulk modulus on the pore pressure is further evident. Considering the relation of

1 = 2 i i i v E K 
, by combining Eqs. ( 10) and ( 11), the effective Poisson's ratio of the fluid-filled porous material can be found independent of the pore pressure, i.e., = ii vv . and can be described by Darcy's law [START_REF] Whitaker | Flow in porous media I: A theoretical derivation of Darcy's law[END_REF]. The fluid front r , i.e., the interface between fluid-filled micro-scale pores and empty ones, see the blue line in Fig. 6(a), moves from the inner boundary 1 r to the outer boundary 1 R of the cylinder unit cell during diffusion. The variation of the fluid front against diffusion time can be obtained as [START_REF] Conrath | Radial capillary transport from an infinite reservoir[END_REF]Liu et al., 2016 (b)) 10) and ( 11) with 0 i  . These two layers form the new skeleton of the macro-scale porous structure which can be regarded as a composite structure with two phases [START_REF] Wang | Tailoring the moduli of composites using hollow reinforcement[END_REF].

2 1 1 0 11 ln = 22 r r t r r t              , ( 12
Under the inner and outer loads, as shown in 16 [START_REF] Christensen | Mechanics of Composite Materials[END_REF]. The homogeneous Young's modulus of this composite skeleton can be obtained as

  1 0 0 1 E E E       , ( 13
)
where 0 E and 0 E are Young's moduli of the dry outer layer and fluid-filled inner layer, respectively, and  is the volume fraction of the inner fluid-filled layer with

    2 2 2 2 1 1 1 = r r R r   . ( 14
)
Here, r is the fluid front position, and 1 r  and 1  is depending on the arrangement of the macropores. Similar to Eq. ( 13), the homogeneous bulk modulus can be expressed as

  1 0 0 1 K K K       , ( 15
)
where 0 K and 0 K are bulk moduli of the dry and fluid-filled layers, respectively. As the effective Poisson's ratio of the fluid-filled layer is independent of fluid pressure, i.e., 00 = vv , the homogeneous Poisson's ratio can also be obtained as 10 = vv .

The effective elastic properties of the macroscale dry porous structure containing partial filled micro-scale pores (see Section 3.1), i.e., 1

K , 1 v , 1
E and 1 M , can be obtained from Eqs. ( 3)-( 5), (9) by setting subscript 1 i  , with elastic constants of the homogeneous solid matrix obtained from Eqs. ( 13) and ( 15). It should be mentioned that the pore-load modulus of macro-scale dry porous structure 1 M can also be calculated through the theoretical model given by Liu et al. (2016 (c)). To unify, here we use the equivalent model of Eq. ( 9) combined with the homogeneous parameters.

Similarly, as shown in Section 3.2, the effective properties of the macro-scale fluidfilled porous structure containing partially filled micro-scale pores (i.e., 1

K , 1 v , 1 E )
can be calculated through Eqs. ( 10) and (11) with a subscript 1 i  . (3), ( 5) and ( 9) with 1 i  ) are also included as solid lines for comparison. It can be found that the theoretical predictions can consist with the FEM simulations well for both cases. With increasing volume fraction of the fluid-filled layer, all moduli increase. This is because the modulus of the fluidfilled layer is larger than the dry layer, and Cite as: Liu, M., Wu, J., Gan, Y., Hanaor, D. A., & Chen, C. Q. (2019). Multiscale modeling of the effective elastic properties of fluid-filled porous materials. International Journal of Solids and Structures,162,[36][37][38][39][40][41][42][43][44] 18 stiffens the composite structure.

Evolution of the effective properties during diffusion

When fluid is injected into the doubleporosity structure, there are two typical stages. Firstly, the macro-scale pores are filled rapidly, incurring the stiffening of the structure. Secondly, pressurized fluid in the macro-scale pore diffuses into the microscale pores imbedded in the matrix, driven by the pressure gradient, resulting in complex changes of the overall elastic behavior.

Relative to the second stage, the first one is much faster. Here we focus on the evolution of the effective properties during the fluid diffusion in the second stage. Two typical injection conditions, i.e., constant pressure and constant injection rate, are considered. These conditions are common in practical engineering problems, such as CO2 storage and hydraulic fracture [START_REF] Coninck | Carbon dioxide capture and storage: issues and prospects[END_REF][START_REF] Liu | A micromechanical analysis of the fracture properties of saturated porous media[END_REF][START_REF] Detournay | Mechanics of hydraulic fractures[END_REF]. Following the filling of the macro-scale pores, the pressure can be held constant by controlling the fluid injection rate or alternatively, the injection rate can be held constant resulting in a change of pressure as fluid is redistributed during the diffusion process.

For the case of constant injection rate, based on Boyer's Law for compressible fluids, a simple assumption can be employed to calculate the variable fluid pressure f P as  

00 ff P V P V V t      , (16)
where, following the initial filling of macropores, 0 3)-( 5), (9) (i=0)

Set:

Calculate: from Eqs. ( 10), (11) (i=0)

Calculate: from Eq. ( 14)

Calculate: from Eqs. ( 13), ( 15)

Calculate: from Eqs. ( 3), ( 5), (9) (i=1)

Calculate: from Eqs. ( 10), (11) (i=1)

Calculate: from Eq. ( 12)

Set:

Calculate: from Eq. ( 17

  of overall mechanical properties of fluidfilled porous materials have been proposed by considering pressurized fluid as an compressible inclusion. To evaluate 5 poroelasticity, by introducing the Biot modulus as a parameter, the effect of fluid pressure on the effective response can be considered. However, as a phenomenological theory, the effect of the geometric distribution of pores is not captured by the Biot theory.

  homogenization techniques, and to further address the evolution of effective properties induced by fluid diffusion in the skeleton of porous materials. In this paper, a two-dimensional (2D) micromechanical model, which consists of pores at two scales, is employed to develop a multiscale framework to analyze the overall elastic responses of the fluid-filled porous materials. This model allows for the coupling of fluid pressure and the deformation of the solid matrix at each pore scale and for fluid diffusion from macro-scale to micro-scale pores imbedded in a non-continuum skeleton driven by a pressure gradient. Based on this micromechanical model, a theoretical model is proposed to predict the effective elastic properties of porous materials with/without pressurized pore fluid and is validated by finite element method (FEM) simulations. With this model, time-dependent elastic properties are obtained for constant and variable fluid pressure cases. elastic solid. At the macroscale, see Fig. 1(a), the cylindrical pores are arranged in a well ordered 2D triangular lattice, and all pores are filled with pressurized fluid. A unit cell comprising a thick-wall cylinder from the macro-scale porous structure is shown in Fig. 1(b), and micro-scale cylindrical pores are distributed in the cylinder wall. The micro-scale pore structure within the solid skeleton is shown in Fig. 1(c), with a similar triangular distribution. It should be noted that the considered macro-scale pores are much larger than the micro-scale pores and much smaller than the dimensions of the overall structure. Accordingly, the porous material can be treated as a continuum solid in representing its overall mechanical behavior. Moreover, we assume the presence of fine diffusion pathways within the solid matrix, facilitating the transport of fluid from macroto micro-scale pores. Additionally, fluid flow between pores at the same scale is assumed to be restricted. For simplicity, the pores within each level are assumed with the same size.

Fig. 1 .

 1 Fig. 1. Schematic representation of a porous material with double porosity: (a) Marco-scale pores filled with pressurized fluid; (b) a unit cell containing pores at two scales; (c) the solid skeleton containing micro-scale pores. Only the triangular distribution of pores is shown as an example.

  Since the ambient pressure and the compressibility of air are much smaller than the elastic stiffness of the solid matrix[H1], their effects on the deformation of the solid skeleton are negligible. The effective bulk modulus can be obtained by analyzing the deformation of porous samples under equibiaxial external load, as shown in Fig. 2(a). It is not easy to directly solve the stress and deformation fields of this structure, on account of the interaction of microstructures (i.e., the ordered pores arranged in the matrix). However, Gor et al. (2015) pointed out that the overall deformation of a plate with many pores can be approximately represented by the deformation of a unit cell comprising a homogeneous cylinder. Our previous works (i.e., Liu et al. 2016 (a) and (c)) further showed that each unit cell is affected by its neighboring cylinders, implying that the outer boundary conditions of the cylinder are related to the geometrical arrangement of the pores. This single-pore model will be extended here to calculate the bulk modulus of the dry porous material.

Fig. 2 .

 2 Fig.2. Schematics of the superposition method used to determine the bulk modulus of porous materials: (a) a porous sample subjected to external equi-biaxial tension, which can be superposed by (b) a porous sample subjected to external equi-biaxial tension and the uniform pore pressure, and (c) a porous sample subjected to uniform pore pressure. The blue dotted lines represent the outer boundaries of the cylindrical unit cell.

  o

.

  deformation as shown in Fig. 2(b), i.e., the same uniform loads are applied on the outer boundary of the sample and the pore surfaces; II. Pore-load deformation, as shown in Fig. 2(c), i.e., the internal pressure is only applied on the pore surface, and the outer boundary of the sample is free. For the uniform deformation case, the outer boundary of the unit cylinder can For the pore-load deformation case, Liu et al. (2016 (a) and (c)) gave the outer boundary condition as io p    , where  is a geometrical factor related to the interaction of neighboring pores and depending on the arrangement of pores, i.e.,  = 1/3 for porous materials with a triangular lattice. For materials with a square lattice pore arrangement, considering the symmetry of both the structure and the loading conditions, the applicability of the single-pore model is also valid, with  = 11 refer to the micro-and macro- scale pores, respectively.

  the homogenized elastic constants of the non-continuum skeleton. Furthermore,

  Fig. 3. Comparison of FEM simulations (symbols) and theoretical predictions (lines) of the effective elastic constants of porous materials with triangular (TR) and square (SQ) lattice as a function of porosity: (a) the normalized effective bulk modulus; (b) the normalized effective Young's modulus; (c) the effective Poisson's ratio.

  Cite as:Liu, M., Wu, J., Gan, Y., Hanaor, D. A., & Chen, C. Q. (2019). Multiscale modeling of the effective elastic properties of fluid-filled porous materials. International Journal ofSolids and Structures, 162,[36][37][38][39][40][41][42][43][44] The FEM simulated porosity dependent effective elastic constants of the porous material are shown in Figs. 3(a)-(c) for the bulk modulus, Young's modulus and Poisson's ratio, respectively. For each case, pore distributions of triangular and square lattices are considered. The corresponding theoretical predictions given by Eqs. (3)-(5) with different geometrical arrangement factors are included as solid lines for the purpose of comparison. For all three cases, one can find that the proposed theoretical model agrees well with the FEM results with different pore distributions. In particular, the opposing relationships between the effective Poisson's ratio and porosity can be qualitatively predicted by employing the geometric factor. The theoretical model of Eqs. (3)-(5) will be used to the further analysis of the effects of fluid pressure and diffusion on the overall elastic properties of fluid-filled porous materials.

Fig. 4 .L

 4 Fig. 4. Schematic representation of the deformation of a fluid-filled porous material subjected to inner fluid pressure and external load: (a) Porous structure with empty pores (pressure free); (b) The reference configuration of porous materials including pore pressure; (c) The deformed configuration of pressurized porous materials subjecting external load.   j x L and   j yL ( j = 0, 1 and 2) refer to the dimensions of the porous material under different states at x-and ydirections, respectively.

  dimension of the porous materials without fluid, f P is the fluid pressure in the porous materials, and the pore-load modulus, i M , is given by Liu et al. (2016 (a) and (c)) as

  is seen that the effective Young's modulus of the fluid-filled porous material is linearly dependent on the pore 13 pressure.

Fig. 5 .

 5 Fig. 5. Comparison of FEM simulations (symbols) and theoretical predictions (lines) of the effective elastic moduli of fluid-filled porous materials as a function of pore pressure: (a) the normalized effective Young's modulus; (b) the normalized effective bulk modulus.

  scale, in which  is the viscosity of the fluid, 0  the porosity of the solid skeleton with microscale pores, k the permeability of the porous skeleton, and f P the pressure of the macro-scale pores. The fluid front position can be predicted quantitatively by Eq. (12). It should be noted that here we consider the diffusion of fluid in the skeleton within micro-scale pores. The effect of the microstructure (e.g., the size and spatial distribution of the micro-pores) on the diffusion process can be reflected by homogenized macroscopic parameters, i.e., porosity and permeability. For example, a quantitative relation can be established to correlate the permeability and the distribution and size of micro-pores (Sobera[START_REF] Sobera | Multiscale modeling of the effective elastic properties of fluid-filled porous materials[END_REF].

Fig. 6 .

 6 Fig. 6. Schematics of the pressure-driven fluid diffusion from macro-scale pores to micro-scale pores and its induced heterogeneity: (a) A cylindrical unit cell with pressurized fluid filled macro-scale pore; (b) The interface between fluid-filled micro-scale pores and empty ones; (c) The illustration of parallel model used to homogenize the porous material with partially pressurized micro-scale pores.

Fig. 6

 6 Fig. 6(c), the parallel model can be employed to homogenize this two-phase composite



  , and cannot be completely arbitrary. In addition, the relation between 1

Fig. 7 .E

 7 Fig. 7. Comparison of FEM simulations (symbols) and theoretical predictions (lines) of the effective elastic moduli of heterogeneous porous materials as a function of volume fraction of the fluid-filled porous layer: (a) φ1 = 0.40; (b) φ1 = 0.10.

  Fig. 8. Flow chart of the iterative algorithm used for determining the variation of effective properties of fluid-filled porous materials.

Fig. 9 .t

 9 Fig. 9. Variation of normalized effective elastic moduli as functions of the normalized diffusion time under injection conditions of: (a) constant pressure; (b) constant injection rate and (c) constant volume.

Acknowledgments

This work is supported by the Science Challenge Project of China (Grant No. TZ2018007), and the National Natural Science Foundation of China (No. 11472149 and11732007). M.L. acknowledges the support from the Endeavour Research Fellowship founded by Australian Government. Y.G. acknowledges the support from The University of Sydney SOAR Fellowship and the Australian Research Council Discovery Project DP170102886.


Keywords: Sitona, Pisum sativum L., crop damage, on-farm approach

The study measures intensity of attacks by Sitona lineatus L. on pea crops in farmers' fields, and the resulting damage. Insect attacks were measured on the aerial parts of the plant at the vegetative stages, and on roots at the start of flowering of the crop. Attack intensity varied between and within plots. Early sown crops were at greater risk of severe damage. There was a good correlation between leaf damage at the beginning of the crop cycle and damage to nodules at the start of flowering in unsprayed crops. Aerial spraying did not completely prevent nodule damage in crops severely injuried by insects at the early vegetative stage. The effects of the attacks on nitrogen nutrition of the pea crop were also assessed; these results are discussed with regard to the consequences for crop management.

INTRODUCTION

The area devoted to spring peas (Pisum sativum L.) increased considerably in France during the 80's, going from less than 100,000 hectares in 1980 to almost 700 000 in 1992. This increase has been accompanied by increases in diseases and pests, especially the pea weevil Sitona lineatus L.. The effects of Sitona lineatus larvae on pea yields is well documented (see for example Bardner et al., 1979). It is considered to be one of the three most dangerous insect species for the crop in France [START_REF] Bournoville | Les insectes ravageurs des protéagineux[END_REF]. The work of [START_REF] Cantot | Quantification des populations de Sitona lineatus L. et de leurs attaques sur pois protéagineux (Pisum sativum L.)[END_REF] led to a method for estimating the infestation sizes at both the adult and larval stages. [START_REF] Cantot | Quantification des populations de Sitona lineatus L. et de leurs attaques sur pois protéagineux (Pisum sativum L.)[END_REF] also observed a relationship between the infestation size of first generation adults and the infestation size of the following generation larvae under semi-field conditions (artificial infestation under cages in a field). He identified effects of the larvae on certain productivity factors, the number of seeds and the yield per plant, under semi-field or controled conditions (hydroponic cultivation of the peas in an air-conditioned room) [START_REF] Cantot | Action larvaire de Sitona lineatus L. sur quelques facteurs de production du pois protéagineux (Pisum sativum L.)[END_REF]. The threshold for these effects is 10 eggs per plant for isolated plants. Experiments comparing the productivity of pea crops with and without insecticide treatments showed that the yield of treated plots could be enhanced, depending on treatment type and infestation size (Bardner & Fletcher, 1979;Bardner et al., 1979;[START_REF] King | Experiments for the control of pea and bean weevill (Sitona lineatus) in peas, using granular and liquid insecticides[END_REF][START_REF] Vulsteke | Lutte contre le sitone du pois[END_REF]. However, the true scale of weevil damage to peas within France is not really known, because the attacks can vary greatly between fields within an area. This fundamental information is needed in order to develop strategies for controlling Sitona. We therefore undertook this study to obtain information about on this variability and on the resulting damage.

MATERIALS AND METHODS

Farm fields in the Seine-et-Marne district of France were studied for three years (34 in 1988, 27 in 1989, and 17 in 1990). All the fields were within an area of about 1000 km², and were chosen to represent the variabilty of the cropping systems used in the area. The soils were loamy, and the main features of the climate are shown in table 1. We studied only a small plot of about 250 m² in each field in 1988 and 1989. The variability of the attacks within each field (sizes from less than 5 to over 20 hectares) was also assessed in 1990. The pea cultivars were of the afila type : Solara in 1988 and[START_REF] Cantot | Action larvaire de Sitona lineatus L. sur quelques facteurs de production du pois protéagineux (Pisum sativum L.)[END_REF]and Solara and Ascona in 1990. The sowing period of the pea crops was divided into two parts by rainfall in each year of the study. Some fields were sown early (first sowing period), others were sown about two weeks later (second sowing period; dates are given in table 2).

Measurements

The severity of infestation (Sitona first generation adults) was determined when the peas were at the "3-4 leaf" stage. The method used was derived from [START_REF] Cantot | Quantification des populations de Sitona lineatus L. et de leurs attaques sur pois protéagineux (Pisum sativum L.)[END_REF] and is based on the number of notches on the first pair of stipules of each plant. Plants were scored on a scale of 0 (no infestation) to 5 (severe infestation) (see figure 1); a score of 5 was used for foliar damage rather more severe than that observed by [START_REF] Cantot | Quantification des populations de Sitona lineatus L. et de leurs attaques sur pois protéagineux (Pisum sativum L.)[END_REF]. [START_REF] Cantot | Quantification des populations de Sitona lineatus L. et de leurs attaques sur pois protéagineux (Pisum sativum L.)[END_REF] showed that a rapid assessment by eye corelated well with the number of notches for grades 0 to 4. We therefore scored 40 to 100 plants per plot by eye. The score for a plot was the median value of the marks per plant (0 to V). The median was used rather when the mean because it is not possible to calculate the mean of qualitative scores. The values for plots were recorded directly by fast visual assessment in 1990, without using the marks on an individual plant basis. We evaluated the variability of the attacks in this third year, by scoring a number of plots equal to the number of hectares of the field rather than a plot about 250 m²/field.

In 1990 we counted notches on leaves 5-8 at the 7-8 leaf stage. These notches in plots which had scored 0 at the 3-4 leaf stage indicated that new weevils were present on this part of the field. If notches were observed on plants which had scored more than 0 at the 3-4 leaf stage, the adult pest had not died or left between the two testing dates.

Larvae are generally present on the nodules from the first nodules formation until nodules senesce during pod filling in most areas of the world where the weevil and peas coexist. The start of flowering is the "best" point during crop development for sampling the larvae of Sitona lineatus L., because it is generally considered to be the time when most, if not all, larvae that will complete development are present. We assessed the damages to nodules due to larvae at this stage by counting the number of plants in samples of 40 to 80 plants per plot in which the nodules were destroyed by the larvae on the uppermost 10 cm of the primary root (dig to remove the plants and wash with water). According to [START_REF] Tricot | -Establishment and distribution of pea primary root nodules (Pisum sativum L.) as affected by light shading[END_REF], this corresponds to the zone where most of the nodules of the primary root are located at this stage of development. Every plot was measured in this way in 1989, and 11 plots in 1990. The mean dry weight of nodules per plant on the same part of the root was determined. Roots were washed in water, the nodules were removed from the roots, dried for 24 hours at 80°C and weighed.

We also measured the nitrogen content of the aerial parts of the pea crops at the same time. Dry matter was measured in 6 samples of about 0.5 m² per plot, and the nitrogen content was determined by the Kjeldhal method.

RESULTS

Abundances of adults at the 3-4 leaf stage

Table 3 shows the results of counts at the stage 3-4 leaf for each of the three years.

The attacks were minor in all fields in 1988, and were more serious in many fields in 1989 and 1990. Attack variability was assessed in 1990. The patterns within fields were very variable: in some cases the infestation was evenly spread and in others it varied across the field. We used numeric frequency distributions, and set up two classes. The fields which had at least 3 different scores among the tested plots were classified as heterogeneous (H), the others were classified as homogeneous (h). Table 4 shows the distributions in 1990 (h or H), for the two sowing periods (first or second sowing period) and the size of the field (3 classes). Field size had no effect on pest distribution type (χ² test, P>0.05), but the sowing period had a highly significant effect (χ² test, P<0.01). Distribution was heterogeneous in 8 of the 9 cases when sowing was early, whereas the 8 fields in which the sowing was late had a homogeneous distribution of mild attacks (never higher than I). The pattern of infestation was not assessed in 1988 and 1989 but the early sown fields had more high marks than the late sown peas in 1989 (table 5). This is consistent with the observations for 1990. The number of early-sown fields was too small in 1988 for conclusions to be drawn.

Changes in the attacks between the 3-4 leaf and 7-8 leaf stages in 1990

Twenty-three of the 81 early-sown were sprayed with insecticide (lindane, 400g/ha, or cypermethrin, 25g/ha, or lambda-cyhalotrine, 8g/ha) by the farmers at crop emergence or around the 3-4 leaf stage. The damage at the 7-8 leaf stage was scored, and classified according to the damage at the 3-4 leaf stage and insecticide treatment (table 6). All unsprayed plots scoring II or more at the 3-4 leaf stage, and many scoring I, suffered further attack. Many of the sprayed crops also suffered further attack, especially on the plots scoring IV or V at the 3-4 leaf stage. No crop sown in the second period was sprayed. About half of the 78 plot were damaged on the upper leaves, but the damage was always scored as less than 5 notches per leaf.

Damage at the start of flowering

We compared the seriousness of attacks by Sitona adults at the 3-4 leaf stage to larval damage of the following generation at the flowering onset stage in 1989 and 1990 (table 7).

Only one of the fields was sprayed with insecticide before the start of flowering. When the score was low at the 3-4 leaf stage, the fraction of plants with destroyed nodules was generally low. When the score was high, the fraction of damaged nodules was higher. Figure 2 shows the relationship between the fraction of plants in which nodules were destroyed and the dry matter of nodules per m². The greatest fractions of plants affected in 1989 were in plots of plants with low nodule dry matter. The nodule dry matter was very variable in the plots in which few plants were damaged by Sitona. This is presumably due to variability in nodulation between plots for reasons other than Sitona infestation (see for example the review by [START_REF] Sprent | Environmental effects on nitrogen fixation[END_REF]. For this reason, the correlation coefficient of a linear regression between fraction of plants with nodules destroyed and nodule dry weight is only r = 0.55. Some of the plots observed in 1990 were in the same field, and are identified by the same letter. The fraction of plants within a field having destroyed nodules is a good measure crop damage. As in 1989, differences between fields lead to a low correlation coefficient (r = 0.77).

The relationship between the nodule dry matter per m² and the amount of N accumulated in the aerial parts per m² at the start of flowering is shown in figure 3. The general distribution of the points for each year shows that low nodule weights correspond to low amounts of accumulated nitrogen. Points corresponding for plots with high nodule damage are identified. They show the great effect of Sitona feeding on pea crop N nutrition. The correlation coefficients were respectively 0.63 in 1989 and 0.64 in 1990.

DISCUSSION

The intensity of attacks by sitona on legumes is known to vary from year to year (see for example [START_REF] Vulsteke | Lutte contre le sitone du pois[END_REF]. The increase between 1988 and 1990 in the area we studied is consistent with previous observations in other areas of northern France. It has sometimes been suggested that this increase could be due to the increase in pea crop areas.

We disagree with this assumption for three reasons. First, according to [START_REF] Hamon | -Flight periodicity and infestation size of Sitona lineatus[END_REF], the fraction of legumes in the total cultivated area has no effect on the seriousness of Sitona attacks. Second, the area devoted to spring pea crops increased between 1988 and 1990 in France, but not in the area of our study. And third, the area under peas decreased slightly in 1991, yet the intensity of Sitona attacks continued to increase. There were very few attacks in 1992 and 1993, although the pea crop area again increased. Thus, frther studies are necessary to dissociate the effect of the pea crop area from the effect of climatic pattern on attack intensity. In one such study, Schotzko & O'Keefe (1988) showed that the climate during the winter influenced the ability of adults to move in spring.

Between-fields and within-field differences in attack intensity were described by [START_REF] Hamon | -Flight periodicity and infestation size of Sitona lineatus[END_REF]. Like these authors, we do not have complete information on the cropping histories of the fields or on the geographic position of legumes in the neighbouring fields during the previous year. Nevertheless, the data available suggest that these factors are not the only ones responsible for the observed differences between and within fields (not shown). There seems to be less risk of a severe attack on late-sown fields than on early-sown fields. This is consistent with the flight periodicity of Sitona observed by [START_REF] Hamon | -Flight periodicity and infestation size of Sitona lineatus[END_REF]. They describe little trivial flight; once individuals have located a crop, they move through it by walking. Such a pattern of infestation could also be responsible for the differences within a field, making it difficult for a farmer to get a rapid, accurate assesment of the attack intensity within a field. Our results on the effect of the sowing period are consistent with those of [START_REF] Hamon | -Flight periodicity and infestation size of Sitona lineatus[END_REF] for beans. They compared two sowing dates over several years in Rothamsted, and found that the latest sowing dates gave the highest yields. [START_REF] Vulsteke | Lutte contre le sitone du pois[END_REF] and [START_REF] Decoin | Pois de printemps: les mesures de protection à prendre[END_REF] also showed differences in attacks on early-and late-sowing peas in experimental situations. These differences are also seenin farmers fields sown at different periods in the same area. It concerns both only the mean attack intensity and its distribution, regardless of the field size. This effect of the sowing period was attributed by [START_REF] Proctor | An experiment to determine the effects of date of sowing on the yield and quality of harvesting peas 1953-56[END_REF] to differences in the sowing bed structure: late sowings result in the finest seedbed soil structures, which are worse conditions for Sitona. We have shown elsewhere (Doré, 1992-a) that this is not valid for our model area. About the half the plots sown during the second period in 1989 had rough structured seedbeds, but did not suffer from heavy attacks.

We have confirmed relationship between adult attacks at the 3-4 leaf stage and damage to nodules at the start of flowering seen by [START_REF] Cantot | Quantification des populations de Sitona lineatus L. et de leurs attaques sur pois protéagineux (Pisum sativum L.)[END_REF] for farmers' fields, at least when no insecticide is used. Three plots had little nodule damage though the mark at the 3-4 leaf was high (score IV) (table 6). We saw no evidence that this was due to low adult fertility, starving, feeding on crops other than pea, or variable winter quiescence (Schotzko & O'Keefe, 1988).

Where adults were not killed by treatment, adult movement by walking seems to be plausible explanation in few cases. Treatment did not kill all the adults under farm conditions when the infestation was heavy, confirming previous studies by [START_REF] Mcewen | The effects of aldicarb, triazophos and benomyl plus zineb on the incidence of pests and pathogens and on the yields and nitrogen uptakes of leafless peas (Pisum sativum L.)[END_REF], [START_REF] King | Experiments for the control of pea and bean weevill (Sitona lineatus) in peas, using granular and liquid insecticides[END_REF], and [START_REF] Bardner | Chemical control of the pea and bean weevil, Sitona lineatus L., and subsequent effects on the yield of field beans Vicia faba L[END_REF].

It is difficult to compare the levels of attacks we detected to those observed by others, because the pea varieties (afila or not) and the methods used to estimate damage vary from one study to another. The numbers of notches per plant have been reported: [START_REF] Vulsteke | Lutte contre le sitone du pois[END_REF] give figures from 1.35 to 35 notches per plant, but the number of leaves involved is not given; [START_REF] Mcewen | The effects of aldicarb, triazophos and benomyl plus zineb on the incidence of pests and pathogens and on the yields and nitrogen uptakes of leafless peas (Pisum sativum L.)[END_REF] give 17.2 to 71.5 notches on the stipules of afila plants, but there is no information on the stage of the plants when scored; [START_REF] Bardner | Chemical control of the pea and bean weevil, Sitona lineatus L., and subsequent effects on the yield of field beans Vicia faba L[END_REF] give 4.18 to 8.29 notches on the last leaf developed on non-afila types; the highest values given by [START_REF] Cantot | Quantification des populations de Sitona lineatus L. et de leurs attaques sur pois protéagineux (Pisum sativum L.)[END_REF] were 20 notches per plant on the first leaf and 12 notches per plant on the stipules of the first leaf, on non-afila types. We scored plants on some plots as 5; it was impossible to count the notches on the stipules of the first leaf at the 3-4 leaf stage on these plants because they were too numerous. These plots, therefore, suffered heavier attacks than those described by [START_REF] Mcewen | The effects of aldicarb, triazophos and benomyl plus zineb on the incidence of pests and pathogens and on the yields and nitrogen uptakes of leafless peas (Pisum sativum L.)[END_REF], which had been sprayed with insecticide against Sitona and had a significantly higher yield than non-treated plot. We cannot directly compare the yields obtained between fields because they depend on factors other than Sitona attacks. Nevertheless, heavy attacks decreased nodule weight in farmers' fields, with some effect on plant nitrogen accumulation. We have shown (Doré, 1992-a) that poor nitrogen nutrition due to heavy Sitona attacks or other factors is directly responsible for low numbers of peas, and for poor pea crop yields. We also showed that the amount of nitrogen accumulated in the pea crop greatly influenced the N supply in the soil for a subsequent wheat crop (Doré, 1992-b).

CONCLUSION

This survey provides a fairly detailed assessment of the effect of Sitona on the pea crop in the French Paris Basin. Comparisons between our data and the abundances noted by other authors indicate that the attacks can be sufficiently heavy and cause enough damage to justify action against these insects. However, designing strategies to kill Sitona will not be easy. The first problem will be to determine the real intensity of an attack on a field. Attacks in farmers' fields are very heterogeneous. The second problem will be to kill, if necessary, the adults. Aerial treatment will not always be satisfactory if the attack is serious. Seed treatment should help to solve this problem [START_REF] Blondel | Pois et haricots de conserve: 14 nouveaux produits homologués[END_REF][START_REF] Decoin | Pois de printemps: les mesures de protection à prendre[END_REF]. Finally, the strategy for fighting Sitona should include management of the whole crop, as shown by the effect of sowing date on the risk of heavy attack. [START_REF] Hamon | -Flight periodicity and infestation size of Sitona lineatus[END_REF] discussed sowing crops later to avoid attack, and concluded that it had little chance of being accepted by the farmers, because late sowing exposes the crop to attack by other pests, and delays the sowing date of the succeeding wheat crop. This latter argument is not valid in our study, as the pea crop is harvested in July or August, and wheat is sown in October. In any case, we have shown elsewhere (Doré, 1992a), that the sowing date should also be considered to reduce the risk of packing the soil structure in th same area. ACKNOWLEDGMENTS: We thank INRA and UNIP for financial support.

Table 1 :Figure 1 :

 11 Figure 1: Damage assesment on leaves by counting the notches on one plant (n = number of notches on the stipules of the first leaf)

Figure 2 :

 2 Figure 2: Relationship between the fraction of plants in which nodules have been destroyed (X-axis) and the nodule dry matter per m²(g.m -2 , Y-axis)-a-27 plots in 1989 and b -11 plots in 1990 (a same letter corresponds to different parts of the same field)

Figure 3 :

 3 Figure 3: Relationship between the nodule dry matter per m² (g.m -2 ,X-axis) and the amount of N accumulated in aerial dry matter (g.m -2 , Y-axis) -a -27 plots in 1989 and b -11 plots in 1990