# 3D Modelling of Mass Transfer into Bio-Composite

Marouane Kabbej, Valérie V. Guillard, Helene Angellier-Coussy, Caroline Wolf, Nathalie Gontard, Sebastien Gaucel

► To cite this version:

Marouane Kabbej, Valérie V. Guillard, Helene Angellier-Coussy, Caroline Wolf, Nathalie Gontard, et al. 3D Modelling of Mass Transfer into Bio-Composite. Polymers, 2021, 13 (14), pp.2257. 10.3390/polym13142257. [https://hal.inrae.fr/hal-03319281](https://hal.inrae.fr/hal-03319281)

HAL Id: hal-03319281

https://hal.inrae.fr/hal-03319281

Submitted on 12 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

3D Modelling of Mass Transfer into Bio-Composite

Marouane Kabbej ID, Valérie Guillard *ID, Hélène Angellier-Coussy ID, Caroline Wolf, Nathalie Gontard and Sébastien Gaucel ID

IATE, University Montpellier, CIRAD, INRAE, Institut Agro, 34060 Montpellier, France;
marouane.kabbej@umontpellier.fr (M.K.); helene.coussy@umontpellier.fr (H.A.-C.);
wolf.caroline1@gmail.com (C.W.); nathalie.gontard@inrae.fr (N.G.); sebastien.gaucel@inrae.fr (S.G.)

* Correspondence: valerie.guillard@umontpellier.fr

**Abstract:** A three-dimensional model structure that allows considering interphase layer around permeable inclusions is developed to predict water vapor permeability in composite materials made of a matrix Poly(3-HydroxyButyrate-co-3-HydroxyValerate) (PHBV) including Wheat Straw Fiber (WSF) particles. About 500 two-phase structures corresponding to composites of different particles volume fractions (5.14 – 11.4 – 19.52 %v/v) generated using experimental particles’ size distribution have permitted to capture all the variability of the experimental material. These structures have served as a basis to create three-phase structures including interphase zone of altered polymer property surrounding each particle. Finite Element Method (FEM) applied on these structures has permitted to calculate the relative permeability (ratio between composite and neat matrix permeability P/Pm). The numerical results of the two-phase model are consistent with the experimental data for volume fraction lower than 11.4 %v/v but the large upturn of the experimental relative permeability for highest volume fraction is not well represented by the two-phase model. Among hypothesis made to explain model’s deviation, the presence of an interphase with its own transfer properties is numerically tested: numerical exploration made with the three-phase model proves that an interphase of 5 µm thick, with diffusivity of Di ≥ 1 × 10⁻¹⁰ m²·s⁻¹, would explain the large upturn of permeability at high volume fraction.

**Keywords:** 3D numerical modelling; three-phase model; interphase; Finite Element Method; water vapor permeability; composite

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Kabbej, M.; Guillard, V.; Angellier-Coussy, H.; Wolf, C.; Gontard, N.; Gaucel, S. 3D Modelling of Mass Transfer into Bio-Composite. *Polymers* **2021**, *13*, 2257. https://doi.org/10.3390/polym13142257

Academic Editor: Dagmar R. D'hooge

Received: 9 June 2021

Accepted: 3 July 2021

Published: 9 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

1. Introduction

In the last 10 years, a serious interest has been devoted to the modelling of mass transfer properties of composite materials and more particularly of gas and vapor permeability. Permeability is a key functional property that determine the functionality of a given material and thus its application (membranes, packaging, etc.) [1,2]. Predictive modelling of such properties is thus of primary importance to design new tailored materials corresponding to the application requirements. Among innovative materials, composites obtained with the addition of at least two non-miscible constituents with different properties, offer new promising possibilities in terms of applications because of the synergy between components that creates unique macroscopic properties that would not be achievable otherwise from the individual constituents [3,4]. Incorporation of particles in a continuous phase as polymer matrix permits to significantly modulate the mass transfer properties of the resulting composite and makes this strategy very attractive for different fields of applications (barrier materials, membrane gas separation, etc.).

A lot of studies were carried out on the modeling of mass transfer properties in (nano/micro) composites containing either impermeable particles for the application of barrier materials as reviewed, for example, by Cui et al. (2015) or Zid et al. (2018) [5,6] or permeable particles in the field of separation membranes [7]. In the following of this state
of the art, we will focus principally on references dealing with mass transfer modelling in
system where permeable particles are dispersed in a continuous polymer matrix.

Historically, first models developed were analytical relationship for estimation of the overall macroscopic diffusivity or permeability of a composite from structural parameters (volume fraction, aspect ratio, etc.) and mass transfer properties of the individual component [8–10]. Based on analogies with electrics and thermo-mechanic theories, models based upon Maxwell, Bruggeman, Pal or Lewis-Nielsen were proposed [10,11]. Developed for regular dispersion of homogeneous-size particles of idealized geometry (e.g., spheres), they were applied with more or less success to predict permeability of binary systems, mostly in the field of mixed-matrix membranes [7,12]. Their main drawbacks were lack of representativity of the real 3D structure of the composite, assumption of spherical shape of filler which renders the particle geometry insignificant [13], a validity limit restricted to low filler volume fraction (diluted system) [7] and necessity to have access to permeability of the disperse particles which is not always easy to accurately determine [14]. In order to improve their reliability, efforts have been made to complexify the analytical models to consider, for instance, more complex particles geometry such as oblong prolate particles instead of spheres [13,15] or to extend range of validity to high particles volume fraction [16].

In parallel, to overcome the limitations of analytical approaches, some authors have proposed numerical models that attempt to consider the real 3D structure of the composite system and that solve the mass transfer using generally, Finite Element Method (FEM). We can cite, as recent examples in that field, the work of Monsalve-Bravo et al. [17] who model the three-dimensional (3-D) transport problem in full-scale mixed-matrix membranes using FEM or the work of Sharifzadeh and co-authors (2019) who simulated gas diffusion behavior in 3D composite filled with permeable spherical particles randomly dispersed [18]. The two aforementioned studies confirmed that the permeability is positively correlated to filler volume fraction and particle size. However, they did not consider the heterogeneity of the particle size, which is often the case in real bio-composites.

At the opposite of simplifications made in sphere-based models, some authors tried to consider directly the true internal 3D microstructure as Jiang and co-authors (2020) who simulated with a 3D FEM modelling approach, water diffusion in jute/PLA composite using a real structure observed by X-ray tomography and mass transfer properties of the individual components [19]. The numerical results were in good agreement with experimental measurements. Nevertheless, since the simulations were only performed on a unique and small volume of material observed by X-ray tomography, the representativeness of the macroscopic diffusion in the material in its wholeness were questionable.

In the previous cited works whether analytical or numerical ones, the presence of an interphase at the particle/matrix interface which exhibit mass transfer properties differing substantially from those of the bulk matrix is always neglected. Role of this interphase is of paramount importance which can lead to a percolating interphase network inside of the composite, as highlighted by Qiao and Brinson (2009) [20] or Petsi and Burganos (2012) [21] in mixed-matrix membranes and by Zid et al. [22] on nanocomposites with impermeable fillers. Aforementioned studies of Qiao and Brinson and Petsi and Burganos proposed a two-dimensional numerical model to study the impact of interphase zones on the overall properties of the composites (mechanical properties for Qiao and Brinson and effective permeability for Petsi and Burganos). In their study, Zid et al. proposed 3D finite element model to predict mass diffusion in (nano) composites but in their case, particles, with ideal ordered distribution, did not participate to the overall diffusion. Some attempts have also been made to upgrade analytical approaches considering three-phase media—see among others the work of Petropoulos et al. (2015) [23]—with a double binary formula application that consider first particles with interphase (surrounding zones) as pseudo-particles of effective permeability P<sub>E</sub> dispersed in the bulk matrix. P<sub>E</sub> is identified using a standard analytical formula for binary medium. Same standard analytical formula is then applied to the virtual binary composite containing the pseudo particles. To eliminate
limitations imposed by the sphere equations that physically limits the maximal volume
fraction investigated (random packing of congruent spheres imposed a non-negligible
lattice volume), change in particle shape from spherical to cubic was proposed and found
not significant, at least in the polymer-gas permeability area [16]. However, this approach is
restricted to regular dispersion of homogenous particle-size and could not take into account
heterogeneity in interphase thickness, particles size and distribution into the matrix.

The objective of the present paper is to propose a three-dimension, three-phase numer-
ical model to compute effective permeability of bio-composite structures where particles
are permeable and widely contributed to the overall mass transport. This model aims
at deciphering the role of each phase contributing to the overall transport (matrix, inclu-
sion) and, in particular, the role of the third phase (interphase). To go beyond state of
the art, and compensate weaknesses of previous numerical and analytical approaches,
numerical computations of the effective permeability are carried out on model structures
stemming from experimental observation of the composite—poly(3-hydroxybutyrate-co-3-
hydroxyvalerate) (PHBV) as continuous phase and wheat straw fibers (WSF) as dispersed
phase—where particles, idealized as spheroids, are randomly distributed in a representa-
tive volume element (RVE). This approach allows to represent all the complexity of the
real bio-composite while simplifying it to keep sufficiently low computation time which
makes possible exploration of a lot of structures and filler volume fractions. Existence of an
interphase layer around the inclusions is considered: interphase layer was built from the
two-phase structure by allowing this layer to freely overlap with interphase layers around
neighboring particles. The effect of the diffusivity and thickness of the interphase layer on
the effective permeability calculated is discussed, as well as relevance of diffusivity values
chosen for individual components (matrix and particle). Conclusions of the present work is
of paramount importance in the field of material science where modelling of permeability
in composite material considering at the same time 3D structure, permeable particles and
presence of an interphase and percolating interphase network for high filler fraction was
never considered before.

## 2. Materials and Methods

2.1. *Experimental Parameters for Water Vapor Transfer*

To instantiate and validate the numerical models, a real case study of water vapor
transfer into a composite consisting of Poly(3-HydroxyButyrate-co-3-HydroxyValerate)
(PHBV) matrix and wheat straw fiber (WSF) particles was considered, based on previous
experimental investigations performed in our laboratory. Some experimental outputs
coming from this previous work [24] were used as input parameters for the numerical
model. These inputs are explained below and recalled in Table 1.

* Effective moisture diffusivity D [m²·s⁻¹] value of each phase matrix and particle.

* The boundary concentrations of water vapor [mol·m⁻³] in the PHBV matrix in contact with dry air and with humid air (relative humidity of 95%). These concentrations were determined using the experimental water vapor sorption isotherm of PHBV film at 20 °C.

* Water vapor partition coefficient K = 54.52, calculated as the slope of the linear relation between the water concentration in PHBV matrix and WSF particles, obtained from experimental water vapor sorption isotherm at 20 °C for matrix and particles.
<table><caption>**Table 1.** Summary of the experimental parameters for water vapor transfer used in this work.</caption><thead><tr><th>Sample</th><th>Diffusivity <sup>a</sup> × 10<sup>-12</sup><br>[m<sup>2</sup>·s<sup>-1</sup>]</th><th>Permeability × 10<sup>-13</sup><br>[mol·m·m<sup>-2</sup>·s<sup>-1</sup>·Pa<sup>-1</sup>]</th><th>Upper Boundary<br>Concentration <sup>a</sup><br>[mol·m<sup>-3</sup>]</th><th>Lower Boundary<br>Concentration <sup>a</sup><br>[mol·m<sup>-3</sup>]</th><th>Partition Coefficient <sup>a</sup></th></tr></thead><tbody><tr><td>PHBV matrix</td><td>2.615 ± 0.56</td><td>8.29 ± 3.96 <sup>b</sup></td><td>337.14</td><td>0</td><td>54.52</td></tr><tr><td>WSF particle</td><td>18.39 ± 4.93</td><td>1664 ± 451 <sup>c</sup></td><td>-</td><td>-</td><td></td></tr></tbody></table>

## 2.2. 2D Image Analysis

Image processing and analyzing was performed on 2D mosaic images of 2869 Wheat Straw Fibre (WSF) particles by using MATLAB (images provided by Wolf and Berthet from co-authors' laboratory [27,28]). Real particles were assimilated to ellipses whose shape descriptors (major *a*<sub>f</sub> and minor *b*<sub>f</sub> axis) were measured.

## 2.3. 3D Structure Generation

The MATLAB code developed for generating our 3D structures was based on Tschopp MATLAB code [29]. The code generates 3D microstructures composed of a population of non-overlapping ellipsoid particles heterogeneously distributed in size and orientation, within a periodic RVE.

## 2.4. Mathematical Modelling and Geometry

### 2.4.1. 3D Structure Generation

Two-phase system. The two-phase composite structure is generated in a cuboid shaped representative volume element (RVE) defined by $(x, y, z) \in [0, L_x] \times [0, L_y] \times [0, L_z]$, where $z$ is the overall diffusion direction. The RVE is supposed periodic along its vertical faces, to represent an infinite repetitive structure along $x$ and $y$ axis. In the present work, particles are considered as elongated spheroids, i.e., ellipsoids of revolution over the first (major) axis $a_p$, where the second axis $b_p$ and the third axis $c_p$ are equal and lower than the major axis (Figure 1).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Modelling a particle as a spheroid in 3D space. (a) $a_p$: major axis, $b_p$: minor axis, $c_p$: third axis (b) $\theta_{p1}$: Azimuth angle. $\theta_{p2}$: elevation angle.</figcaption></figure>

Structure generation required the RVE size and the particles volume fraction $\varphi_p$ as inputs. First, geometric parameters of particles (major axis $a_p$ and aspect ratio $\alpha_p = a_p/b_p$) are randomly (non-uniform distribution) generated until the target volume fraction of particles is reached, using the stop criterion $|\varphi_{p,\text{target}} - \varphi_p| \le 0.01\%v/v$. Second, these particles are sequentially positioned in the periodic RVE, by decreasing first axis following the steps below.

* Step 1. The position (center coordinates $x_p, y_p, z_p$) and orientation (azimuth $\theta_{p1}$ and elevation $\theta_{p2}$ angle) are randomly drawn using uniform distributions.

<footnote><sup>a</sup> Obtained from dynamic sorption experiments (DVS, Dynamical Vapor Sorption system, Surface Measurement System, London, UK) [24];</footnote>

<footnote><sup>b</sup> Correspond to the average of three experimental sets of measures made in the same laboratory and directly measured from gravimetric experiment (Modified ASTM procedure) [24–26]. <sup>c</sup> Calculated as the product of experimental particle diffusivity by experimental particle solubility obtained from dynamic sorption experiments (DVS, Dynamical Vapor Sorption system, Surface Measurement System, London, UK) [24].</footnote>
* Step 2. The non-overlapping of the particle with the horizontal faces of the RVE ($z = 0$ and $z = L_z$) and with the existing particles is tested.

If the non-overlapping tests were successful, the particle was added to the structure and then the next particle is considered for tests 1–2. If at least one non-overlapping condition was not satisfied, then a new position was drawn, orientation being unchanged, and step 2 was performed again on the updated particle. This last sequence was repeated until the particle was added to the structure.

It should be noted that if a particle intercepted one of the vertical faces of the RVE then, the particle section outside from the RVE was shifted to the opposite face, in order to ensure the periodicity of the RVE (Figure 2). Organizational chart summarizing the structure generation algorithm could be found in Supplementary Materials (Figure S1).

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Representation of the RVE showing: (a) the boundaries conditions at the six composite faces and at the interface matrix–particle for a two–phase structure; (b) the boundaries conditions at the matrix–interphase and interphase–particle interfaces for a three–phase structure.</figcaption></figure>

Three–phase system. The three-phase composite structure was built by adding an interphase of controlled but fixed thickness around each particle of the already generated two-phase structures. The interphase volume was thus different for each particle keeping constant the initial particle volume. Three different interphase thicknesses $e_i$ were considered $e_i = 1 - 2.5 - 5 \ \mu m$.
### 2.4.2. Governing Equations

Mass transfer in the micro-composite is described by Fick's second law of diffusion in both the matrix $\Omega_m$ and the particle $\Omega_p$ domains for the two-phase system (Figure 2a) and also in the interphase domain $\Omega_i$ for the three-phase system (Figure 2b). In stationary regime and in the absence of mass source, this law is expressed by the following partial differential:

$$ \mathrm{div}(\vec{J}_k) = 0 \qquad (1) $$

where

$$ \vec{J}_k(x, y, z) = -D_k \nabla c_k(x, y, z) \quad \text{in} \quad \vec{J}_k(x, y, z) = -D_k \nabla c_k(x, y, z) \qquad (2) $$

where $D_m, D_p$ and $D_i$ [$\mathrm{m}^2 \cdot \mathrm{s}^{-1}$] are the diffusivity coefficient in the matrix, the particle and the interphase domains respectively and they are considered constant (not concentration/temperature/time dependent). $k$ stands for the domain considered, matrix (m), particle (p) or interphase (i). $\vec{J}_m, \vec{J}_p$ and $\vec{J}_i$ are thus the molar surface flux vector [$\mathrm{mol} \cdot \mathrm{m}^{-2} \cdot \mathrm{s}^{-1}$] depending on $c_m, c_p$ and $c_i$ [$\mathrm{mol} \cdot \mathrm{m}^{-3}$] that are the concentration of the water vapor in the matrix, the particle and the interphase domain respectively.

### 2.4.3. Boundary Conditions

Periodic boundary conditions were imposed on the vertical side boundaries of the RVE, which consist to impose equality of concentration and flux on the so-called (source) and (destination) boundaries (Figure 2). The periodic boundary conditions allowed simulating an infinite repetitive structure. Constant concentrations (cte) were imposed on the upper face $\partial\Omega^{m,upper}$ and the lower face $\partial\Omega^{m,lower}$ of the matrix (Equations (3) and (4) and Figure 2a). Particle and interphase domains are assumed to not overlap the upper and lower faces of the RVE.

$$ c_m = cte \neq 0 \quad \text{at} \quad \partial\Omega^{m,upper} \qquad (3) $$

$$ c_m = 0 \quad \text{at} \quad \partial\Omega^{m,lower} \qquad (4) $$

Two-phase system. There are discontinuities in the concentration profile $c_m \neq c_p$ at the matrix-particle interface $\partial\Omega^{m,p}$. In that respect, matrix and particle concentrations are considered linearly dependent at the matrix-particle interface $\partial\Omega^{m,p}$ by the dimensionless partition coefficient $K = c_p/c_m$. To get continuous flux at the matrix-particle interface $\partial\Omega^{m,p}$, a special type of boundary condition using the stiff-spring method [30] was applied:

$$ \left(-\vec{J} \cdot \vec{n}\right)_m = M(c_p - K \cdot c_m) \quad \text{at} \quad \partial\Omega^{m/p} \qquad (5) $$

$$ \left(-\vec{J} \cdot \vec{n}\right)_p = M(K \cdot c_m - c_p) \quad \text{at} \quad \partial\Omega^{p/m} \qquad (6) $$

where $M$ is a (non-physical) velocity [$\mathrm{m} \cdot \mathrm{s}^{-1}$] large enough to let the concentration differences in the brackets approach zero, thereby satisfying $K = c_p/c_m$. This boundary condition gives a continuous flux across the interfaces provided that $M$ is sufficiently large. In all simulations $M$ was taken equal to $1000 \mathrm{m} \cdot \mathrm{s}^{-1}$.

Three-phase system. The interphase was considered to have the same sorption properties than the matrix, continuity of concentration at the matrix-interphase interface $\partial\Omega^{m,i}$ was applied:

$$ c_m = c_i \quad \text{at} \quad \partial\Omega^{m,i} \qquad (7) $$
At the interphase-particle interface $\partial\Omega^{i,p}$, the same condition than previously applied for the two-phase system was adopted with partition coefficient $K = c_p/c_i$. The boundary conditions using the stiff-spring method in the three-phase system are thus:

$$ \left( -\vec{J} \cdot \vec{n} \right)_i = M(c_p - K \cdot c_i) \quad \text{at} \quad \partial\Omega^{i/p} \qquad (8) $$

$$ \left( -\vec{J} \cdot \vec{n} \right)_p = M(K \cdot c_i - c_p) \quad \text{at} \quad \partial\Omega^{p/i} \qquad (9) $$

where $M$ is the same than above.

### 2.4.4. Effective Permeability Evaluation

The solution of the boundary value problem yielded the molar concentration field $c_i(x, y, z)$ [mol·m$^{-3}$·s$^{-1}$] and the molar surface flux vector $\vec{J}(x, y, z)$ [mol·m$^{-2}$·s$^{-1}$] of the permeating species at the discretization points of each domain $\Omega_k$ ($k = m, p, i$). The molar flux $\phi_z$ [mol·s$^{-1}$] (along z-axis) across the upper ($z = L_z$), middle ($z = L_z/2$) and lower ($z = 0$) cross sections was calculated by:

$$ \phi_z(z) = \int_0^{L_x} \int_0^{L_y} J_z(x, y, z) dx dy \qquad (10) $$

where $J_z$ is the molar surface flux (along z-axis) at the discretization points. Then, the effective permeability [mol·m·m$^{-2}$·s$^{-1}$·Pa$^{-1}$] of the micro-composite was finally given by:

$$ P = \frac{\phi_z}{L_x \times L_y} \cdot \left( \frac{L_z - 0}{p_{upper} - p_{lower}} \right) \qquad (11) $$

where $L_x \times L_y$ is the surface of the faces [m$^2$], $L_z$ is the RVE thickness [m], $p_{upper}$ and $p_{lower}$ are the water vapor pressure [Pa] imposed on the upper ($z = L_z$) and the lower ($z = 0$) faces of the RVE respectively and $\phi_z$ is the molar flux [mol·s$^{-1}$].

### 2.5. Numerical Simulations

The 3D boundary value problem of diffusion was solved by using the numerical Finite Element Method (FEM) using COMSOL Multiphysics 5.5 software. An unstructured mesh consisting of tetrahedral elements was used for the discretization of the composite geometry.

For simulation, the Transport of Diluted Species physics interface of the Chemical Reaction Engineering module and the COMSOL CAD Import module, were used. The simulations were performed in DELL computer with Intel Xeon E-2176M Processor (2.7 GHz) and 32 Gb of Ram.

The entire computational procedure (from structure generation to simulation) was driven within the MATLAB environment and was partially automated via the COMSOL LiveLink for MATLAB module (from the step of importing geometry data into COMSOL).

## 3. Results and Discussion

The 3D modelling approach presented above was applied to the prediction of the water vapor permeability into bio-composite materials made of PHBV as continuous phase and WSF as dispersed and permeable phase. Experimental data corresponded to the average of three experimental sets of measures made in the same laboratory for three increasing volume fractions ($\phi_p = 5.14 - 11.4 - 19.52 \%v/v$) [24–26]. For the sake of clarity, in the following, the results are presented and discussed in terms of relative permeability, i.e., the ratio between the composite permeability and that of the neat matrix ($P/P_m$). Experimental results showed an increase of the relative permeability with increasing fiber particles volume fraction due to the hydrophilic nature of particles which are much more permeable than the PHBV matrix considered [24].
Preliminary calculations made using the analytical solution the most representative of the system studied here, i.e., the Maxwell–Wagner–Sillar equation [13] for dilute dispersion of ellipsoids, confirmed that the system was much more complex than a simple binary composite. Indeed, the sharp upturn of the experimental curve for the highest volume fraction (19.52 %v/v) was not captured by the analytical approaches (see Figure S2 in Supplementary Materials). Among hypotheses usually made explaining such a discrepancy, presence of an interphase is the most frequently quoted [20,21]. A three-phase numerical model was thus developed to explore the role of this interphase on the effective macroscopic permeability of the composite. To meet this objective, the strategy applied was first to build the two-phase model and then to add an interphase layer around the inclusions of the two-phase structures to obtain three-phase systems. In order to generate simplified but representative 3D structures, the first step was to determine the particle morphology, size and shape distributions of particles in the experimental composite.

### 3.1. From Particle Morphology to 3D Structure Generation

2D images containing a total number of 3594 WSF particles were analyzed to obtain geometrical descriptors of the particle morphology. Each particle was described by the parameters of its inertia ellipse, major axis $a_f$ and aspect ratio $a_p$ (ratio of major axis to minor axis $a_p = a_p/b_p$), as it allowed to replace the actual particle by an ellipse while preserving centroid, area and orientation. It is worth noting that the very fine WSF particles ($a_f < 5$ µm) were discarded, reducing the total number of particles from 3594 to 2869. Indeed, these very fine particles represent a negligible volume fraction, about 0.118% (volume calculated considering real fibers as spheroids) and their contribution to mass transfer is negligible. Moreover, for numerical simulations, the absence of very fine particles allows to avoid using extremely fine meshes that increase the computation time. The empirical distributions (number frequency) of major axis and aspect ratio of the final set of particles showed monomodal shape (Figure 3) with mean $E = 13.15$ µm and standard deviation $SD = 15.75$ µm for the major axis, $E = 1.9$ and $SD = 1.17$ for the aspect ratio, respectively.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison between experimental distribution and fitted distributions of (a) major axis and (b) aspect ratio obtained from 2869 particles.</figcaption></figure>

The distributions of these shape descriptors were fitted with three theoretical distribution laws, truncated normal, truncated log-normal and truncated exponential laws (Figure 3). Figure 3 shows that the fitted truncated log-normal distribution is the closest to the real distribution. The identified mean and standard deviation of the truncated log-normal distribution after fitting were $E = 13.4$ µm and $SD = 8.14$ µm for the major axis, $E = 1.97$ and $SD = 0.75$ for the aspect ratio. The expressions of the mean and standard deviation of the truncated log-normal distribution are given in Appendix A. Experimental means of both major axis and aspect ratio are well described by the truncated log-normal distribution while standard deviations are a bit underestimated. This latter fact is explained by the presence of a very few number of long particles which are difficult to capture with continuous distributions.
In order to generate 3D structures, the particles were modelled as elongated spheroids, i.e., ellipsoids of revolution along the first axis, the longer one. Major axis $a_f$ and aspect ratio $a_p$ were randomly generated according to the truncated log-normal distributions. In absence of 3D experimental assessment of the particle morphology, assumption of spheroid shape allowed to mimic elongated particles, contrarily to spheres, while keeping smooth boundary which ease the meshing and the numerical solving, contrarily to cylinders or ribbons. No constraints were added for the orientation of particles. The method used (see Section 2) allowed to generate hundreds of composite structures, in moderate time, with particle volume fraction ranging from 0 to 38 %v/v (equivalent to the mass fraction of 48 %wt). Figure 4 illustrates some structures obtained with volume fractions equal to experimental ones ($\varphi_p = 5.14 - 11.4 - 19.52 \%v/v$).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Examples of 3D composite structures (two-phase system) generated for particles volume fraction corresponding to true composite materials (a) $\varphi_p = 5.14 \%v/v$ (122 particles), (b) $\varphi_p = 11.4 \%v/v$ (214 particles) and (c) $\varphi_p = 19.52 \%v/v$ (611 particles). Examples of 3D composite structures (three-phase system) built by adding interphase to the two-phase structure (b) with different interphase thicknesses (d) $e_i = 1$ µm, (e) $e_i = 2.5$ µm and (f) $e_i = 5$ µm.</figcaption></figure>

This approach based on extracting information from real images to generate size distribution was rarely used in the field of mass transfer study, since in most of numerical studies in this field, the generated size distributions are parameterized arbitrarily as in [31,32].

## 3.2. Simulations of Two-Phase Model

### 3.2.1. Selection of Mesh and RVE Sizes

First, the influence of the mesh size on the numerical accuracy of the solution was studied. Preliminary simulations performed on a hundred structures of RVE size of $100 \times 100 \times 300 \ \mu\text{m}^3$ for each composite ($\varphi_p = 5.14 - 11.4 - 19.52 \%v/v$) revealed that the best convergence of solutions was obtained with the mesh element size of [0.45 – 10.5] µm
(predefined in COMSOL) which was then applied in all the further simulations. Indeed, the relative permeability stopped varying approximatively from the mesh element size of [0.45 – 10.5] µm (see Figures S3–S5 in Supplementary Materials). Generally, the mesh was more refined in the regions that require a higher resolution, such as near the matrix-particle interface (Figure 5a) or the matrix-interphase and the interphase-particle interface (Figure 5b).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Tetrahedral mesh presentation of 3D composite structures corresponding to particles volume fraction of $\varphi_p = 11.4 \%v/v$ (240 particles): (a) two-phase system and (b) three-phase system.</figcaption></figure>

Then, influence of the RVE size was investigated. For the present study, fixed RVE thickness ($L_z = 300$ µm) was imposed which was consistent with the thickness of the material represented. It is indeed important to represent the water flux on the entire thickness of the material to capture all the complexity of the system in the direction of the monodirectional flow (with side effects on the upper and lower faces where particles did not outcrop as experimentally observed); therefore, $L_z = 300$ µm was kept as constant value. Further simulations performed on about 100 structures for each RVE size (100 × 100 × 300; 150 × 150 × 300; 200 × 200 × 300 µm³) and each particle volume fraction ($\varphi_p = 5.14 – 11.4 – 19.52 \%v/v$) revealed that the RVE size did not significantly impact the results of mass flux and permeability calculation. For example, for 19.52 $\%v/v$ the permeability [$mol·m·m⁻²·s⁻¹·Pa⁻¹$] found were $7.87 ± 0.29 × 10⁻¹³$ for a RVE of dimensions 100 × 100 × 300 µm³, $7.80 ± 0.21 × 10⁻¹³$ for 150 × 150 × 300 µm³ and $7.88 ± 0.21 × 10⁻¹³$ for 200 × 200 × 300 µm³. Therefore, a minimal RVE of (100 × 100 × 300 µm³) was chosen for all further simulations (Figure 4), permitting to achieve reasonable computational time.

In total, no fewer than 500 two-phase structures of RVE size of (100 × 100 × 300 µm³) were generated using experimentally observed particle size distribution: 185 structures for $\varphi_p = 5.14 \%v/v$, 199 structures for $\varphi_p = 11.4 \%v/v$ and 158 structures for $\varphi_p = 19.52 \%v/v$.

### 3.2.2. Selection of the Number of Structures to Analyze

For each filler load, the variability of numerical vapor permeability was analyzed, depending on the number of structures considered. As expected, accumulation of structures led to a stabilization of the mean permeability and a decrease of its standard deviation. Regarding the numerical variability obtained on water vapor permeability, it was found that numerical standard deviation became acceptable (i.e., less than 10%) for accumulation of results obtained on at least 10 structures, whatever the filler load. In the following, permeability was computed for the whole set of structures generated in order to completely
ensure stable results and also to capture all the numerical variability and get rid of potential outliers caused by minority structures. By multiplying the number of computed structures, the variability of the material structure could be integrated more easily, leading to a better characterization and prediction of the material properties compared to works relying on actual 3D structures, e.g., costly 3D tomographic pictures [19].

### 3.2.3. Numerical Results of the 2-Phase Model

As expected, the numerical permeability values obtained for the 2-phase system were not able to capture the upturn of permeability ratio for high volume fraction, even if good prediction of the of the experimental relative permeability ratio were obtained for the composites of $\phi_p = 5.14 \%v/v$ and $\phi_p = 11.4 \%v/v$ (Figure 6).

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison between experimental and numerical relative permeability $P/P_m$ (calculated using the 2-phase model).</figcaption></figure>

Three different hypotheses could be made to explain this deviation:

1. the physical properties of the fiber particle and especially its diffusivity value would be modified once embedded into the polymer matrix compare to the one measure on the native component,

2. the diffusivity of the polymer matrix, measured before fiber particles addition, would be modified after fiber particles addition and therefore not well representative of what occurs in the composite material,

3. the presence of an interphase, third compartment with its own physical properties, at the interface matrix/particle would influence the overall permeability into the composite.

All these hypotheses were tested in the following, especially on the composite of $\phi_p = 19.52 \%v/v$ where largest deviations were noted.

### 3.2.4. Modification of Particles Diffusivity Values in the Two-Phase Model

The first hypothesis tested in the 2-phase model was the modification of the diffusivity of the fiber particles $D_p$ once embedded in the polymer matrix. Figure 7, obtained on composite ($\phi_p = 19.52 \%v/v$), showed that the relative permeability increased when $D_p$ increased from $1 \times 10^{-14} \text{ m}^2\cdot\text{s}^{-1}$ to $1 \times 10^{-11} \text{ m}^2\cdot\text{s}^{-1}$ and reached a plateau when $D_p \ge 1 \times 10^{-11} \text{ m}^2\cdot\text{s}^{-1}$. However, this increase was not enough to bring the numerical results closer to the experimental data for the composite of $\phi_p = 19.52 \%v/v$ (Figure 8).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Effect of particle diffusivity $D_p$ on numerical relative permeability $P/P_m$ (calculated using the 2-phase model) into the composite of $\varphi_p = 19.52 \%v/v$ (hypothesis tested on 10 structures).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison between experimental and numerical relative permeability $P/P_m$ (calculated using the 2-phase model) for different particle diffusivity values varying from $D_p = 1 \times 10^{-14} \text{ m}^2 \cdot \text{s}^{-1}$ to $D_p = 1 \times 10^{-9} \text{ m}^2 \cdot \text{s}^{-1}$ (hypothesis tested on 10 structures).</figcaption></figure>

In other words, modification of the diffusivity of the fiber particles $D_p$ in the two-phase model was not a valid hypothesis to explain deviation of the 2-phase model.

### 3.2.5. Modification of Matrix Diffusivity Values in the Two-Phase Model

The second hypothesis tested in the 2-phase model was the modification of diffusivity of the polymer matrix induced by fiber particles addition. As shown in Figure 9, the numerical relative permeability $P/P_m$ sharply increased when $D_m$ increased. The two-phase numerical model is highly sensitive to matrix diffusivity $D_m$. The optimal matrix diffusivity, i.e., the one leading to the best fit of the experimental data, was $D_m = 3.08 \times 10^{-12} \text{ m}^2 \cdot \text{s}^{-1}$ for $\varphi_p = 5.14 \%v/v$, $D_m = 3.92 \times 10^{-12} \text{ m}^2 \cdot \text{s}^{-1}$ for $\varphi_p = 11.4 \%v/v$ and $D_m = 10.07 \times 10^{-12} \text{ m}^2 \cdot \text{s}^{-1}$ for $\varphi_p = 19.52 \%v/v$ (Figure 10). These values of the diffu-
sivity were, as expected, close to the experimental one for $\phi_p = 5.14 \%v/v$ and $11.4 \%v/v$ ($D_m = 2.615 \times 10^{-12} \text{ m}^2\cdot\text{s}^{-1}$), as the 2-phase model relatively well predicted the experimental relative permeability for these two particle loads. On the contrary, the composite of $\phi_p = 19.52 \%v/v$ would require a matrix diffusivity about three times higher than that currently used in the 2-phase model. That means that the experimental relative permeability could be well represented by the two-phase model for $\phi_p = 19.52 \%v/v$ when $D_m$ is multiplied by 3. Relevance of this increase by three orders of magnitude in the polymer matrix is questionable but could be justified by modification of polymer crystallinity when high volume fraction of particles is added in the neat matrix. Impact of crystallinity on diffusivity of small-molecule penetrants in semicrystalline polymer is something well described in the literature [33] even if never confirmed for PHBV polymeric matrix. In the specific case of the PHBV matrix under study in the present work, Berthet et al. (2015) [26] measured by WAXD a decrease of crystallinity from 68.1 to 46.8%, between neat matrix and composite containing 20 $\%w/w$ (e.g., $11.4 \%v/v$) of WSF (same particle and same polymer than in the present work). This decrease of crystallinity rate means an increase of amorphous zones and higher mobility of the polymer in the continuous phase, which is generally ascribed to higher diffusivity values. For example, Trifol et al. (2020) [34] measured a water vapor diffusivity value multiplied by 2 between amorphous PLA and semicrystalline PLA with 35% of crystallinity. Considering the same effect (decrease of 35% of crystallinity rate leads to $D_m$ multiplied by 2) for our PHBV/WSF composite ($\phi_p = 11.4 \%v/v$) would lead to a $D_m$ multiplied by 1.5. Extrapolating this to the $\phi_p = 19.52 \%v/v$ composite, supposing that same drop of crystallinity rate would be observed between $\phi_p = 11.4 \%v/v$ and $\phi_p = 19.52 \%v/v$ than between $\phi_p = 11.4 \%v/v$ and neat matrix, would lead to a $D_m$ value multiplied by 3. The hypothesis of strong modification of diffusivity of the polymer matrix induced by fiber particles addition is thus completely credible. Berthet et al. (2015) [26] also noted that in parallel to the decrease of crystallinity, addition of fibers induced a decrease of polymer molecular weight (219,340 to 169,906 g mol$^{-1}$), which is also in favor of higher polymer mobility (less entanglements) and thus increased of $D_m$ value in composite. A modification of the crystal size was also highlighted by the authors: the size of crystals was increased from 1.03 to 1.18 nm in the presence of WSF. This higher crystal size is also in favor of an increase of diffusivity in the continuous phase because of lower tortuosity in the less tight crystal structure.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Effect of matrix diffusivity $D_m$ on numerical relative permeability $P/P_m$ (calculated using the 2-phase model) into the composite of $\phi_p = 19.52 \%v/v$ (hypothesis tested on 10 structures).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison between experimental and numerical relative permeability *P*/*P*<sub>m</sub> (calculated using the 2-phase model): adjustment of numerical model to experimental data by fitting the matrix diffusivity *D*<sub>m</sub> (fitting performed on 10 structures for each composite).</figcaption></figure>

To sum up, hypothesis of a modification of diffusivity of the polymer matrix induced
by fiber particles addition for high filler fraction would be realistic (for φ<sub>p</sub> = 19.52 %v/v
it is necessary to multiply D<sub>m</sub> by 3) and could be explained by a concomitant decrease of
crystallinity rate, a decrease of polymer molecular weight and an increase of crystal size in
the PHBV/WSF composite.

It is however difficult to decipher the individual role of each effect (crystal size,
molecular weight, crystallinity rate), all the more that increase of crystal size in composite
hinder fiber/matrix adhesion [26] and would lead to interfacial phenomenon such as
creation of an interphase with its own transfer properties at the interface particle/matrix.
It is thus particularly relevant to investigate in the following the role of this interphase.

3.3. *Simulations of Three-Phase Model*

To build the three-phase model, an interphase of controlled but fixed thickness was
added around each particle of the two-phase structures previously generated, reflecting
a fixed degree of lack of adherence between inclusions and polymer. Interphase layers
were allowed freely to overlap with interphase layers around neighboring particles. Three
different thicknesses were explored: eᵢ = 1 – 2.5 – 5 µm, 5 µm being a worst-case scenario
following literature analysis on that topic. Indeed, relevant thickness range of 0–2 µm
was found in similar composite systems [35]. Also, a micro-thermal analysis on a glass
fiber/epoxy composite revealed the presence of a zone of 4 µm thickness around the
glass fibers with higher molecular mobility that is generally ascribed to higher diffusivity
values [36].

A total of 69 structures for (φₚ = 5.14 %v/v), 93 structures for (φₚ = 11.4 %v/v) and
39 structures for (φₚ = 19.52 %v/v) were used for simulations (see Figure 4d–f for some
examples). Simulations of the three-phase model were performed with several modalities
for the diffusivity Dᵢ of the interphase and its thickness eᵢ. It must be noted that the standard
deviation on the calculated relative permeability increases as the thickness and diffusivity of
the interphase increase (Figure 11): this could be ascribed to the fact that less structures were
analyzed for high φₚ and high φᵢ values. Indeed, numerical constraints on structures with
high φₚ and high φᵢ have limited the number of structures analyzed and thus the number of
permeability results considered. The equivalent graphs of Figure 11 for φₚ = 5.14 %v/v and
φₚ = 11.4 %v/v were given in the Supplementary Materials (Figures S6 and S7).
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Evolution of the numerical relative permeability of the composite ($\varphi_p = 19.52 \%v/v$) as a function of the diffusivity of interphase and for different thicknesses of the interphase. The numerical results (symbols) correspond to the average of the relative permeability of 14 structures ($e_i = 1 \ \mu\text{m}$), 19 structures ($e_i = 2.5 \ \mu\text{m}$) and 6 structures ($e_i = 5 \ \mu\text{m}$). Lines represent the experimental relative permeability for $\varphi_p = 19.52 \%v/v$: solid line for mean value, dashed lines for mean ± standard deviation.</figcaption></figure>

Figure 11 showed that the numerical relative permeability is positively correlated to $D_i$ and $e_i$. Compared to the 2-phase results ($e_i = 0$), the presence of an interphase can lead to an increase of the relative permeability, for $D_i > D_m$, or a decrease of the relative permeability, for $D_i < D_m$. For a thin interface, $e_i = 1 \ \mu\text{m}$, the diffusivity of the interface, $D_i$, had a strong effect on the relative permeability for $D_i < D_m$ but very limited impacts for $D_i > D_m$. For larger thicknesses ($e_i = 2.5 - 5 \ \mu\text{m}$), the results exhibit a greater impact of $D_i$ on the relative permeability for both $D_i < D_m$ and $D_i > D_m$.

The increase of $e_i$ led, as expected, to an increase of the volume fraction of the interphase and so to an increase (resp. decrease) of the water vapor flux for $D_i > D_m$ (resp. $D_i < D_m$), e.g., for $\varphi_p = 19.52 \%v/v$, the volume fraction of the interphase was $\varphi_i = 5.58 \pm 1.37 \%v/v$ for $e_i = 1 \ \mu\text{m}$, $\varphi_i = 16.03 \pm 3.2 \%v/v$ for $e_i = 2.5 \ \mu\text{m}$ and $\varphi_i = 29.24 \pm 5.49 \%v/v$ for $e_i = 5 \ \mu\text{m}$. This volume fraction, $\varphi_i$, can become equivalent or even higher than that of the particles, $\varphi_i \ge \varphi_p$. Even if the relationship between volume fraction of the interface and relative permeability is not straightforward (Figures S8–S10 in Supplementary Materials), the general trend observed is globally an increase of relative permeability for $D_i > D_m$ when the interphase volume fraction increased. For high volume fraction of the interphase, the interphases are strongly connected and some continuous pathways appeared along the main diffusion axis (z axis), which can act as percolation pathways and have an additional increasing impact on the mass transfer.

This typical impact of the presence of interphases on mass transfer within composites has already been observed by Zid et al. [22] for nanocomposite with impermeable fillers and by Petsi and Burganos [21] in mixed matrix membranes. Zid et al. [22] noticed that the interphase layers, depending on their diffusivity (weakly or highly diffusive) can be either beneficial or totally detrimental to the nanocomposite overall barrier properties. They particularly highlighted the impact of the effect of continuous diffusion paths, which may occur between overlapping interphases, that are particularly critical for the barrier performance in the case of highly diffusive interphases. This is perfectly in line with what has been observed for high filler volume fraction in the present work. Petsi and Burganos [21] went even further in their conclusions obtained on mixed matrix membranes by investigating the impact of constant interphase thickness (hypothesis made in the present work) or thickness that is proportional to the particle size. They concluded from the numerical exploration of their system that when the interphase thickness remains constant
during particle decrease, the permeability increased thanks to the increased role of the interphase layer. It may be anticipated from this result that role of interphase layer would be higher in the case of composite with smaller particles.

Indeed, a more detailed analysis of the results showed that the number of particles but also the spatial distribution impacted the relative permeability (Figures S8–S10 in Supplementary Materials). Indeed, as the interphase thickness was considered constant for all the particles, the total volume of the interphase is related to the dimensions of the spheroid particles, and to the position of particles as intersection of interphases can occur. Consequently, two different structures having the same particle volume fraction but different number of particles or different particle spatial distributions can lead to significantly different interphase volume fractions or interphase shape, and so different impacts on relative permeabilities, e.g., for $\varphi_p = 19.52 \%v/v$ and $e_i = 2.5 \ \mu\text{m: (i)} n_p = 255$, $\varphi_i = 15.42 \%v/v$ and (ii) $n_p = 253$, $\varphi_i = 13.75 \%v/v$ (Figure S9).

Figure 11 revealed that, for $\varphi_p = 19.52 \%v/v$, the three-phase model with thickness of $e_i = 1 \ \mu\text{m}$ was not sufficient to bring the numerical results closer to the experimental data whatever the $D_i$ value investigated between $1 \times 10^{-14}$ and $1 \times 10^{-6} \ \text{m}^2 \cdot \text{s}^{-1}$. Good fit of the three-phase model on experimental relative permeability for this composite could be achieved only with thicker interphase (either $e_i = 2.5 \ \mu\text{m}$ or $e_i = 5 \ \mu\text{m}$) and by considering interphase diffusivity $D_i \ge 1 \times 10^{-9} \ \text{m}^2 \cdot \text{s}^{-1}$ for 2.5 µm of interphase thickness and $D_i \ge 1 \times 10^{-10} \ \text{m}^2 \cdot \text{s}^{-1}$ for 5 µm of interphase thickness. According to these results, the interphase would be as diffusive as or even more diffusive than the particles $D_i \ge D_p$.

Considering the same interphase thickness and diffusivity for each particle volume fraction would not be the right hypothesis to represent the increase in relative permeability in the composite as a function of the particles volume fraction (Figure 12). According to Figure 12, the two-phase model would be enough for $\varphi_p = 5.14 \ %v/v$ for a good model fitting of experimental relative permeability value. Considering an interphase of 2.5 µm and $D_i = 1 \times 10^{-10} \ \text{m}^2 \cdot \text{s}^{-1}$ would slightly improve prediction for $\varphi_p = 11.4 \ %v/v$ even if the prediction of the two-phase model was enough and included in the experimental variability of results obtained at this volume fraction. On the contrary, it would be absolutely necessary to consider an interphase of 5 µm (with $D_i = 1 \times 10^{-10} \ \text{m}^2 \cdot \text{s}^{-1}$) for the $\varphi_p = 19.52 \ %v/v$. However, reliability of this hypothesis still needs to be experimentally validated.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison between experimental and numerical relative permeability: Effect of the particles volume fraction $\varphi_p$ and effect of the diffusivity $D_i$ and thickness $e_i$ of the interphase.</figcaption></figure>

To sum up, the presence of an interphase layer around the inclusions, globally more diffusive than the particles could explain the high upturn of permeability curve for high par-
ticle volume fraction ($\varphi_p = 19.52 \%v/v$). A minimal interphase thickness would be required to form a continuous diffusion paths, which may occur between overlapping interphases.

## 4. Conclusions

In this paper, a 3D three-phase model was developed in order to predict water vapor permeability in micro-composite (PHBV matrix/WSF particles). This 3D study was essential to improve understanding of structure/material transfer relationships in bio-composite materials containing permeable particles. The developed two-phase model led to good prediction of the experimental relative permeability for the composites of $\varphi_p = 5.14 \%v/v$ and $\varphi_p = 11.4 \%v/v$. Nevertheless, for the composite of $\varphi_p = 19.52 \%v/v$, the results revealed that the mass transfer properties cannot be predicted directly from the mass transfer properties of each phase (matrix and particle) considering a binary system. Among the different hypothesis tested to explain the gap between model prediction and experimental results, an increase of matrix diffusivity ($\times 3$) in the composite with $19.52 \%v/v$ of volume fraction would lead to a good model fitting. Alternatively, presence of an interphase of 5 µm between matrix and particles with its own diffusivity $1 \times 10^{-10} \text{ m}^2 \cdot \text{s}^{-1}$ would permit to perfectly fit the experimental data for the highest volume fraction. Thorough analysis of published literature on PHBV based composites has permitted to conclude that probably both phenomena concomitantly occur. Results of this work highlighted that presence or not of an interphase, and thus choice of the model to be used would be strongly dependent on the particle volume fraction and its size distribution. This plenty justified the choice of developing a 3D, three-phase model based on experimental particles' size distribution to well decipher the determinants of permeability change in composite of high volume fraction. In addition, such approach has permitted to multiply the number of computed structures and thus, to take into consideration the variability of the material structure leading to a better characterization and prediction of the material properties. The 3D FEM model proposed here was necessary to capture all the complexity of the bio-composites studied and to better construe impact of interphase, matrix and particle properties on macroscopic water vapor permeability.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/polym13142257/s1. Figure S1. Numerical process to generate the 3D microstructure of the RVE performed in MATLAB. Figure S2. Comparison between experimental relative permeability $P/P_m$ with calculated ones by using Maxwell-Wagner-Sillar equation. Figure S3. Evolution of the numerical relative permeability of water vapor as a function of the mesh element size for the composite of $\varphi_p = 5.14 \%v/v$. Mesh tests are performed with 10 structures. Figure S4. Evolution of the numerical relative permeability of water vapor as a function of the mesh element size for the composite of $\varphi_p = 11.4 \%v/v$. Mesh tests are performed with 10 structures. Figure S5. Evolution of the numerical relative permeability of water vapor as a function of the mesh element size for the composite of $\varphi_p = 19.52 \%v/v$. Mesh tests are performed with 8 structures. Numerical results corresponding to some mesh element sizes are not available due to mesh errors encountered on some simulations. Figure S6. Evolution of the numerical relative permeability of water vapor as a function of the diffusivity and thickness of the interphase for the composite of $\varphi_p = 5.14 \%v/v$: Comparison between experimental and numerical results. The numerical results (bullets) correspond to the average of the relative permeability of 20 structures ($e_i = 1 \mu\text{m}$), 23 structures ($e_i = 2.5 \mu\text{m}$) and 26 structures ($e_i = 5 \mu\text{m}$). The volume fraction of the interphase was $\varphi_i = 1.85 \pm 0.34 \%v/v$ for $e_i = 1 \mu\text{m}$, $\varphi_i = 5.71 \pm 1.24 \%v/v$ for $e_i = 2.5 \mu\text{m}$ and $\varphi_i = 15.19 \pm 3.36 \%v/v$ for $e_i = 5 \mu\text{m}$. Figure S7. Evolution of the numerical relative permeability of water vapor as a function of the diffusivity and thickness of the interphase for the composite of $\varphi_p = 11.4 \%v/v$. The numerical results (bullets) correspond to the average of the relative permeability of 36 structures ($e_i = 1 \mu\text{m}$), 36 structures ($e_i = 2.5 \mu\text{m}$) and 21 structures ($e_i = 5 \mu\text{m}$). The volume fraction of the interphase was $\varphi_i = 3.32 \pm 0.9 \%v/v$ for $e_i = 1 \mu\text{m}$, $\varphi_i = 9.93 \pm 2.69 \%v/v$ for $e_i = 2.5 \mu\text{m}$, and $\varphi_i = 22.47 \pm 6.65 \%v/v$ for $e_i = 5 \mu\text{m}$. Figure S8. Evolution of the numerical relative permeability of water vapor as a function of the diffusivity of the interphase for the composite of $\varphi_p = 19.52 \%v/v$ and $e_i = 1 \mu\text{m}$: Comparison between experimental and numerical results. The numerical results corresponded to the relative permeability of 14 structures. $\varphi_i$ and $np$ are the interphase volume fraction and the number of particles respectively. Figure S9. Evolution of
the numerical relative permeability of water vapor as a function of the diffusivity of the interphase for the composite of $\phi_p = 19.52 \%v/v$ and $e_i = 2.5 \ \mu m$: Comparison between experimental and numerical results. The numerical results corresponded to the relative permeability of 19 structures. $\phi_i$ and $np$ are the interphase volume fraction and the number of particles respectively. Figure S10. Evolution of the numerical relative permeability of water vapor as a function of the diffusivity of the interphase for the composite of $\phi_p = 19.52 \%v/v$ and $e_i = 5 \ \mu m$: Comparison between experimental and numerical results. The numerical results corresponded to the relative permeability of 6 structures. $\phi_i$ and $np$ are the interphase volume fraction and the number of particles respectively.

**Author Contributions:** Conceptualization, V.G. and N.G.; methodology, M.K. and S.G.; software, M.K. and S.G.; validation, M.K., S.G. and V.G.; formal analysis, M.K.; investigation, M.K. and C.W.; resources, V.G. and N.G.; data curation, M.K. and C.W.; writing—original draft preparation, M.K.; writing—review and editing, V.G., S.G. and H.A.-C.; visualization, M.K. and H.A.-C.; supervision, V.G. and S.G.; project administration, N.G.; funding acquisition, V.G. and N.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 773375.

**Data Availability Statement:** Raw simulation data used in the article could be uploaded from the publicly available dataset: https://doi.org/10.15454/D0GZDQ, accessed on 7 July 2021.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature

Abbreviations

RVE Representative Volume Element

FEM Finite Element Method

PHBV Poly(3-HydroxyButyrate-co-3-HydroxyValerate)

WSF Wheat Straw Fibers

2D, 3D Two and Three Dimension

E Mean

SD Standard Deviation

PDF Probability Density Function

CDF Cumulative Distribution Function

cte Constant value

Latin symbols

$L_x, L_y, L_z$ RVE length along x-axis, y-axis and z-axis (m)

$a_p, b_p, c_p$ Major, minor and third axis of the particle (m)

$x_p, y_p, z_p$ Center coordinates of the particle (m)

$D_k$ Diffusivity of water vapor in the phase k ($m^2 \cdot s^{-1}$)

$P_k$ Permeability of water vapor in the phase k ($mol \cdot m^{-2} \cdot s^{-1} \cdot Pa^{-1}$)

$P$ Permeability of water vapor in the composite ($mol \cdot m^{-2} \cdot s^{-1} \cdot Pa^{-1}$)

$c_k$ Concentration of water vapor in the phase k ($mol \cdot m^{-3}$)

$\vec{J}_k$ Molar surface flux vector of water vapor in the phase k ($mol \cdot m^{-2} \cdot s^{-1}$)

$p_{upper} - p_{lower}$ Water vapor pressure differential across the film (Pa)

$K$ Partition coefficient: concentration ratio between particles and matrix at equilibrium ($c_p/c_m$) (-)

$M$ Velocity (non-physical property) ($m \cdot s^{-1}$)

$e_i$ Interphase thickness (m)

$n_p$ Number of particles (-)
Greek symbols

$\alpha_p$ Aspect ratio: ratio between major and minor axis of the particle ($a_p/b_p$) (-)

$\varphi_k$ Volume fraction of the phase k: ratio between the volume of the phase k and the composite volume ($\%v/v$)

$\phi_z$ Molar flux (along z-axis) of water vapor across a composite face (mol·s⁻¹)

$\theta_{p1}$ Azimuth angle: angle between the x-axis and the orthogonal projection of the semi-major axis onto the xy-plane (degree°)

$\theta_{p2}$ Elevation angle: angle between the semi-major axis and its orthogonal projection onto the xy-plane (degree°)

Subscripts

m Matrix

p Particle

i Interphase

Appendix A

Probability Density Function (PDF) of log-normal distribution $f(x)$ of random variable $x$ with parameters $\mu$ and $\sigma$ was defined as

$$f(x) = \frac{1}{x\sigma\sqrt{2\pi}}e^{-\frac{(\ln(x)-\mu)^2}{2\sigma^2}} \quad x \in ]0, +\infty[ \qquad (A1)$$

The mean $E$, variance $V$ and standard deviation $SD$ of log-normal distribution were:

$$E = e^{\mu + \frac{\sigma^2}{2}}; \quad V = (e^{\sigma^2} - 1)e^{2\mu + \sigma^2}; \quad SD = \sqrt{V} \qquad (A2)$$

PDF of truncated log-normal distribution $g(x)$ of random variable $x$ with parameters $\mu$, $\sigma$, $x_{min}$ and $x_{max}$ was defined as:

$$g(x) = \frac{f(x)}{\Phi\left(\frac{\ln(x_{min}) - \mu}{\sigma}\right) - \Phi\left(\frac{\ln(x_{max}) - \mu}{\sigma}\right)} \quad x \in [x_{min}, x_{max}] \qquad (A3)$$

where $f(x)$ is the PDF of log-normal distribution and $\Phi(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-x^2/2} dx$ is the Cumulative Distribution Function (CDF) of the standard normal distribution. The mean $E_{trunc}$, variance $V_{trunc}$ and standard deviation $SD_{trunc}$ of truncated log-normal distribution were:

$$E_{trunc} = E\frac{T_2}{T_1} \qquad (A4)$$

$$V_{trunc} = E^2\left(\frac{T_3}{T_1} - \frac{T_2^2}{T_1^2}\right) + V\frac{T_3}{T_1} \qquad (A5)$$

$$SD_{trunc} = \sqrt{V_{trunc}} \qquad (A6)$$

with

$$T_1 = \Phi\left(\frac{\ln(x_{max}) - \mu}{\sigma}\right) - \Phi\left(\frac{\ln(x_{min}) - \mu}{\sigma}\right) \qquad (A7)$$

$$T_2 = \Phi\left(\frac{\ln(x_{max}) - \mu}{\sigma} - \sigma\right) - \Phi\left(\frac{\ln(x_{min}) - \mu}{\sigma} - \sigma\right) \qquad (A8)$$

$$T_3 = \Phi\left(\frac{\ln(x_{max}) - \mu}{\sigma} - 2\sigma\right) - \Phi\left(\frac{\ln(x_{min}) - \mu}{\sigma} - 2\sigma\right) \qquad (A9)$$

References

1. Angellier-Coussy, H.; Guillard, V.; Guillaume, C.; Gontard, N. Role of packaging in the smorgasbord of action for sustainable food consumption. *Agro Food Ind. Hi Tech* **2013**, *24*, 15–19.
2. Guillard, V.; Gaucel, S.; Fornaciari, C.; Angellier-Coussy, H.; Buche, P.; Gontard, N. The Next Generation of Sustainable Food Packaging to Preserve Our Environment in a Circular Economy Context. *Front. Nutr.* **2018**, *5*, 1–13. [CrossRef] [PubMed]

3. Berthet, M.-A.; Angellier-Coussy, H.; Chea, V.; Guillard, V.; Gastaldi, E.; Gontard, N. Sustainable food packaging: Valorising wheat straw fibres for tuning PHBV-based composites properties. *Compos. Part A Appl. Sci. Manuf.* **2015**, *72*, 139–147. [CrossRef]

4. Berthet, M.A.; Angellier-Coussy, H.; Guillard, V.; Gontard, N. Vegetal fiber-based biocomposites: Which stakes for food packaging applications? *J. Appl. Polym. Sci.* **2016**, *133*. [CrossRef]

5. Cui, Y.; Kumar, S.; Rao Kona, B.; van Houcke, D. Gas barrier properties of polymer/clay nanocomposites. *RSC Adv.* **2015**, *5*, 63669–63690. [CrossRef]

6. Zid, S.; Zinet, M.; Espuche, E. Modeling Diffusion Mass Transport in Multiphase Polymer Systems for Gas Barrier Applications: A Review. *J. Polym. Sci.* **2018**, *56*, 621–638. [CrossRef]

7. Monsalve-Bravo, G.M. Modeling permeation through mixed-matrix membranes: A review. *Processes* **2018**, *6*, 172. [CrossRef]

8. Fricke, H. The maxwell-wagner dispersion in a suspension of ellipsoids. *J. Phys. Chem.* **1951**, *57*, 934–937. [CrossRef]

9. Barrer, R.M.; Petropoulos, J.H. Diffusion in heterogeneous media: Lattice and parallelepipeds in a continuous phase. *J. Appl. Phys.* **1961**, *12*, 691–697. [CrossRef]

10. Petropoulos, J.H. A comparative study of approaches applied to the permeability of binary composite polymeric materials. *J. Polym. Sci. Polym. Phys. Ed.* **1985**, *23*, 1309–1324. [CrossRef]

11. Maxwell, J.C. A Treatise on Electricity and Magnetism 1873. Available online: https://www.nature.com/articles/007478a0 (accessed on 7 July 2021).

12. Gonzo, E.E.; Parentis, M.L.; Gottifredi, J.C. Estimating models for predicting effective permeability of mixed matrix membranes. *J. Memb. Sci.* **2006**, *277*, 46–54. [CrossRef]

13. Rafiq, S.; Maulud, A.; Man, Z.; Ibrahim, M.; Mutalib, A.; Ahmad, F.; Khan, A.U.; Khan, A.L.; Ghauri, M.; Muhammad, N. Modelling in Mixed Matrix Membranes for Gas Separation. *Can. J. Chem. Eng.* **2015**, *93*, 88–95. [CrossRef]

14. Thoury-Monbrun, V.; Angellier-Coussy, H.; Guillard, V.; Legland, D.; Gaucel, S. Impact of Two-Dimensional Particle Size Distribution on Estimation of Water Vapor Diffusivity in Micrometric Size Cellulose Particles. *Materials* **2018**, *11*, 1712. [CrossRef]

15. Rafiq, S.; Man, Z.; Maulud, A.; Muhammad, N.; Maitra, S. Separation of CO2 from CH4 using polysulfone/polyimide silica nanocomposite membranes. *Sep. Purif. Technol.* **2012**, *90*, 162–172. [CrossRef]

16. Papadokostaki, K.G.; Minelli, M.; Doghieri, F.; Petropoulos, J.H. A fundamental study of the extent of meaningful application of Maxwell's and Wiener's equations to the permeability of binary composite materials. Part II: A useful explicit analytical approach. *Chem. Eng. Sci.* **2015**, *131*, 353–359. [CrossRef]

17. Monsalve-bravo, G.M.; Smart, S.; Bhatia, S.K. Simulation of multicomponent gas transport through mixed-matrix membranes. *J. Memb. Sci.* **2019**, *577*, 219–234. [CrossRef]

18. Sharifzadeh, M.; Zamani Pedram, M.; Ebadi Amooghin, A. A new permeation model in porous filler-based mixed matrix membranes for CO₂ separation. *Greenh. Gases Sci. Technol.* **2019**, *9*, 719–742. [CrossRef]

19. Jiang, N.; Li, Y.; Li, D.; Yu, T.; Li, Y.; Xu, J.; Li, N.; Marrow, T.J. 3D finite element modeling of water diffusion behavior of jute/PLA composite based on X-ray computed tomography. *Compos. Sci. Technol.* **2020**, *199*, 108313. [CrossRef]

20. Qiao, R.; Catherine Brinson, L. Simulation of interphase percolation and gradients in polymer nanocomposites. *Compos. Sci. Technol.* **2009**, *69*, 491–499. [CrossRef]

21. Petsi, A.J.; Burganos, V.N. Interphase layer effects on transport in mixed matrix membranes. *J. Memb. Sci.* **2012**, *421–422*, 247–257. [CrossRef]

22. Zid, S.; Zinet, M.; Espuche, E. 3D numerical analysis of mass diffusion in (nano) composites: The effect of the filler-matrix interphase on barrier properties. *Model. Simul. Mater. Sci. Eng.* **2020**, *28*, 75003. [CrossRef]

23. Petropoulos, J.H.; Papadokostaki, K.G.; Doghieri, F.; Minelli, M. A fundamental study of the extent of meaningful application of Maxwell's and Wiener's equations to the permeability of binary composite materials. Part III: Extension of the binary cubes model to 3-phase media. *Chem. Eng. Sci.* **2015**, *131*, 360–366. [CrossRef]

24. Wolf, C.; Guillard, V.; Angellier-Coussy, H.; Silva, G.D.G.; Gontard, N. Water vapor sorption and diffusion in wheat straw particles and their impact on the mass transfer properties of biocomposites. *J. Appl. Polym. Sci.* **2016**, *133*, 1–12. [CrossRef]

25. David, G.; Gontard, N.; Angellier-Coussy, H. Mitigating the Impact of Cellulose Particles on the Performance of Biopolyester-Based Composites by Gas-Phase Esterification. *Polymers* **2019**, *11*, 200. [CrossRef]

26. Berthet, M.-A.; Angellier-Coussy, H.; Machado, D.; Hilliou, L.; Staebler, A.; Vicente, A.; Gontard, N. Exploring the potentialities of using lignocellulosic fibres derived from three food by-products as constituents of biocomposites for food packaging. *Ind. Crops Prod.* **2015**, *69*, 110–122. [CrossRef]

27. Wolf, C. Multi-Scale Modelling of Structure and Mass Transfer Relationships in Nano- and Micro- Composites for Food Packaging. Ph.D. Thesis, University of Montpellier, Montpellier, France, 2014. Available online: https://tel.archives-ouvertes.fr/tel-02178975 (accessed on 7 July 2021).

28. Berthet, M.-A. Biocomposites Durables pour l'Emballage Alimentaire, À Partir de Sous-Produits de l'Industrie Alimentaire: Relations Procédé-Structure-Propriétés. Ph.D. Thesis, University of Montpellier, Montpellier, France, 2014.

29. Tschopp, M.A. 3-D Synthetic Microstructure Generation with Ellipsoid Particles; Technical Report; Defense Technical Information Center: Fort Belvoir, VA, USA, 2017. [CrossRef]
30. COMSOL Multiphysics. Separation Through Dialysis. 2008. Available online: https://www.comsol.ch/forum/thread/attachment/16637/dialysis_comsol-2156.PDF (accessed on 7 July 2021).

31. Zid, S.; Zinet, M.; Espuche, E. 3D Mass Diffusion in Ordered Nanocomposite Systems: Finite Element Simulation and Phenomenological Modeling. *J. Polym. Sci. Part B Polym. Phys.* **2019**, *57*, 51–61. [CrossRef]

32. Liang, M.; Feng, K.; He, C.; Li, Y.; An, L.; Guo, W. A meso-scale model toward concrete water permeability regarding aggregate permeability. *Constr. Build. Mater.* **2020**, *261*, 120547. [CrossRef]

33. Hedenqvist, M.; Gedde, U.W. Diffusion of small-molecule penetrants in semicrystalline polymers. *Prog. Polym. Sci.* **1996**, *21*, 299–333. [CrossRef]

34. Trifol, J.; Plackett, D.; Szabo, P.; Daugaard, A.E.; Baschetti, M.G. Effect of Crystallinity on Water Vapor Sorption, Diffusion, and Permeation of PLA-Based Nanocomposites. *ACS Omega* **2020**, *5*, 15362–15369. [CrossRef]

35. Kim, J.-K.; Hodzic, A. Nanoscale characterisation of thickness and properties of interphase in polymer matrix composites. *J. Adhes.* **2003**, *79*, 383–414. [CrossRef]

36. Joliff, Y.; Rekik, W.; Belec, L.; Chailan, J.F. Study of the moisture/stress effects on glass fibre/epoxy composite and the impact of the interphase area. *Compos. Struct.* **2014**, *108*, 876–885. [CrossRef]