Romain Desbats

François Bermond

Sabine Compigne

Stéphane Nicolle

Philippe Vezin

Philippe Vezin Lumped

Lumped element model of THORNT and PMHS abdomen under seatbelt and impactor loading

Keywords: abdomen response, biofidelity, lumped element model, THOR dummy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Abdominal injuries represent only a small proportion of all injuries caused by car crashes compared to those of the head, chest or limbs, but their proportion increases considerably with regard to serious to severe injuries (AIS 3 or more): 8% of AIS3+ injuries, 16.5% of AIS4+ injuries and 20.5% of AIS5+ injuries [START_REF] Elhagediab | Patterns of Abdominal Injury in Frontal Automotive Crashes[END_REF]. Besides, it is also shown that rear occupants have a 1.7 times higher risk of suffering an injury to the abdomen than front occupants [START_REF] Martin | Specificities of Rear Occupant Protection: Analysis of French Accident Data[END_REF]. Therefore, to provide the best protection to all car occupants, assessing the risk of abdominal injury using a biofidelic crash test dummy is important, especially under seatbelt loading known to be the most common source of AIS3+ injuries to the digestive system [START_REF] Elhagediab | Patterns of Abdominal Injury in Frontal Automotive Crashes[END_REF]- [START_REF] Frampton | An In-depth Study of Abdominal Injuries Sustained by Car occupants in Frontal Crashes[END_REF].

Abdomen biofidelity of crash test dummies under seatbelt loading is not included in biofidelity targets. The biofidelity of the THOR dummy proved to be limited under seatbelt and impactor loadings, as shown in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] and [START_REF] Compigne | Proposal for a Modified THOR Lower Abdomen Including Abdominal Pressure Twin Sensors[END_REF] where force-penetration responses of the THOR-NT (reprocessed data from [START_REF] Hanen | Contribution to the Improvement of Crash Test Dummies in Order to Decrease Abdominal Injuries in Road Accidents[END_REF]) were compared to PMHS corridors. Two impactor conditions and two seatbelt conditions were considered in [START_REF] Compigne | Proposal for a Modified THOR Lower Abdomen Including Abdominal Pressure Twin Sensors[END_REF]. For the impactor cases, the dummy response is above the PMHS corridors after 80 mm whereas for the seatbelt cases, the dummy response is below the PMHS corridors. Finite element models, as well as lumped element models, are tools to help the analysis of human surrogates abdomen response such as PMHS and dummy. The first lumped element model developed to simulate human behaviour under impact is the thorax model from [START_REF] Lobdell | Impact Response of the Human Thorax[END_REF]. Reference [START_REF] Kent | Thoracic Response to Dynamic, Non-Impact Loading from a Hub, Distributed Belt, Diagonal Belt, and Double Diagonal Belts[END_REF] used this model to simulate an impactor loading case, as well as a seatbelt loading case to the thorax. A simplified version of the same model was used in [START_REF] Parent | Thoracic Biofidelity Assessment of the THOR MOD KIT ATD[END_REF], where the mechanical parameters were optimised to match a given thoracic response. A lumped element model was used to describe the abdomen response under seatbelt tests in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF]. The model consisted of a spring in parallel with a damper and a mass at the front for the dummy case only. The mass value was tuned to fit THOR dummy abdomen behaviour. The PMHS version of this model was enhanced in [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF]: a second stage was added to the model, linked with a mass to the first one. This allowed modelling the force between the subject and the back support of the test bench during a seatbelt test. R. Desbats is Ph.D. student in Biomechanics, F. Bermond (email: francois.bermond@ifsttar.fr, Tel. +33 4 72 14 23 78) is researcher, S. Nicolle is Associate Professor and P. Vezin is senior researcher at Université de Lyon -IFSTTAR. S. Compigne is technical manager at Toyota Motor Europe NV/SA.

Lumped element model of THOR-NT and PMHS abdomen under seatbelt and impactor loading

Romain Desbats, François Bermond, Sabine Compigne, Stéphane Nicolle, Philippe Vezin

The objective of this study is to continue the study from [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] with more loading conditions and to quantify the biofidelity of the THOR dummy through a comparative analysis of the responses to seatbelt and impactor loadings of the PMHS and the THOR dummy. A common lumped element model will be used to reproduce the PMHS and THOR responses by optimising the model's mechanical parameters. The use of a lumped element model allows to identify mechanical properties more quickly than a FE model, due to a shorter computation time and the reduced number of parameters. The differences between the PMHS and the THOR responses could be explained and quantified by the variations between parameters values.

II. METHODS

Test data

Seatbelt

The seatbelt tests selected from the literature and applied to the lumped element model are taken from [START_REF] Kent | Thoracic Response to Dynamic, Non-Impact Loading from a Hub, Distributed Belt, Diagonal Belt, and Double Diagonal Belts[END_REF][START_REF] Parent | Thoracic Biofidelity Assessment of the THOR MOD KIT ATD[END_REF][START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] for the PMHS tests and from [START_REF] Compigne | Proposal for a Modified THOR Lower Abdomen Including Abdominal Pressure Twin Sensors[END_REF] (replication of tests from [START_REF] Foster | High-Speed Seatbelt Pretensioner Loading of the Abdomen[END_REF]) and [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] for the THOR dummy. These are fixedback tests where the belt is placed on the abdomen at the umbilicus level and retracted toward the back of the subject thanks to pretensioners (or by a piston for one of two condition in [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF]). Three conditions from the literature are common to PMHS subjects and the dummy: B and C conditions from [START_REF] Foster | High-Speed Seatbelt Pretensioner Loading of the Abdomen[END_REF] (different pretensioners) and condition 2 from [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF].

For the PMHS tests, the average penetration across all the subjects for a given condition was applied to the model. The average force was compared to the model response. Two tests were performed for each condition in [START_REF] Compigne | Proposal for a Modified THOR Lower Abdomen Including Abdominal Pressure Twin Sensors[END_REF]. The average penetration of those two tests was computed, differentiated to obtain velocity and then differentiated a second time to determine acceleration. The three signals were filtered in order to avoid noise that would distort the model computations. Average penetration, velocity and acceleration were applied to the model and the average force was compared to the model response. One test per condition was performed in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF]. Fig. 1 shows those input conditions applied to the lumped element model overlaid together. It can be seen that all the input conditions do not have the same time duration and for some of them the data is not available over the full time range of the test.

Impactor

Six impactor conditions were applied to the model. For the PMHS subjects, two conditions were taken from [START_REF] Cavanaugh | Lower Abdominal Tolerance and Response. 30th Stapp Car Crash Conference Proceedings[END_REF] and [START_REF] Hardy | Abdominal Impact Response to Rigid-Bar, Seatbelt, and Airbag Loading[END_REF]. The tests with a 64 kg mass from [START_REF] Cavanaugh | Lower Abdominal Tolerance and Response. 30th Stapp Car Crash Conference Proceedings[END_REF] have been excluded. For the THOR dummy two conditions were taken from [START_REF] Compigne | Proposal for a Modified THOR Lower Abdomen Including Abdominal Pressure Twin Sensors[END_REF]. There is only one common configuration for the PMHS and the dummy, the 6 m/s (nominal velocity) 32 kg test condition. For each configuration, the penetration data from the different tests were averaged and the corresponding initial velocity was taken as input condition 0 v . Table I shows the considered input velocities. This mass is fixed to model a seatbelt loading and let free to model an impactor test. An important modification of the structure of the model was carried out in order to represent both seatbelt and impactor loading cases.

Indeed, the original model from [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] had the springs and dampers in parallel (Kelvin-Voigt model) and therefore could model a seatbelt test but not an impactor test. Since the extremity of the model has a non-zero initial velocity, the fact of having a damper in parallel created a non-zero initial reaction force. This is in contradiction with the test data. Therefore, it was chosen to transform the front stage of the model into a spring in series with a damper (Maxwell model). Fig. 2 shows the models for seatbelt and impactor loading cases.

Model equations and resolution

Seatbelt

The variables of the model can be seen on Fig. 2

( ) ( ) ( ) (b) 0 (a) 1 2 2 1 2 1 1 1 1 1 x x C x x K x x C x C x K x - ⋅ + - ⋅ - = - ⋅ - ⋅ + ⋅ = ⋅ - (1)
Equations 1(a) and 1(b) are two coupled differential equations where the unknowns are 1

x and 2

x and where 2

x is used to compute F from Equation 2 (obtained from isolating the mass m ) and 1

x is used to compute the force between the back of the subject and the test bench, as shown on Equation 3.

( )

2 x x K x m F - ⋅ + ⋅ = & &
(2)

1 1 1 1 back x C x K F & ⋅ + ⋅ = (3)

Impactor

The model for impactor case is detailed on Fig. 2(b). The initial impactor velocity, ( )

