# On the Equivalent In-Plane Permeability

Francisco Chinesta, Arnaud Poitou, Sylvain Chatel, Juan Antonio Garcia, Llanos Gascón

► To cite this version:

Francisco Chinesta, Arnaud Poitou, Sylvain Chatel, Juan Antonio Garcia, Llanos Gascón. On the Equivalent In-Plane Permeability. 13th ESAFORM Conference on Material Forming, 2010, Brechia, Italy. pp.651-654, 10.1007/s12289-010-0854-5. [hal-01008592](https://hal.science/hal-01008592)

HAL Id: hal-01008592

https://hal.science/hal-01008592

Submitted on 23 Apr 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
# ON THE EQUIVALENT IN-PLANE PERMEABILITY

F. Chinesta¹*, A. Poitou¹, S. Chatel², J.A. Garcia³, LL. Gascon³

¹GEM, Centrale de Nantes, BP 92101, F-44321 Nantes cedex 3, France

²EADS France Innovation Works, Allée du Chaffault, 44340 Bouguenais, France

³Universidad Politécnica de Valencia, Camino de Vera s/n, E-46071 Valencia, Spain

**ABSTRACT:** Simulation of RTM processes is usually performed in 2D. Thus, only a mesh of the middle plane is needed with the associated degrees of freedom savings with respect to fully 3D modeling. However such modeling needs the definition of an equivalent in-plane permeability representing the ignored dimension (the thickness). The definition of such permeability is not a trivial task because each ply in the thickness direction can be anisotropic, being the principal anisotropy direction different from one ply to the neighbor plies. In this work we propose a novel fully 3D modeling whose computational cost is equivalent to a 2D solution. It allows addressing properly the equivalent in-plane permeability issue.

**KEYWORDS:** RTM, Permeability, Separated representations, Model reduction

## 1 INTRODUCTION

In general the simulation of RTM processes assumes a 2D flow model. The most usual model results of combining the Darcy's law and the flow incompressibility:

$$ 
\begin{cases} 
\mathbf{v} = -\mathbf{K} \cdot \nabla p \\ 
\nabla \cdot \mathbf{v} = 0 
\end{cases} 
$$

that results in the second order BVP:

$$ \nabla \cdot (\mathbf{K} \cdot \nabla p) = 0 $$

The main issue in defining this model concerns the definition of the permeability tensor **K**. Different techniques exist, but in that follows we are assuming that an averaged permeability has been determined for each type of reinforcement architecture.

First we are considering a laminate composed of many layers with the same principal directions, whose permeability tensor writes

$$ \mathbf{K} = \begin{bmatrix} \alpha & \beta & 0 \\ \beta & \gamma & 0 \\ 0 & 0 & \zeta \end{bmatrix} $$

Now, we consider the flow in a plate mould, being the permeability principal directions parallel to the mould walls. In this case the permeability tensor becomes diagonal. The simplest scenario concerns $\alpha = \gamma$ and $\beta = 0$. In this case the flow is unidirectional. When the ply is placed after performing a certain in-plane rotation to the principal permeability directions, the $\beta$

component of the permeability tensor is no more zero. In that case when we impose a constant injection pressure, being the pressure in the opposite mould boundary zero, the resulting pressure field is depicted in Fig. 1.

<figure><img src="image_2.png" /><figcaption>**Figure 1: Pressure distribution**</figcaption></figure>

By computing the velocity fields associated with the pressure fields depicted in Fig. 1 we can notice significant deviation with respect to the unidirectional flow.

Now, we are analyzing another 2D situation, in which a laminate composed of different layers is considered. The principal directions of all the plies correspond to the coordinate axes, but the permeability is changing from one layer to other. Figure 2 depicts $\alpha(z)$, being

<footnote>* Francisco Chinesta: Centrale Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, Francisco.Chinesta@ec-nantes.fr</footnote>
$β = 0$, $γ$ and $ζ$ constants. Figure 2 depicts the evolution of $α(z)$.

<figure><img src="image_3.png" /><figcaption>**Figure 2: Evolution of the permeability in the thickness**</figcaption></figure>

Figure 3 depicts the pressure distribution related to the choice of the permeability. The velocity fields associated to this pressure distribution is depicted in Fig. 4

<figure><img src="image_4.png" /><figcaption>**Figure 3: Pressure distribution**</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4: Velocity field**</figcaption></figure>

We can conclude that if the principal directions of the permeability of each layer are aligned along the coordinate axes, even if the permeability evolves in the thickness direction, the flow remains unidirectional. The situation is radically different if the permeability also evolves in the x-direction. When we assume an evolution of $α(x,z)$ as depicted in figure 5, the pressure field that results, illustrated in figure 6, implies the velocity field depicted in figure 7, which exhibit noticeable deviations with respect to an unidirectional flow.

<figure><img src="image_6.png" /><figcaption>**Figure 5: Evolution of the permeability**</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6: Pressure field**</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7: Velocity field**</figcaption></figure>

These simple examples have proven that in general laminates, involving many layers with different reinforcement orientations the situation could be quite complex and the assumption of a simple average of the permeability of the different layers could be a crude approximation that deserves a deeper analysis. For this purpose a fully 3D simulation could be the best alternative, but an accurate representation of the thickness direction could become too expensive from the computational view point.

An appealing alternative consists in applying a separated representation of the pressure field involved in the porous media flow modelling. This kind of approximation allows computing 3D solutions with a computational cost characteristic of 2D simulations, that is, of standard RTM simulations. In the next section we
summarized the main elements of such technique, deeply described in many of our former works [1-2], and that is here particularized to the treatment of the Darcy's model defined in plate moulds. Its extension to shell geometries is quite direct.

## 2 PGD IN PLATE DOMAINS

In what follows we are illustrating the construction of the Proper Generalized Decomposition of a model defined in a plate domain $\Xi = \Omega \times I$ with $\Omega \subset \mathbb{R}^2$ and $I = [0, H]$:

$$ \nabla \cdot (\mathbf{K} \cdot \nabla p) = 0 \quad (1) $$

We consider that the laminate is composed of $P$ different anisotropic plies each one characterized by a well defined permeability tensor $\mathbf{K}_i(x, y)$-it is assumed constant in the ply thickness-. Moreover, without a loss of generality, we assume the same thickness for the different layers of the laminate, that we denotes by $h$. Thus, we can define a characteristic function representing the position of each layer:

$$ \chi_i(z) = \begin{cases} 1 & z_i \le z \le z_{i+1} \\ 0 & \text{otherwise} \end{cases}, \quad i=1,\dots,P \quad (2) $$

where $z_i = (i-1) \times h$. Now, the laminated permeability can be given in the following separated form:

$$ \mathbf{K}(x, y, z) = \sum_{i=1}^{i=P} \mathbf{K}_i(\mathbf{x}) \cdot \chi_i(z) \quad (3) $$

where $\mathbf{x} = (x, y) \in \Omega$.

The weak form of Eq. (1) writes:

$$ \int_{\Xi} \nabla p^* \cdot (\mathbf{K} \cdot \nabla p) d\Xi = 0 \quad (4) $$

with the test function $p^*$ in an appropriate functional space. The solution $p(x, y, z)$ is searched under the separated form:

$$ p(\mathbf{x}, z) \approx \sum_{i=1}^{i=N} X_i(\mathbf{x}) Z_i(z) \quad (5) $$

In what follows we are illustrating the construction of one such decomposition. For this purpose we assume that at iteration $n < N$ the solution is already known:

$$ p^n(\mathbf{x}, z) = \sum_{i=1}^{i=n} X_i(\mathbf{x}) Z_i(z) \quad (6) $$

and that at the present iteration we look for the solution enrichment:

$$ p^{n+1}(\mathbf{x}, z) = p^n(\mathbf{x}, z) + R(\mathbf{x}) \cdot S(z) \quad (7) $$

The test function involved in the weak form is searched under the form:

$$ p^*(\mathbf{x}, z) = R^*(\mathbf{x}) \cdot S(z) + R(\mathbf{x}) \cdot S^*(z) \quad (8) $$

By introducing Eqs. (7) and (8) into (4) it results:

$$ \begin{aligned} & \int_{\Xi} \left( \begin{pmatrix} \nabla R^* \cdot S \\ R^* \cdot \frac{dS}{dz} \end{pmatrix} + \begin{pmatrix} \nabla R \cdot S^* \\ R \cdot \frac{dS^*}{dz} \end{pmatrix} \right) \cdot \mathbf{K} \cdot \begin{pmatrix} \nabla R \cdot S \\ R \cdot \frac{dS}{dz} \end{pmatrix} d\Xi = \\ & = - \int_{\Xi} \left( \begin{pmatrix} \nabla R^* \cdot S \\ R^* \cdot \frac{dS}{dz} \end{pmatrix} + \begin{pmatrix} \nabla R \cdot S^* \\ R \cdot \frac{dS^*}{dz} \end{pmatrix} \right) \cdot \mathbf{Q}^n d\Xi \end{aligned} \quad (9) $$

where $\nabla$ denotes the plane component of the gradient operator $\nabla^T = (\partial/\partial x, \partial/\partial y)$ and $\mathbf{Q}^n$ denotes the flux at iteration $n$:

$$ \mathbf{Q}^n = \mathbf{K} \cdot \sum_{i=1}^{i=n} \begin{pmatrix} \nabla X_i(\mathbf{x}) \cdot Z_i(z) \\ X_i(\mathbf{x}) \cdot \frac{dZ_i(z)}{dz} \end{pmatrix} \quad (10) $$

Now, as the enrichment process is non-linear we propose to search the couple of functions $R(\mathbf{x})$ and $S(z)$ by applying an alternating direction fixed point algorithm. Thus, assuming $R(\mathbf{x})$ known, we compute $S(z)$, and then we update $R(\mathbf{x})$. The process continues until reaching convergence. The converged solutions allow defining the next term in the finite sums decomposition:

$$ R(\mathbf{x}) \rightarrow X_{n+1}(\mathbf{x}) \text{ and } S(z) \rightarrow Z_{n+1}(z). $$

We are illustrating each one of the just referred steps:

### 1. Computing $R(\mathbf{x})$ from $S(z)$:

When $S(z)$ is known the test function reduces to:

$$ p^*(\mathbf{x}, z) = R^*(\mathbf{x}) \cdot S(z) \quad (11) $$

and the weak form (9) reduces to:

$$ \begin{aligned} & \int_{\Xi} \left( \begin{pmatrix} \nabla R^* \cdot S \\ R^* \cdot \frac{dS}{dz} \end{pmatrix} \cdot \mathbf{K} \cdot \begin{pmatrix} \nabla R \cdot S \\ R \cdot \frac{dS}{dz} \end{pmatrix} \right) d\Xi = \\ & = - \int_{\Xi} \begin{pmatrix} \nabla R^* \cdot S \\ R^* \cdot \frac{dS}{dz} \end{pmatrix} \cdot \mathbf{Q}^n d\Xi \end{aligned} \quad (12) $$

Now, as all the functions involving the coordinate $z$ are known, they could be integrated in $I = [0, H]$. Thus, if we consider:

$$ \mathbf{K} = \begin{pmatrix} \mathbf{k} & \mathbf{k} \\ \mathbf{k}^T & \kappa \end{pmatrix} \quad (13) $$
$$ \text{with } \mathbf{k} = \begin{pmatrix} k_{xx} & k_{xy} \\ k_{yx} & k_{yy} \end{pmatrix}, \quad \mathbf{k} = \begin{pmatrix} k_{xz} \\ k_{yz} \end{pmatrix} \text{ and } \kappa = k_{zz}, \text{ then} $$

we can define:

$$ {}^x\mathbf{K} = \begin{pmatrix} \int_{z=0}^{z=H} k S^2 dz & \int_{z=0}^{z=H} \mathbf{k} \frac{dS}{dz} S dz \\ \int_{z=0}^{z=H} \mathbf{k}^T \frac{dS}{dz} S dz & \int_{z=0}^{z=H} \kappa \left(\frac{dS}{dz}\right)^2 dz \end{pmatrix} \quad (14) $$

and

$$ {}^x\mathbf{Q}^n = \sum_{i=1}^{i=N} \left( \begin{pmatrix} \int_{z=0}^{z=H} k S Z_i dz & \int_{z=0}^{z=H} \mathbf{k} S \frac{dZ_i}{dz} dz \\ \int_{z=0}^{z=H} \mathbf{k}^T \frac{dS}{dz} Z_i dz & \int_{z=0}^{z=H} \kappa \frac{dS}{dz} \frac{dZ_i}{dz} dz \end{pmatrix} \begin{pmatrix} \nabla R X_i(\mathbf{x}) \\ X_i(\mathbf{x}) \end{pmatrix} \right) \quad (15) $$

that allows writing equation (12) into the form

$$ \begin{aligned} & \int_{\Omega} \left( \begin{pmatrix} \nabla R^* \\ R^* \end{pmatrix} \cdot \begin{pmatrix} {}^x\mathbf{K} & \begin{pmatrix} \nabla R \\ R \end{pmatrix} \end{pmatrix} \right) d\Omega = \\ & \qquad = - \int_{\Omega} \left( \begin{pmatrix} \nabla R^* \\ R^* \end{pmatrix} \cdot {}^x\mathbf{Q}^n \right) d\Omega \end{aligned} \quad (16) $$

that defines an elliptic 2D problem defined in the middle plane of the plate.

## 2. Computing $S(z)$ from $R(x)$:

When $R(x)$ is known the test function writes:

$$ p^*(\mathbf{x}, z) = R(\mathbf{x}) \cdot S^*(z) \quad (17) $$

and the weak form (9) reduces to:

$$ \begin{aligned} & \int_{\Xi} \left( \begin{pmatrix} \nabla R \cdot S^* \\ R \cdot \frac{dS^*}{dz} \end{pmatrix} \cdot \begin{pmatrix} \mathbf{K} & \begin{pmatrix} \nabla R \cdot S \\ R \cdot \frac{dS}{dz} \end{pmatrix} \end{pmatrix} \right) d\Xi = \\ & \qquad = - \int_{\Xi} \left( \begin{pmatrix} \nabla R \cdot S^* \\ R \cdot \frac{dS^*}{dz} \end{pmatrix} \cdot \mathbf{Q}^n \right) d\Xi \end{aligned} \quad (18) $$

Now, as all the functions involving the in-plane coordinates $\mathbf{x} = (x, y)$ are known, they could be integrated in $\Omega$. Thus, using the previous notation, we can define:

$$ {}^z\mathbf{K} = \begin{pmatrix} \int_{\Omega} (\nabla R) \cdot (\mathbf{k} \cdot (\nabla R)) d\Omega & \int_{\Omega} (\nabla R) \cdot \mathbf{k} R d\Omega \\ \int_{\Omega} \mathbf{k} \cdot (\nabla R) R d\Omega & \int_{\Omega} \kappa(R)^2 d\Omega \end{pmatrix} $$

$$ {}^z\mathbf{Q}^n = \sum_{i=1}^{i=N} \left( \begin{pmatrix} \int_{\Omega} \nabla R \cdot (\mathbf{k} \cdot \nabla X_i) d\Omega & \int_{\Omega} \nabla R \cdot \mathbf{k} X_i d\Omega \\ \int_{\Omega} \mathbf{R} \mathbf{k} \cdot \nabla X_i d\Omega & \int_{\Omega} \kappa R X_i d\Omega \end{pmatrix} \begin{pmatrix} Z_i \\ \frac{dZ_i}{dz} \end{pmatrix} \right) \quad (20) $$

and

$$ {}^z\mathbf{Q}^n = \sum_{i=1}^{i=N} \left( \begin{pmatrix} \int_{\Omega} \nabla R \cdot (\mathbf{k} \cdot \nabla X_i) d\Omega & \int_{\Omega} \nabla R \cdot \mathbf{k} X_i d\Omega \\ \int_{\Omega} \mathbf{R} \mathbf{k} \cdot \nabla X_i d\Omega & \int_{\Omega} \kappa R X_i d\Omega \end{pmatrix} \begin{pmatrix} Z_i \\ \frac{dZ_i}{dz} \end{pmatrix} \right) \quad (20) $$

that allows writing equation (18) into the form

$$ \begin{aligned} & \int_{i} \left( \frac{S^*}{dS} \right) \cdot \begin{pmatrix} {}^z\mathbf{K} & \begin{pmatrix} S \\ \frac{dS}{dz} \end{pmatrix} \end{pmatrix} dz = \\ & \qquad = - \int_{i} \left( \frac{S^*}{dS} \right) \cdot {}^z\mathbf{Q}^n dz \end{aligned} \quad (21) $$

that defines a one-dimensional BVP.

## 3 CONCLUSIONS

This paper comes back to a recurrent issue in the numerical modelling of RTM flows, the one related to the pertinence of using 2D flow models making use of an averaged permeability of the different layers involved in the laminate. In complex situations significant deviations could be found, these deviations are being studied at present, and could justify the use of a fully 3D modelling. However, 3D simulations are reputed expensive from the point of view of the computational resources required for addressing complex scenarios. The use of separated representations as the ones involved in the proper generalized decompositions -PGD- could be an appealing alternative for addressing 3D models with a cost characteristic of 2D simulations. The application of the PGD on the fully 3D simulation of flows encountered in RTM processes constitutes a work in progress.

## REFERENCES

[1] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, J. Non-Newtonian Fluid Mech., 139: 153-176, 2006.

[2] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations, J. Non-Newtonian Fluid Mech., 144: 98-121, 2007.