0 0 = = t x v &
, is imposed to the model. Equations 4(a) to 4(d) give the equations of motion of the model. The interaction force between the impactor and the subject is computed according to Equation 5. The differential equations systems are solved with a Runge-Kutta method programmed in the Scilab software.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (d) (c) (b) 0 (a) 3 1 1 3 1 1 3 b 3 1 1 3 1 1 1 2 1 1 2 2 2 i x x C x x K x m x x C x x K x x C x -M x x C x x K x x K x m m - ⋅ - - ⋅ - = ⋅ - - ⋅ + - ⋅ + - ⋅ - = ⋅ - ⋅ + - ⋅ - = - ⋅ = ⋅ + - (4) ( ) x m m F & & ⋅ + = i (5)
Cubic non-linear element as used in [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] proved to be difficult to use in this study: a cubic stiffness has little effect on small compression values but great effect on higher compression values due to the nature of the cubic relation. It was therefore difficult to fit all the test data with such elements. Good non-linear components would include more than one parameter in order to amplify differently low and high compression values. It was therefore decided to use only linear components in the model, if they are sufficient to fit the test data, so as not to add parameters to the model if not necessary. However, it was not possible to fit the THOR dummy impactor test data with linear elements. It was therefore decided to replace the element K with a non-linear spring that shows a strain dependence as seen in Equation 6. The nonlinear stiffness browses successively two regimes of compression and is characterised by three parameters: a linear stiffness 0 K ; a compression limit d ; and a power p (imposed as an integer number). For low compression values, when (

)

d x x << -2
, the term between curly brackets in Equation 6 reduces to unity and the response of the spring is linear, ( )

2 0 spring x x K F - ⋅ =
. At moderate to high compression, a stiffening sets in and the restore force becomes higher than the 0 K linear contribution.

( )

                  - + ⋅ - ⋅ = p d x x x x K F 2 2 0 spring 1 (6)

Determination of model parameters

Masses determination

Due to the abdomen representation described above, the masses of the model had to be estimated for the human subjects. It was chosen to use the THUMS model described in [START_REF] Shigeta | Development of Next Generation Human FE Model Capable of Organ Injury Prediction[END_REF] to measure those masses. The THUMS model has a representation of the abdomen internal geometry and densities for a 50 th percentile male , where BM is the average body mass of the PMHS subjects of the considered test condition. The masses for the THOR dummy were measured on the dummy v1.0 FE model (see [START_REF] Untaroiu | Evaluation of a Finite Element Model of the THOR-NT Dummy in Frontal Crash Environment[END_REF] for v0.0). The part densities were taken from [START_REF] Gesac | Test Support for Finite Element Modeling of THOR Crash Test Dummy[END_REF]. The mass of the two foam blocks of the abdomen together is approximately 0.3 kg. Therefore, M was set to 0.2 kg and m was set to 0.1 kg. The mass m should also account for a portion of the dummy jacket, but it is not possible to estimate this contribution.

The mass b m was computed as b m = 78.3 kg m M --= 78 kg, since the total mass of the dummy is equal to 78.3 kg according to [START_REF] Gesac | Thor NT Advanced Frontal Impact Dummy Internet: www.gesacinc[END_REF].

Optimisation loop

The response of the model for a specific loading condition is computed with a Scilab program. Based on manually adjusted initial values, the parameters values are optimised for each test configuration with the Scilab "optim" function until the model response (the force signal) matches the test data, minimising the criteria

( ) ( ) ( ) ∑ = - = N i i F i F f 1 2 test model
with N the number of data points. The goodness of fit is then estimated by calculating cross-correlation coefficients and the variation range of the parameters is determined.

Goodness of fit assessment

The assessment of the goodness of fit between the model response and the reference test data is obtained by = m front flesh = M abdomen content back with rest of the body calculating amplitude ratio, shape factor and phase shift as described in [START_REF] Xu | Repeatability Evaluation of the Pre-Prototype NHTSA Advanced Dummy Compared to the Hybrid III[END_REF]. Equations 7(a) to 7(c) define the three ratios with ( ) t x being the signal to compare to the reference ( )

t y and ( ) ∫ +∞ ∞ - ⋅ = t t x x d 2 the norm of x .
The integrals of digital signals are computed as:

( ) ∑ ∫ - = - +∞ ∞ - + ⋅ = ⋅ 1 1 1 2 d N i i i x x T t t x
with T the sampling rate (numerical integration by trapezoidal rule). A perfect fit between the two signals would lead to an amplitude ratio and a shape factor of 1 and a phase shift of 0.

( ) ( )

( ) ( ) (c) minimal is d that so shift phase (b) d factor shape (a) ratio amplitude y x t h t y t x h y x t t y t x y x ⋅ ⋅ + ⋅ = ⋅ ⋅ ⋅ = = ∫ ∫ ∞ + ∞ - ∞ + ∞ - (7)

Sensitivity analysis

In order to compare the model parameters for different conditions, it was necessary to estimate the range in which the optimised parameters can vary without significantly changing the model response result. The optimised values of the model parameters (spring stiffnesses and damping parameters) and the masses values were subjected to variations by steps of 10% of their value, ranging from 0.1 to 2 times the initial value. For each parameter value the three correlation coefficients mentioned above between the model response and the test data were computed. The percentage of variation of the three coefficients was computed at each step. A variation of one of the coefficients exceeding 10% of the coefficient value computed with the optimised parameter was considered as the limit of the range of variation of the parameter. An exception was made for the phase shift coefficient. In order to avoid the artefacts due to small phase values, the value of phase shift used in these computations was divided by 10 ms and the coefficient ms 10 shift phase 1was considered instead of the phase shift in ms, in order to have a coefficient equal to 1 in case of a perfect fit. For the seatbelt conditions, the initial negative slope of the model force response for some tests is an artefact due to the negative slope in the x & data (see Fig. 1(b)). This may be a consequence of the filtering of penetration data to compute the input velocity. The THOR dummy response shows sharp peaks and an abrupt initial slope that the model can-not fit exactly, regardless of the potential time offset in the force signal from test data. For the impactor conditions, The non-linearity of the dummy response can be seen on Fig. 5(c For both seatbelt and impactor case the model allows to fit the test data, the response being mainly in the standard deviation corridors for the PMHS, the amplitude and shape ratios stay between 0.8 and 1.1 and the phase shift between ±1.2 ms. This means that this model is suitable to analyse the human and dummy abdomen response.

III. RESULTS

Model fit to test data

Parameters values

Figure 6 shows the parameters identified for each common PMHS/dummy condition along with their range of variation for the seatbelt and impactor cases, respectively. The numerical values are reported in For the seatbelt case, the first thing to be noticed on Fig. 6(a) is that the K values for the PMHS subjects are much higher than for the THOR dummy. The C values are higher for the dummy than for the PMHS subjects (Fig. 6(c)), except for configuration 2 from [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] where they are equal. The 1 C values (Fig. 6(g)) are higher for the dummy compared to the PMHS, except for the configuration 2 from [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF], but since ranges of variation overlap for most of the test conditions, it is concluded that the 1 C values are the same for the PMHS and the dummy.

Additionally, the 1 K values presented on Fig. 6(e) are similar between the PMHS and the dummy, but large ranges of variation exist for all of the conditions. This high range of variation of the parameter 1 K is due to the fact that 1

x is significantly lower than x for most of the conditions. Therefore, a given variation of 1 K would have less effect on the abdomen force than the same variation of K . Regarding the results for configuration 2 from [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF], the fact that C is equal for the PMHS and the dummy and that 1 C is lower for the dummy, can be explained by the fact that the penetration velocity is the highest for this condition (Fig. 1(b)).

For the impactor case, since a non-linear spring was needed to fit the THOR test conditions, the effective stiffness of these conditions is compared to the K values of the other conditions on Fig. 6 x & for seatbelt,

IV. DISCUSSION

Interpretation of the results

The parameters values between seatbelt and impactor cases can be compared on Fig. 6. The K values for the PMHS subjects (Fig. 6 K values match between seatbelt and impactor cases for the THOR dummy, whereas values are much lower in the impactor case for the PMHS (Fig. 6(e) and (f)). However, the large range of variation (relative to the optimised value) for the PMHS in the impactor case allows 1 K to take any value to fit the data. 1 C

values are similar for PMHS and the THOR dummy (Fig. 6(g) and (h)) between the seatbelt and impactor cases and some large ranges of variation exists.

In order to understand the meaning of the parameters value comparison, one has to note that the influence of the parameters is not the same whether they are placed in series or in parallel in the model. For instance, a spring in parallel will have no effect when its stiffness is equal to 0, whereas a spring in series has no effect when its value tends toward infinity. It is the same for damping parameters. This explains the important range These findings correlate with the test data reported in [START_REF] Compigne | Proposal for a Modified THOR Lower Abdomen Including Abdominal Pressure Twin Sensors[END_REF]. For the impactor case the dummy force response was higher than the PMHS responses and for the seatbelt case the dummy force response was lower than the PMHS responses. This is explained by the previous findings. An impactor loading creates high abdomen penetration, therefore the dummy force response is higher than the PMHS one due to the higher influence of K (the spring non-linearity enhances this phenomenon since high penetration increases the spring force contribution). On the other hand, a seatbelt loading creates less penetration (compared to the impactor loading) and a sharp velocity peak, therefore the dummy force response is lower than the PMHS one due to the lower influence of C .

The fact that the K values are extremely different between the dummy and the PMHS subjects only for the seatbelt case is also explained by the PMHS mainly viscous behaviour for the seatbelt condition, whereas for the impactor condition the spring and damper contribution are more balanced. This may be due to the fact that for the seatbelt condition, the displacement and velocity of the extremity of the model is imposed. The force profile to fit having a similar shape than the velocity profile explains that the K value has less influence. The results would perhaps be different if the force were imposed. In the impactor case, only the initial velocity of the model's extremity is imposed, which imposes less constraints on the model. This also explains why the PMHS K values are higher in the seatbelt case as compared to the impactor case.

This spring-driven behaviour of the dummy and the damper-driven behaviour of the PMHS highlighted by the model are due to the nature of the constitutive material of the human subjects and the dummy. The human abdomen is made of soft tissues that include a large proportion of water, therefore causing a viscous behaviour whereas the dummy abdomen is made of foam and therefore has an elastic behaviour, although not linear and strain rate dependent. The influence of the dummy spine and pelvis (more important than for the human) explains the need for a non-linear spring in the impactor case. The optimised values of the d parameter of the non-linear spring (42 mm and 52 mm) represent the level of compression after which the non-linearity appears. This corresponds to the depth the abdomen foam has to be compressed before the pelvis is involved (measured at approximately 45 mm on the dummy).

Comparison with other studies

References [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] and [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] also developed lumped element models of the abdomen for seatbelt loading. It should be noted that in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF], the penetration of the seatbelt into the abdomen was estimated as the average of the left and right belt strand displacement, although in [START_REF] Foster | High-Speed Seatbelt Pretensioner Loading of the Abdomen[END_REF] and [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] the belt penetration in the abdomen was measured directly.

The stiffness parameters found in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] were from 12000 N/m to 17000 N/m (except for the subject PRT035

showing a very low stiffness) and the damping parameters from 600 N/m.s to 900 N/m.s. The lumped element model used for the THOR dummy in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] had an extra mass compared to the model used for PMHS and the spring was non-linear. It is difficult to compare those values with the K and C values of this study since the first stage of the present model, where the displacement is imposed, consists of two components in series and not in parallel as in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF]. It was concluded in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] that the THOR abdomen was too stiff (approximately two times)

and not viscous enough (approximately four times) compared to the human subjects. The same conclusion can be drawn in the present study for K and C , with the distinction that the components are in series in the present model, therefore it is more correct to say that the abdomen response is more influenced by the spring and less by the damper for the dummy compared to the PMHS. It was chosen to compare the parameters of the THOR dummy and the human subjects based on the same lumped element model until it was necessary to use a non-linear spring for the dummy under impactor loading.

Reference [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] applied a lumped element model to PMHS response under seatbelt loading. But as the front stage of the model had K and C in parallel it is difficult to compare their values with those from the present model. Non-linear elements were used for K , 1 K and 1 C , which also makes it difficult to compare them with the values from this study. Moreover, in [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] the force from test data was applied to the model and the penetration given by the model was compared with the test data. This is a different approach from the one used in this study, where the penetration was applied to the model, as in [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] and [START_REF] Kent | Thoracic Response to Dynamic, Non-Impact Loading from a Hub, Distributed Belt, Diagonal Belt, and Double Diagonal Belts[END_REF], which is believed to be more representative of a seatbelt loading with a piston or pretensioners. The two different methods could potentially lead to different results with the same parameters, especially in cases where the elements of the model are non-linear.

There is also a difference in the manner of determining the value of the parameters for given reference data. In [START_REF] Trosseille | Abdominal Response to High-Speed Seatbelt Loading[END_REF] and [START_REF] Lamielle | 3D Deformation and Dynamics of the Human Cadaver Abdomen under Seatbelt Loading[END_REF] the parameters were identified from the test data as sloped in a force-penetration diagram for a stiffness parameter, for example. In this study an optimisation procedure was used, as in [START_REF] Lobdell | Impact Response of the Human Thorax[END_REF] and [START_REF] Parent | Thoracic Biofidelity Assessment of the THOR MOD KIT ATD[END_REF]. A sensitivity analysis of the optimised parameters was also performed to predict the effect of the variation of the parameters on the response of the model. would account for those limitations.

Limitations

As described above, the second stage of the model contributes less to the global penetration response of the model. This causes the 1 K and 1 C parameters to have less influence on the response and more variability than the K and C parameters. This limits the relevance of a two stages model.

The fact that for the PMHS subjects in the impactor case the model penetration is higher than the penetration from test data raises a limitation of the model. The penetration is not imposed on the model, unlike for the seatbelt case, and since the model parameters are optimised to fit the force from test data only, the penetration response of the model does not necessarily fit the test data. This could be due to the fact that the model does not take into account the abdomen depth, therefore allowing the predicted penetration to be higher than the test penetration. This could also be due to the Maxwell model used for the first stage of the model, where the C damping element does not unload. This can be seen on

Guidelines for improving the THOR dummy

Although there are differences in the model parameters values between the seatbelt and impactor case, the results show that the THOR dummy abdomen is more elastically and less viscously deformable compared to the human abdomen. This behaviour is mainly explained by the first stage of the model, i.e. related to parameters K and C since the second stage contributes less to the response as explained above.

Based on these observations, some improvements can be applied to the dummy at the material and structural levels. The material properties of the lower abdomen foam should be modified to favour the viscous behaviour instead of the elastic behaviour. The need of a non-linear spring for the dummy model at high penetrations, whilst it was not necessary for the PMHS, means that stiffer parts (such as the pelvis) contribute non-biofidelically more to the response than soft parts. This confirms the design changes of this region implemented on THOR Mod Kit, which shortened the pelvis flesh at antero-superior iliac spines by around 20 mm. Furthermore, the difference in moving mass between the dummy abdomen and the PMHS has an influence and therefore the dummy abdomen mass should be increased.

V. CONCLUSION

A lumped element model has been designed in order to reproduce the abdomen force-penetration response of human subjects and the THOR dummy under seatbelt and impactor loading. This model proved to be capable of fitting the different test data and allows an explanation of the differences between the abdomen response of human subjects and a crash test dummy in terms of mechanical parameters. The viscous contribution to the dummy abdomen deformation has to be increased relative to the elastic contribution in order to improve its biofidelity. To achieve this goal, the material and structural modifications of THOR abdomen proposed in this study will be tested on THOR Metric FE model released by NHTSA in 2015 (v2.1). The ultimate goal will be to propose a new abdomen prototype.

VIII. APPENDICES

Fig. 1 .

 1 Fig. 1. Input conditions for seatbelt loading conditions common to the PMHS and the THOR dummy. B: B condition from Foster et al. 2006 C: C condition from Foster et al. 2006 2: configuration 2 from Trosseille et al. 2002

Fig. 2 .

 2 Fig. 2. Simplified abdomen model for seatbelt and impactor loading cases.

(Fig. 3 .

 3 Fig. 3. Masses repartition of the THUMS model.

Figures 4

 4 Figures 4 and 5 show the model response fit to the test data under seatbelt and impactor loading respectively, for the common conditions between PMHS subjects and the dummy. This model response is obtained with the optimised mechanical parameters values reported in Tables A.I. and A.II., respectively (see Appendix). Nonnormalised data have been used for the PMHS penetration and force response. The standard deviation corridors are indicated for the PMHS force data. The goodness of fit coefficients are reported at the bottom of each subfigure. Force-penetration graphs of those results are presented in Fig. A.1. (see Appendix).

  Figures A.1. (c) and (d) for the PMHS data, it appears that the model predicts a higher penetration than measured in the test data.

F 1 K

 1 (b). The effective stiffness for the non-linear spring is computed as taken from Equation6, the force given by the spring. ff e K is plotted against the relevant elongation, 2 x x -. For a linear spring (PMHS case), K K e = ff . For clarity, only optimised values are plotted, without range of variation. The ff e K values for the THOR dummy are lower than the K values of the PMHS subjects on most of the penetration range. The C values are higher for the THOR dummy compared to the PMHS subjects (Fig. 6(d)) except for 10 m/s condition where the PMHS have a higher C value than the other conditions. The same trend is noted for the 1 C parameter (Fig. 6(h)), although there is a large range of variation for the 10 m/s condition. The 1 K values presented on Fig. 6(f) seem higher for the dummy compared to the PMHS, but due to relatively large ranges of variation of this parameter it can take the same values for all the conditions. The same observation can be made for the high range of variation of than in the seatbelt case. The fact that the 10 m/s condition has a higher C value compared to the other PMHS condition can be due to the fact that this is the highest velocity of the test conditions. It stands out that the second stage of the model contributes less to the global penetration response of the model since the displacement and velocity of the second stage, ( 1 x and 1

  impactor), are lower than those of the front stage ( -). This leads to focussing on the K and C parameters for the analysis.For the conditions that are not common to the PMHS and the dummy, the optimised parameters are presented in Fig. A.2. (see Appendix). The same trends and order of magnitudes as explained above can be noticed.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Fit of the model response to test data for seatbelt case. B: B condition from Foster et al. 2006 C: C condition from Foster et al. 2006 2: configuration 2 from Trosseille et al. 2002

(

  (a) and (b)) are very different between the seatbelt and impactor cases. Regarding the THOR dummy, the K values for the seatbelt case are the same order of magnitude as the effective stiffness values of the impactor case. The origin value of the THOR curves on Fig. 6(b) is the 0 K value. The linear stiffness of the dummy modelled with a non-linear spring is lower than the stiffness of the other conditions. The C values (Fig. 6(c) and (d)) match well between the seatbelt and impactor cases (except for the 10 m/s condition). 1

  of variation toward high values for K and C concerning PMHS and THOR respectively. A high value for such a component in series means it deforms little. Therefore a higher value means it deforms even less, which does not necessarily change the global result. The model in series has a characteristic time K C = τ . The model behaviour is elastic for times inferior to τ and viscous for times superior to τ . The characteristic times have been computed with the optimised K and C values. It can be concluded that the penetration response of the abdomen is mainly caused by the instantaneous deformation of the spring K for the THOR dummy and by the long-term deformation of the damper C for the PMHS. For the seatbelt case, this is highlighted by Figures A.3. (a) and A.3. (b) (see Appendix), which show the prominence of the K spring and C damper responses on the abdominal penetration of the dummy and the PMHS, respectively. The PMHS characteristic time being low (less than 3 ms) compared to the loading duration, the PMHS response is mainly viscous. On the opposite, the dummy characteristic times are between 20 ms and 80 ms, higher than the loading duration, therefore the response is mainly elastic. Figures A.4. (a) and A.4. (b) (see appendix) show the same phenomenon for the impactor case. It is difficult, however, to compare the characteristic times with a non linear spring in the dummy case.

  Lumped element models applied to impact biomechanics have inherent limitations. They only represent a one-dimensional loading and use simple mechanical parameters. A linear model proved to be sufficient to fit the PMHS and dummy response until non linearity appeared for dummy impactor tests. The structure chosen for the first stage of the present model (Maxwell model) theoretically represents a fluid behaviour. If loaded with an imposed constant force (like at the end of a seatbelt loading, for instance), it would result in an infinitely increasing penetration. A generalised Maxwell model with three elements (standard linear solid)

  Fig. A.4. (b) where the compression of C does not decrease during unloading. It would also be useful to have abdomen depth considered in the lumped element model to account for the individual anthropometry variations between PMHS subjects.

Fig. A. 1 .

 1 Fig. A.1. Force-penetration responses. B: B condition from Foster et al. 2006 C: C condition from Foster et al. 2006 2: configuration 2 from Trosseille et al. 2002 PMHS 6.1 m/s: 32 kg 6.1 m/s condition from Cavanaugh et al. 1986 PMHS 10 m/s: 32 kg 10 m/s condition from Cavanaugh et al. 1986 THOR 3 m/s: 32 kg 3 m/s condition from Compigne et al. 2015 THOR 6.1 m/s: 32 kg 6.1 m/s condition from Compigne et al. 2015

Fig. A. 3 .

 3 Fig. A.3. Displacement results for seatbelt loading conditions. B: B condition from Foster et al. 2006 C: C condition from Foster et al. 2006 2: configuration 2 from Trosseille et al. 2002

Fig. A. 4 .

 4 Fig. A.4. Displacement results for impactor loading conditions. PMHS 6.1 m/s: 32 kg 6.1 m/s condition from Cavanaugh et al. 1986 PMHS 10 m/s: 32 kg 10 m/s condition from Cavanaugh et al. 1986 THOR 3 m/s: 32 kg 3 m/s condition from Compigne et al. 2015 THOR 6.1 m/s: 32 kg 6.1 m/s condition from Compigne et al. 2015

TABLE I

 I

	INITIAL IMPACT VELOCITY VALUES 0 v (m/s) m (kg) i Cavanaugh et al. 1986 6.1 m/s condition 6.1 32 Cavanaugh et al. 1986 10 m/s 10 32 Hardy et al. 2001 6 m/s 6.5 48 Hardy et al. 2001 9 m/s 8.5 48 Compigne et al. 2015 3 m/s 3.0 32 Compigne et al. 2015 6.1 m/s 6.3 32 (a) Abdomen penetration. PMHS THOR (b) Abdomen penetration velocity. kinetic energy (J) 595 1600 1028 1745 141

  (a). It was chosen to impose the x displacement and its derivatives x & and x& & as input conditions. The model response was considered as the force F . They are taken from the test data and imposed on the model. Equations 1(a) and 1(b) are the equations of motion of the system, obtained by isolating the mass M and the point where 2 x is measured (massless node linking the spring and the damper in series), respectively.

Table A .

 A Table A.I. and II. (see Appendix).

TABLE A

 A

				Trosseille et al. 2002 condition 2
			K	C	K	1	1 C	M	m
		superior limit	137326	3178	43387	428	0.4	0.2
		optimised value	124841	2445	39443	389	0.2	0.1
		inferior limit	112357	1956	35499	350	0.02	0.01
				Foster et al. 2006 B condition
	THOR	superior limit	K 125792	C 19356	1 62396 K	C 1 1000	M 0.4	m 0.2
		optimised value	114357	9678	51997	910	0.2	0.1
		inferior limit	102921	6774	41597	819	0.02	0.01
				Foster et al. 2006 C condition
			K	C	K	1	1 C	M	m
		superior limit	104873	4986	55076	2317	0.4	0.2
		optimised value	95339	4155	27538	2106	0.2	0.1
		inferior limit	85805	3324	2754	1896	0.02	0.01
				.I.
		IDENTIFIED PARAMETERS FOR SEATBELT LOADING CONDITIONS, INCLUDING SENSITIVITY ANALYSIS
		(PARAMETERS VALUE RANGE WHICH KEEPS GOODNESS OF FIT PARAMETERS WITHIN ±10% RANGE)
				Trosseille et al. 2002 condition 2
			K (N/m)	C (N/m.s)	1 K (N/m)	1 C (N/m.s)	M (kg)	m (kg)
		superior limit	1350054	1217	210787	3282	10.2	1
		optimised value	900036	1217	105393	2344	6	0.5
		inferior limit	630025	1217	10539	1641	1.8	0.05
				Foster et al. 2006 B condition
			K	C	K	1	1 C	M	m
	PMHS	superior limit	1140043	1107	76428	1527	9.6	1.6
		optimised value	570021	1107	38214	1175	6	1
		inferior limit	342013	1107	3821	940	3	0.3
				Foster et al. 2006 C condition
			K	C	K	1	1 C	M	m
		superior limit	1692911	330	77983	1608	12	1.1
		optimised value	846456	300	38992	804	6	1
		inferior limit	423228	270	3899	80	1.2	0.9

TABLE A .

 A II. IDENTIFIED PARAMETERS FOR IMPACTOR LOADING CONDITIONS, INCLUDING SENSITIVITY ANALYSIS (PARAMETERS VALUE RANGE WHICH KEEPS GOODNESS OF FIT PARAMETERS WITHIN ±10% RANGE)

		(a) Seatbelt test data.				(b) Seatbelt model results.
					Cavanaugh et al. 1986 32 kg 6.1 m/s condition
					K		C	K	1	1 C		M	m	m	b
					(N/m)		(N/m.s)	(N/m)	(N/m.s)	(kg)	(kg)	(kg)
		superior limit			31179		2244	3949	1189		9	2	110
		optimised value			23984		1726	1974	991		6	1	55
		inferior limit			21585		1381	197	793		3	0.1	33
					Cavanaugh et al. 1986 32 kg 10 m/s condition
					K		C	K	1	1 C		M	m	m	b
		superior limit			66570		25217	4149	2746	10.8	2	154
		optimised value			55475		12609	2075	2288		6	1	77
	PMHS	inferior limit			49928		6304	207	1830		1.2	0.1	53.9
						Hardy et al. 2001 48 kg 6 m/s condition
		(c) Impactor test data.	K		C	1 (d) Impactor model results. K 1 C	M	m	m	b
		superior limit			37220		1384	81092	5847		12	2	80.4
		optimised value			37220		1259	40546	2923		6	1	67
		inferior limit			37220		1133	4055	1462		0.6	0.1	60.3
						Hardy et al. 2001 48 kg 9 m/s condition
					K		C	K	1	1 C		M	m	m	b
		superior limit			73832		3338	20284	1873		9.6	2	74.1
		optimised value			67120		2781	10142	1561		6	1	57
		inferior limit			60408		2225	1014	1249		1.8	0.1	45.6
						Compigne et al. 2015 32 kg 3 m/s condition
			K	0	d	p	C	K	1	1 C		M	m	m	b
			(N/m)	(mm)	(no unit)	(N/m.s)	(N/m)	(N/m.s)	(kg)	(kg)	(kg)
		superior limit	5284	52	4	11588	32005	3326		0.4	0.2	101.4
		optimised value	4804	52	4	6817	16002	2772		0.2	0.1	78
	THOR	inferior limit	4323	52	4	4772	1600	2217	0.02	0.01	62.4
					Compigne et al. 2015 32 kg 6.1 m/s condition
			K	0	d	p	C	K	1	1 C		M	m	m	b
		superior limit	6911	42	4	9514	31955	3519		0.4	0.2	109.2
		optimised value	4936	42	4	6796	15977	3199		0.2	0.1	78
		inferior limit	3949	42	4	5437	1598	2879	0.02	0.01	62.4

& & & & & & &

& & & & & & & & & & & & & &

VI. ACKNOWLEDGEMENTS

The authors would like to thank Messrs Philippe Petit and Xavier Trosseille for allowing to use data from Trosseille et al. 2002 and from Lamielle et al. 2008 as well as Mr Craig Foster for providing the data from Foster et al. 2006.


B Abdallah

• F Willot

• D Jeulin

Stokes flow through a Boolean model of spheres: Representative Volume Element

Keywords: Porous media, Stokes flow, FFT methods, Representative Volume Element, Tortuosity, Streamlines

). The latter is applied to large system sizes representative of the microstructure. The overall permeability of the two models as well as the Representative Volume Element (RVE) are estimated as a function of the pore volume fraction. We give numerical estimates for the asymptotic behavior of the permeability in the dilute limit for the solid phase, and close to the percolation threshold of the pores. FFT maps of the velocity field are presented, for increasing values of the pore volume fraction. The patterns of the local velocity field is analysed using various morphological criteria. The tortuosity of the streamlines is found to be much higher than the geometrical tortuosity, for both models. The histograms of the velocity field are computed at increasing pore volume fraction. The covariance of orientation is used to characterize the spatial correlation of the velocity field.

Introduction

Since the pioneering works of [START_REF] Carman | Fluid flow through granular beds[END_REF], the permeability of random media has been a long-standing problem. In some early works, estimates for the permeability of flow around spheres have been derived making use of periodic structures [START_REF] Richardson | The sedimentation of a suspension of uniform spheres under conditions of viscous flow[END_REF][START_REF] Martin | Pressure drop through stacked spheres. Effect of orientation[END_REF] or of concentric sphere models [START_REF] Happel | Viscous flow in multiparticle systems: slow motion of fluids relative to beds of spherical particles[END_REF][START_REF] Kuwabara | The forces experienced by randomly distributed parallel circular cylinders or spheres in a viscous flow at small reynolds numbers[END_REF]. Other studies have focused on dilute, random suspensions of spheres [START_REF] Childress | Viscous flow past a random array of spheres[END_REF][START_REF] Howells | Drag due to the motion of a Newtonian fluid through a sparse random array of small fixed rigid objects[END_REF][START_REF] Hinch | An averaged-equation approach to particle intercations in a fluid suspension[END_REF] or on periodic arrays [START_REF] Zick | Stokes flow through periodic arrays of spheres[END_REF][START_REF] Sangani | Slow flow through a periodic array of spheres[END_REF][START_REF] Cheng | Flow past periodic arrays of spheres at low Reynolds number[END_REF]. Another approach, based on variational formulations, consists in deriving rigorous bounds for permeability [START_REF] Doi | A new variational approach to the diffusion and the flow problem in porous media[END_REF][START_REF] Rubinstein | Flow in random porous media: mathematical formulation, variational principles, and rigorous bounds[END_REF][START_REF] Torquato | Rigorous bounds on the fluid permeability: Effect of polydispersivity in grain size[END_REF]. These bounds take into account two-point correlations functions as well as the so-called "surfacesurface" and "surface-void" correlation functions. Interestingly, earlier bounds made use of three-point correlation functions [START_REF] Prager | Viscous flow through porous media[END_REF][START_REF] Weissberg | Viscous flow through porous media. III. Upper bounds on the permeability for a simple random geometry[END_REF]. More recently, self-consistent schemes for the effective permeability have also been proposed [START_REF] Boutin | Study of permeability by periodic and self-consistent homogenisation[END_REF][START_REF] Koo | Effective-medium theories predicting hydrodynamic transport properties of bidisperse suspensions[END_REF]. We refer to [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF][START_REF] Dormieux | Microporomechanics[END_REF] for an overview on the subject.

Among numerical methods for the computation of the effective permeability, finite-difference schemes [START_REF] Martys | Universal scaling of fluid permeability for sphere packings[END_REF] and, more recently, fast multipole expansions [START_REF] Koo | Effective-medium theories predicting hydrodynamic transport properties of bidisperse suspensions[END_REF] as well as Boltzmann Lattice Methods [START_REF] Belov | Modelling of permeability of textile reinforcements: lattice Boltzmann method[END_REF][START_REF] Lee | Effects of specific surface area and porosity on cube counting fractal dimension, lacunarity, configurational entropy, and permeability of model porous networks: Random packing simulations and NMR microimaging study[END_REF] have been used. Recently, largescale computations have been carried out using a finite element method [START_REF] Karim | Permeability determination of porous media using large-scale finite elements and iterative solver[END_REF]. Interestingly, several Fourier-based algorithms have been devised to treat Stokes flow [START_REF] Wiegmann | Computation of the permeability of porous materials from their microstructure by FFF-Stokes[END_REF][START_REF] Monchiet | A FFT-based method to compute the permeability induced by a Stokes slip flow through a porous medium[END_REF]N'guyen et al, 2013;[START_REF] Bignonnet | FFTbased bounds on the permeability of complex microstructures[END_REF]. Such Fourier-based numerical methods solve periodic Stokes flow problems, which are based on a rigorous homogenization theory [START_REF] Ene | Équations et phénomènes de surface pour lécoulement dans un modele de milieu poreux[END_REF][START_REF] Matheron | Remarques sur la loi de Darcy[END_REF][START_REF] Matheron | Genèse et signification énergétique de la loi de Darcy[END_REF]. The FFT method of [START_REF] Wiegmann | Computation of the permeability of porous materials from their microstructure by FFF-Stokes[END_REF] has been used to predict the effective permeability of open foam structures, using 3D images [START_REF] Redenbach | Modelling a ceramic foam for filtration simulation[END_REF]. However, the local velocity and pressure fields occuring in porous random microstructures has not been extensively studied in 3D. The "Representative Volume Element" (RVE) has been tackled elasticity and conductivity for polycrystals [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF], Boolean models [START_REF] Willot | Elastic behavior of composites containing boolean random sets of inhomogeneities[END_REF] and multiscale materials [START_REF] Willot | Elastic and electrical behavior of some random multiscale highly-contrasted composites[END_REF], making use of Matheron's theory of regionalized variables [START_REF] Matheron | The theory of regionalized variables and its applications[END_REF]. Other recent works include optics [START_REF] Azzimonti | Optical properties of deposit models for paints: full-fields FFT computations and representative volume element[END_REF] and acoustic [START_REF] Peyrega | Estimation of acoustic properties and of the representative volume element of random fibrous media[END_REF]. Results are also provided in [START_REF] Du | On the size of representative volume element for Darcy law in random media[END_REF] for Stokes flow in a 2D hard-core models of discs, however the subject has not been studied to the same extent for permeability in random media. This work focuses on the slow motion of a fluid past a random set of obstacles. It is organized as follows. The Stokes equation and microstructure of interest are given in Sec. ( 2), where we also present the numerical Fourier-based method. The Representative Volume Element (RVE) for the permeability is studied in Sec. (4). The microstructure permeability is computed in the entire range of pore volume fraction in Sec. ( 5). The asymptotic behaviour in the dilute and percolation limit is studied in Sec. ( 6). The streamlines, their tortuosity, the velocity field's histograms and covariances are computed in Sec. (7).

Stokes Flow in porous media

We consider a domain V made of a porous (P ) and solid (S) phase. We note f = V V = |P |/|V | the pore volume fraction and 1 -f = 1 -V V = |S|/|V | the volume fraction of the solid phase. For simplicity, we take V = [-L/2; L/2] 3 , a cubic domain of side L aligned with the axis e k (k = 1, 2, 3) of a Cartesian coordinates system. We are interested in the steady-state flow of an incompressible Newtonian fluid in the pores. The fluid velocity vector u(x) is defined in the pores where it follows the Stokes equation:

µ∆u i (x) = ∂ i p(x), ∂ i u i (x) = 0 (x ∈ P ), ( 1
)
where µ is the fluid viscosity, measured in Pa•s, p(x) is the pressure field at point x, measured in Pa and ∆ is the Laplacian operator. At the pore-solid interface ∂S, the fluid satisfies the no-slip boundary condition:

u(x) = 0 (x ∈ ∂S). (2)
By convention, the domain of definition of the fluid is extended to the solid phase as: u = 0 (x ∈ S).

(3)

Along the boundary of the domain, periodic boundary conditions for the fluid velocity are applied:

u(x ± e k L) ≡ u(x). (4)
A macroscopic pressure gradient of intensity δP along the direction E is applied:

∂ i p(x) = δP E i , |E| = 1, (5)
where • stands for the average over the domain V . The boundary conditions for the pressure field specifies that its fluctuations p * are periodic [START_REF] Ene | Équations et phénomènes de surface pour lécoulement dans un modele de milieu poreux[END_REF]:

p * (x) = p(x) -δP (x • E) = p * (x ± Le k ). ( 6
)
The permeability tensor κ ij , measured in m 2 , is defined at the macroscopic scale by Darcy's law:

u i (x) = - κ ij µ ∂ j p(x) = - κ ij E j µ δP. (7)
In the remaining part of this work, we assume that the porous phase is macroscopically isotropic. By convention, we choose E = e 1 so that the permeability tensor κ is identified with its component κ = κ 11 and:

κ = - µ u 1 (x) δP . ( 8
)
The permeability is independent of the fluid viscosity.

Equations. ( 1)-( 6) are solved numerically using the iterative Fourier-based scheme of [START_REF] Wiegmann | Computation of the permeability of porous materials from their microstructure by FFF-Stokes[END_REF]. The velocity and pressure field are discretized on a 3D grid of voxels. The method consists in rewriting the Stokes flow problem as a set of four Poisson equations which are discretized by centered differences. Artificial forces are added along solid surfaces to enforce the no-slip condition (2). Each component of the local velocity vector is interpreted as the flow of the fluid across voxel faces. The gradient and Laplacian operators are discretized by finite differences between adjacent voxels. The system of equations is solved iteratively using the generalized minimal residual method for linear systems [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. Finite differences are computed in the Fourier domains using the shift operator. The method uses 3D microstructure images, and no meshing is necessary, which allows one to treat arbitrarily complex geometries.

The permeability of the medium is computed numerically by setting δP = 1 [Pa], averaging on the velocity field u 1 (x) and using Eq. ( 8). By convention, we also fix µ = 1 [Pa•s]. In [START_REF] Wiegmann | Computation of the permeability of porous materials from their microstructure by FFF-Stokes[END_REF], convergence is achieved using minimal residual method (MINRES). At each iterations, Fourier transforms are used to invert the Laplacian operator. In the limit of a very high number of iterations, the fluid velocity is zero inside the obstacles. In practice, iterations are stopped when the following two conditions are met:

max x∈S ||u(x)|| ≤ η 1 max x∈P ||u(x)||, (9) ||u(x)|| S ≤ η 2 ||u(x)|| P , (10)
where • S,P are averages over the solid and porous phases and η 1 = 10 -3 , η 2 = 10 -4 .

In the following, we also compare our numerical results with analytical bounds and estimates. We consider the Carman-Kozeny estimates [START_REF] Carman | Fluid flow through granular beds[END_REF][START_REF] Kozeny | Ueber kapillare leitung des wassers im boden[END_REF]:

κ ≈ κ CK = f 3 cγ 2 , ( 11
)
where γ = |∂S|/|V | is the specific surface area, or surface area density, measured in m -1 , and c is an empirical constant. This adjustable parameter takes into account the shape of particles in random media. For the Boolean model, we fix it to c = 5 [START_REF] Happel | Viscous flow in multiparticle systems: slow motion of fluids relative to beds of spherical particles[END_REF]. For comparison purpose, the upper-bound κ UB proposed in [START_REF] Doi | A new variational approach to the diffusion and the flow problem in porous media[END_REF] is considered:

κ UB = 2 3 ∞ 0 dh f 2 γ 2 F SS (h) - 2f γ F SP (h) + F P P (h) , (12)
where F SS , F SP and F P P are the solid-solid, fluid-solid and solid-solid surface correlation functions, respectively. They are readily computed using the method given in [START_REF] Doi | A new variational approach to the diffusion and the flow problem in porous media[END_REF]. The histograms P ui of the components of the velocity field read:

P ui (t) = 1 P P d 3 x δ(t -u i (x)). ( 13
)
They are numerically computed from full-fields FFT maps. The latter also provide the tortuosity τ of the fluid streamlines which is defined by:

τ = L f L . (14)
In the above, L is the distance between two opposite faces in the domain, and L f is the mean length of the fluid paths which connect two points located at opposite faces. For comparison purposes, we also estimate the geometrical tortuosity:

τ Geo = L Geo f L , (15)
defined from the length L Geo f of geodesics in the porous phase. The latter are the shortest paths in the pores that connect two opposite sides of the domain. Another quantity of interest is the static viscous tortuosity τ 0 [START_REF] Johnson | Tortuosity and acoustic slow waves[END_REF][START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF][START_REF] Perrot | Dynamic viscous permeability of an open-cell aluminum foam: Computations versus experiments[END_REF] related to the fluctuations of the fluid local velocity by:

τ 0 = u 2 u 2 , ( 16
)
with u 2 = u i u i and u 2 = u i u i . The tortuosity factor given in [START_REF] Matyka | How to Calculate Tortuosity Easily?[END_REF] is similarly defined as:

τ A,B Matyka = ||u|| u 1 . (17)
Finally, we estimate the correlation of orientations of the local velocity field from [START_REF] Altendorf | Stochastic modeling of a glass fiber reinforced polymer[END_REF] using the following morphological criterion, based on scalar products:

C α (h) = x∈V d 3 h (u(x) • u(x + h)) 2 ||u(x)|| 2 ||u(x + h)|| 2 . ( 18
)
3 Boolean model of spheres

In the remaining part of this work, is considered a Boolean model of spheres [START_REF] Matheron | Élements pour une theorie des milieux poreux[END_REF] of constant diameter D > 0. In this model, the sphere centers are located according to a Poisson point process in V , so that spheres are allowed to overlap. The model is parametrized by the density of spheres centers or equivalently by the spheres volume fraction. We denote 1 -q the spheres volume fraction and q the volume fraction of the embedding medium. Two types of models are considered where spheres are either "solid" or "porous". In model (A) flow occurs solely outside the spheres whereas in model (B) flow occurs inside spheres only. Therefore, we have f = q for model (A) and f = 1 -q for model (B). We denote the permeability of these models by κ A and κ B respectively.

The specific surface area γ is obtained from:

γ = -4 ∂C(h) ∂h h=0 , (19)
where C(h) is the covariance function [START_REF] Matheron | Élements pour une theorie des milieux poreux[END_REF]. In a Boolean model of spheres [START_REF] Matheron | Élements pour une theorie des milieux poreux[END_REF] we have for models (A) and (B), respectively:

C A (h) = q 1+ 3h 2D -h 3 2D 3 , C B (h) = 1 -2q + q 1+ 3h 2D -h 3 2D 3 . ( 20
)
For both models the specific surface area is:

γ = - 6q D logq (21)
Therefore the Carman-Kozeny estimates (Eqs. 11) for models (A) and (B) read, respectively:

κ CK-A = f D 2 180 log 2 (f ) , κ CK-B = f 3 D 2 180(1 -f ) 2 log 2 (1 -f ) . ( 22
)
In the following, the Boolean model is discretized over a 3D grid of L×L×L voxels, with increasing size L = 100, 256, 512 (Fig. 1). The spheres diameter D is equal to D = 15 voxels. This makes sure the particles ressemble spheres, as shown in the enlarged 3D view in Fig. (1). Hereafter, by convention, all lengths are rescaled by the spheres radius R, or equivalently, we take

D = 2R = 2 [m].
Accordingly, microstructures contain an increasing number of spheres when L increases. The entire range of pore volume fraction

f A c ≤ f ≤ 1 for model (A) and f B c ≤ f ≤ 1 for model (B) is considered, where f B c and f A c
are the percolation threshold of the spheres and of its complementary, respectively. The following approximate values are available: f B c ≈ 0.2895 [START_REF] Rintoul | Precise determination of the critical threshold and exponents in a three-dimensional continuum percolation model[END_REF] whereas f A c ranges from 0.0317 (van der Marck, 1996; Priour Jr, 2014) to 0.0540 (Jeulin and Moreaud, 2006) according to different authors.

Note that the spheres are reproduced periodically in all directions without further modification. In particular, some spheres are completely surrounded with the fluid in model A.

Representative Volume Element

The effective permeability is estimated by the apparent permeability κ, computed by averaging on the velocity field (7). This estimation gives the exact permeability only in the limit of very large system sizes L → ∞. In a volume V of finite size, the relative error of the estimate for the permeability is given by:

= 2D κ (V ) √ nκ ( 23
)
where D 2 κ (V ) is the variance of the apparent permeability over n random realizations of volume V . Alternatively, using a single large realization of size V , we compute an estimate of D 2 κ (V ) for all volumes V < V by considering disjoint subdomains V ⊆ V . The asymptotic expansion of D 2 κ (V ) when V V and V → ∞ reads [START_REF] Matheron | The theory of regionalized variables and its applications[END_REF]:

D 2 κ (V ) D 2 κ ∼ A κ 3 1 V - 1 V , (24)
where A κ 3 is the integral range for the permeability [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]). In the above, the term 1/V represents a lower-order correction compared to 1/V . The point variance D 2 κ = u 2 1 -u 1 2 is numerically estimated using the velocity fields u 1 (x). In practice, the term 1/V can be neglected and we have:

D 2 κ (V ) D 2 κ ∼ A κ 3 V . (25)
Likewise, fluctuations are encountered when estimating the pore volume fraction in a volume V . The latter scales as:

D 2 f (V ) f (1 -f ) ∼ A f 3 V , V A f 3 , (26)
where we refer to A f 3 as the integral range for the microstructure. For model (B) [START_REF] Matheron | Élements pour une theorie des milieux poreux[END_REF]: ) the apparent permeability in each subdomain is given by the mean of the velocity field u 1 over each region. Accordingly, when the volume V is small, the variance is computed over a large data set. Conversely, the variance D 2 κ (V ) corresponding to the maximum volume size V = V /8 is computed over a set of 8 values.

A f 3 (f ) = |h|≤2R f 1+ 3h 2D -h 3 2D 3 -f 2 f (1 -f ) d 3 h , (27)
κ UB κ CK-B κ CK-A κ A (L=256) κ A (L=512) κ B (L=256) κ B (L=512) f c =3.2%, 29.0%

Permeability of the Boolean models

The effective permeability is estimated for increasing pore volume fraction f in the range [f A c ; 1] for model (A) and [f B c ; 1] for model (B). Results are plotted in Fig. 2 in lin-log scale. As expected, the permeability is largest for the complementary model of spheres (B), with lower percolation threshold than in model (A). Both models show a convexity change in lin-log scale at about 65% for model (B) and 50% of pore volume fraction for model (A).

Furthermore, the estimate κ CK-A (shown in green in Fig. 2) is close to FFT predictions in the range 0.3 < f < 0.8. Outside of this region, the estimate differs noticeably from the FFT data. As reported in [START_REF] Happel | Viscous flow in multiparticle systems: slow motion of fluids relative to beds of spherical particles[END_REF], the Carman-Kozeny estimate (11) with c = 5 agrees well with experimental data in the range 0.26 < f < 0.8 for particles of various shapes maintained in a fixed position. This indirectly confirms the validity of the FFT computations for Stokes flow. Interestingly, the Carman-Kozeny estimate κ CK-B for model B (red curve) is close to our FFT predictions in a different range of porosity, namely 0.7 < f < 0.9. In this domain, the regions surrounding the spheres are mostly disconnected, and the solid phase is akin to a set of random fixed particules with non-spherical shape. In the dilute regime f → 1, as expected, the Carman-Kozeny estimate fails for all models.

FFT fields are also used to study the representative volume element for the permeability as a function of the pore volume fraction f . The integral range A κ 3 is estimated numerically by fitting the left-hand side of (25) with a scaling law D 2 κ (V ) ∼ 1/V . As expected, the latter is recovered when V is large enough. This is illustrated in Fig. 3 where the variance of the apparent permeability D 2 κ (V ) of model B with f = 0.68 is represented as a function of V . For a requested precision , the corresponding RVE size is obtained using ( 23) and the previously determined integral range.

In Fig. 4 the relative precision is fixed to = 5%. We set V = 3 and plot the RVE lengths A κ and B κ for the permeability in models (A) and (B) respectively, as a function of f . The microstructure RVE lengths A f and B f at the same relative precision are plotted for comparison. For both models, the RVE size for the permeability strongly increases when approaching the percolation threshold. Strong microstructural effects are especially important for model B. Using the discretization D = 15 voxels, a system with 1350 3 voxels is required to estimate the permeability κ B at 5% relative error when f ≈ 0.32. The RVEs for permeability appear to be significantly higher than in elasticity or thermal conductivity [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF][START_REF] Willot | Elastic behavior of composites containing boolean random sets of inhomogeneities[END_REF]. It is also useful to compare these results with numerical data obtained for a Stokes flow around randomly-distributed, non-overlapping discs [START_REF] Du | On the size of representative volume element for Darcy law in random media[END_REF]. The authors studied the convergence of the apparent permeability with respect to the size of the computational domain. Highest convergence rates were obtained at lowest porosity (50%). Conversely, the worst convergence rates were observed at highest porosity (80%). In 3D, the RVE for model A is quite insensitive to the porosity in the range 0.5 < f < 0.8, but increases with f in the dilute regime f > 0.8 (Fig. 4). The behavior observed in 2D for non-overlapping discs is accordingly qualitatively recovered in 3D for the Boolean model, but for a different range of porosity. Interestingly, [START_REF] Du | On the size of representative volume element for Darcy law in random media[END_REF]) link the increase of the RVE with that of the size of the channels located in-between obstacles.

Dilute limit and behaviour near the percolation thresholds

We first consider the behavior of the permeability in models (A) and (B) in the dilute limit f → 1, i.e. for a small volume fraction of obstacles. Numerical FFT data are represented in log-log scale in Fig. (5a) (black and white symbols). The latter are compared to the dilute expansions provided by [START_REF] Rubinstein | Flow in random porous media: mathematical formulation, variational principles, and rigorous bounds[END_REF] and [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]:

κ A ≈ κ A,dil = 2R 2 9(1 -f ) , κ B ≈ κ B,dil = 4R 2 9(1 -f ) log 2 (1 -f ) , f → 1 (28)
which are valid for spherical obstacles and voids, respectively. The first-order corrections above are plotted in V'

D 2 κ (V ) κ A 3 ' Fig. 3 Variance D 2 κ )
of the apparent permeability computed on non-overlapping subdomains of volume V , with respect to V , in number of voxels, in log-log scale. Dashed line: fit ∼ D 2 κ A κ 3 /V resulting from the scaling law (25).

RVE for ǫ = 5%

System size computed for the permeability as a function of the pore volume fraction f , for models (A) and (B). The relative precision is fixed to rel = 5%. Thick lines: RVE lengths A,B f for the microstructure at the same precision. Horizontal black line: length of the most representative employed microstructure with L 3 = 512 3 voxels. Vertical lines: percolation thresholds f A c and f B c . All lengths are given in units of spheres radius.

RVE length ℓ A κ /R, L = 256 RVE length ℓ B κ /R, L = 256 RVE length ℓ A κ /R, L = 512 RVE length ℓ B κ /R, L = 512 f A c = 3.17% f B c = 28.95% RVE length ℓ A f /R RVE length ℓ B f /R
using the results of [START_REF] Rubinstein | Flow in random porous media: mathematical formulation, variational principles, and rigorous bounds[END_REF]. The analytical expansion ( 28) is in good agreement with FFT data in the region 1 -f < 10 -2 . For 1 -f > 10 -2 , the effective permeability in model (A) is significantly lower than the dilute expansion and upper-bound. Similarly, the permeability of model (B) predicted by FFT computations is smaller than expansion (28) for 1 -f > 10 -2 . In the highly-dilute regime 1 -f < 10 -2 , the agreement between expansion (28) and FFT data for model (B) is not as good as for model (A). Remark however that particles in model (B) have highly elongated shapes so discretization effects should be important in this case. The obstacles surface/volume ratio is indeed, on average, |∂γ/∂q| → ∞ when q = 1 -f → 0 (Eq. 19).

κ UB κ B,dil κ A,dil κ A (L=256) κ A (L=512) κ B (L=256) κ B (L=512) 10 -2 10 -1 f-f c A, B
κ UB κ A (L=256) κ A (L=512) κ B (L=256) κ B (L=512) powerlaw fit κ A ~f3.2 Powerlaw fit κ B ~f2.0 (
The behavior of the permeability close to the percolation thresholds f → f Numerical fits indicate:

κ A ∼ 0.28R 2 (f -f A c ) 3.2 , κ B ∼ 0.18R 2 (f -f A c ) 2.0 . ( 29
)
Predictions from the literature give κ A ∼ (f -f A c ) 4 [START_REF] Feng | Transp. properties of continuum systems near the percolation threshold[END_REF] for model (A). Likewise, an exponent equal to 4 was numerically determined by [START_REF] Martys | Universal scaling of fluid permeability for sphere packings[END_REF] for model (A), using spheres of diameter 15 voxels in a grid of 100 3 voxels. In the same condition, with a larger grid we recovered an exponent of about 3.5. In the present work, the range of validity of the fits ( 29) is limited to at most one decade. Computations on larger systems are required to confirm the above results.

Numerical values for the exponents near the percolation threshold are rough estimates, due to the challenging numerical computations but indicate that the exponent in model (B) is significantly lower than in model (A), i.e. close to the percolation threshold, the permeability in model B increases at a faster rate when f increases. In this section, are given and discussed maps of the velocity and pressure fields.

The fields histogram and streamlines tortuosity are also studied.

7.1 Maps and histograms of the velocity field FFT field maps are represented in Tabs. ( 1) and (2) for models (A) and (B) respectively. The maps are 512×512 voxels, 2D sections along the plane (e 1 , e 2 ) where e 1 is oriented right to left and e 2 top to bottom. The norm ||u|| of the velocity vector as well as its components u i are shown at 4 fixed values of the pore volume fraction. The component u 1 is parallel to the applied pressure field whereas u 2 and u 3 are transverse components. To compare the field maps for different values of the pore volume fraction, the colour scale is given and the values are indicated on the corresponding side of the field. Blue and red correspond to the lowest and highest values respectively, with green and yellow in-between, whereas the solid phase is coloured in grey. To highlight the fields patterns, each map has been thresholded at 75% of its maximum values. The maximum and minimum values of the velocity field, prior to the thresholding, is given in small font along the colour scale. At lowest pore volume fraction f A = 11% in model (A), the velocity field is localized at a few bright spots, which appears in yellow in the 2D section in Tab. 1, column 1. As expected, the velocity field is small even in regions where the largest pores appear. These zones are, presumably, either disconnected from the percolating paths, or not well connected to them. At higher pore volume fraction f A = 0.22%, the density of hot spots increases (parallel component u 1 , column 2). Furthermore, the velocity field in these spots is higher. In general, the maximum value of the parallel component u 1 increases with the pore volume fraction from 0.06 (f A = 0.22%) to 3.08 (f A = 11%). The minimal value, which is always negative, also increases with the pore volume fraction, i.e. streamlines going backward are increasingly rare events. At the highest pore volume fraction f A = 0.90, the velocity component u 1 is very low in regions comprising 2 or 3 obstacles close to one another (shown in dark blue). It is also significantly higher in a few large zones of the pores (shown in red). This suggests that, even at high pore volume fraction, only a fraction of the pores participates to the overall permeability. At the pore volume fraction f A = 0.90, the transverse components u 2 and u 3 are much more homogeneous than u 1 .

The field patterns for the parallel component u 1 observed on the FFT maps is confirmed by the field histogram Fig. 6 (a). When f = 90%, the field histogram P u1 is strongly non-symmetric. A local maximum field histogram is observed at u 1 ≡ 0, corresponding to regions located around clusters of obstacles. A second local maximum is observed at higher values of the velocity field, u 1 ≈ 0.6. At lower pore volume fraction, the field histogram is non-symmetric. However, they contain one local maximum only at u 1 = 0. The histogram for the transverse component u 2 , or equivalently, that of u 3 is represented in Fig. 6 (b). As for the parallel component, we observe a sharp peak around u 2 = 0 for the density of probability P u2 . Except for the highest pore volume fraction (f = 0.95) on Fig. 6 (b), the transverse component shows a strongly non-Gaussian distribution of the velocity.

Streamlines and orientation of the velocity field

We now extract streamlines from the 3D velocity fields, defined as paths parallel to the velocity vector at each point. They span the microstructure domain from side (of equation x 1 = -L/2) to the opposite one (of equation x 1 = L/2). The streamlines's tortuosity τ is represented in Fig. (7a) for models (A) and (B) at increasing pore volume fraction and compared to the geometrical tortuosity τ A,B Geo . Each point corresponds to one random configuration and different symbols are used for different system sizes. Dotted straight lines indicate the percolation threshold where τ = ∞. For model (A), in the region close to the percolation threshold, around f A < 30%, the system sizes L = 256 and L = 512 used in this work are too small and the numerical results show strong variations. In any case, we expect the tortuosity to be a decreasing function of f . As shown in Fig. (7a) the streamlines's tortuosity for model (A) is only slightly lower than that of model (B), except when the pore volume fraction approaches the percolation threshold. As expected, lower tortuosity implies higher permeability. Indeed, the permeability of model (A) is higher than that of model (B), at given pore volume fraction f (Fig. 2). Nevertheless the permeability in model (A) is often 2 or 3 times higher than that of (B), whereas 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 f the difference in terms of tortuosity is much smaller. For a given pore volume fraction, the morphological tortuosity, deduced from geodesic paths is much smaller than the tortuosity deduced from fluid flow characteristics. As a consequence the overall permeability κ cannot be predicted directly from the geometrical tortuosity τ Geo .

1 1.2 1.5 1.8 2 2.2 2.5 τ τ A Geo τ B Geo τ A (L=256) τ A (L=512) τ B (L=256) τ B (
The static viscuous tortuosity τ A,B 0 and that defined in Eq. ( 17) τ A,B Matyka are represented in Fig. (7b). They are computed as a function of the porosity f for both models A and B. The tortuosity τ A,B 0 is much higher than the streamlines tortuosity τ A,B , especially in the dilute regime f ≈ 1.

The tortuous behavior of the velocity field at low pore volume fraction is confirmed by the covariance of orientation C α (h) plotted in Fig. 8 for model (A). The latter decreases with h at low pore volume fraction, and stays nearly constant for h ≥ h c , where h c depends on f . The length h c represents the typical size of regions where the velocity field is locally correlated due to viscosity. Beyond h c , the finite covariance C α (h) is a result of macroscopic coupling, i.e. the overall direction of the fluid follows that of the applied pressure gradient. As expected, the limiting value of C α (h) increases with f , because of the tortuous path followed by the fluid near the percolation threshold, and h c increases with f , because the size of the channels followed by the fluid also increases with f . In the limit f = 1, we expect C α (h) ≡ 1, a result consistent with Poiseuille flow.

Conclusion

In this numerical work, we used FFT-based full-fields computations for the Stokes flow to estimate the effective permeability of two simple 3D Boolean random models, where the fluid flows inside spheres or its complementary set.

f = 92% f = 88% f = 84% f = 80% f = 68% f = 56% f = 43% f = 32% f = 24% f = 15% f = 11% f = 96%
Fig. 8 Covariance of the orientation of the velocity field Cα(h) for model (A), for various increasing pore volume fraction f .

Very large system sizes were required to carry out computations, especially near the percolation points. For volume fractions 5% higher than the percolation threshold and a relative precision of 5% on the permeability, the size of the RVE becomes larger than 80 times the spheres's radius. No analytical formula was found to provide meaningful estimates of the effective permeability, for any considered model. Similar tortuosity ratios were observed for the fluid's streamlines, for both models, except in the vicinity of the percolation threshold. This is contrast with the geometrical tortuosity, which is significantly lower than that of the streamlines.

The effective permeability was found to be quite different for the two models, with the highest one recovered for spherical obstacles and corresponding to the lowest percolation threshold. Furthermore, as highlighted by the field maps and histograms, at moderate pore volume fraction, the velocity field is heterogeneous, and very small in regions located around or in-between clusters of obstacles. The velocity histograms were found to be non-Gaussian in general. These results confirm that the permeability, a highly structure-dependent property, is not only affected by the tortuosity but also by the distribution of the channels's size and surface areas.

Fig. 1

 1 Fig. 1 Boolean model of spheres discretized on a grid of 512 3 voxels. Spheres diameter: 15 voxels; Spheres volume fraction: 10%.

  for model (A) f must be replaced by 1 -f in relation (27).The point-variance D 2 κ is numericaly obtained by computing second moments of the velocity field over the voxels grid. The variances D 2 κ (V ) are computed by first subdividing a large volume V into a set of non-overlapping subdomains, each of size V . In practice, cuboidal domains are considered and

Fig. 2

 2 Fig. 2 Effective permeability κ A,B of the Boolean models (A) and (B) (disks and circles) with respect to the pore volume fraction f , at two resolutions. Solid lines: Carman-Kozeny estimate κ CK and upper-bound κ U B .

  Fig. (5a) in blue (model A) and red (model B). The solid line in black represents the upper bound from Doi (

Fig. 4

 4 Fig. 4 Symbols: Normalized RVE lengths A,B κ

Fig. 5

 5 Fig. 5 Effective permeability κ A,B of models A (black disks) and B (white disks) vs. 1 -f , when f → 1 (a) and f -f A,B c when f f A,B c is close to the percolation threshold (b). The FFT results were computed on grids of 256 3 (small disks) and 512 3 voxels (large disks). Solid lines: upper bound κ U B (black). (a): analytical results (28) for models A (red) and B (blue) and (b): powerlaw fits (29) of the FFT data for models A (blue) and B (red).

  A c for model (A) and f → f B c for model (B) is represented in Fig. (5b). In the latter, the percolation thresholds are fixed to f → f A c = 3.17%, f B c = 28.95% and the permeability is plotted as a function of f -f A,B c in log-log scale. We consider volumes discretized on grids of 256 3 and 512 3 voxels and fit the data with power laws in a region of interest |f -f A,B c | < 0.2 close to the percolation threshold. The latter are shown in blue (256 3 voxels) and red (512 3 voxels).

  maps of the velocity field u in model (A): 2D sections of 512 × 512 voxels, in the (e 1 , e 2 ) plane with e 1 oriented from left to right and e 2 top to bottom. First row: norm ||u||, rows 2 -4: components u 1 , u 2 and u 3 resp. at increasing pore volume fraction f = 11%, 22%, 52%, 90% (columns 1 -4 resp.). A macroscopic pressure gradient δP = 1 is applied in the direction e 1 . Values are given by the corresponding color scales. Continuation of Tab. 1: FFT maps of the velocity field u in model (B).

Fig. 6

 6 Fig. 6 (a) Histogram of the parallel component of the velocity field u 1 in model B, for three pore volume fraction, in lin-log plot. (b) Histogram of the transverse component of the velocity vector u 2 in model B for various pore volume fraction, in lin-log plot.

Fig. 7

 7 Fig. 7 (a) Tortuosity τ A,B of the fluid streamlines computed using FFT data for models A (filled disks) and B (circles). Diamonds: geometrical tortuosity τ A,B Geo . (b) Viscuous tortuosity τ A,B 0 vs. porosity f (triangles) for models A (blue) and B (red). Squares: tortuosity τ A,B Matyka

Mines ParisTech, PSL Research University, Center for Mathematical Morphology, 35 rue S t -Honoré, 77300 Fontainebleau, France. Corresponding author: F. Willot, E-mail: francois.willot@mines-paristech.fr, Tel.: +33 1 64 69 47 06, Fax.: +33 1 64 69 47 07

Acknowledgements The research leading to the results presented has received funding from the European Union's Seventh Framework Programme (FP7 / 2007-2013) for the Fuel Cells and Hydrogen Joint Technology Initiative under grant agreement 303429. The authors thank an anonymous referee for useful remarks.