# Three-dimensional reconstruction by dual kriging of the geometry of woven fibrous reinforcements in composite materials

Anna Madra, Piotr Breitkopf, François Trochu

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

## To cite this version:

Anna Madra, Piotr Breitkopf, François Trochu. Three-dimensional reconstruction by dual kriging of the geometry of woven fibrous reinforcements in composite materials. 12e colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01400457

HAL Id: hal-01400457

https://hal.science/hal-01400457

Submitted on 22 Nov 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Three-dimensional reconstruction by dual kriging of the geometry of woven fibrous reinforcements in composite materials

A. Madra<fnref n="1" />,<fnref n="2" />, P. Breitkopf<fnref n="1" />, F. Trochu<fnref n="2" />

**Résumé** — Composite materials with woven fibrous reinforcements are difficult to model numerically because of their complex internal geometry. The design of composite components requires modeling of the fabrication process and of the mechanical behavior. This creates a need for accurate three-dimensional geometrical models in order to improve the predictions of the fabrication process including fiber deformations during draping and closing of the mold. This work provides a solution based on dual kriging with nugget effect to create a parametric model of the surfaces of fiber tows from the voxel data obtained by X-ray microtomography.

**Mots clés** — woven composites, X-ray microtomography, mesh, numerical methods, Big Data.

## 1 Introduction

Composite materials based on a woven reinforcement and polymer matrix have been present in aero-nautics industry for a while now. Despite their superior mechanical performance and flexibility in design when compared to traditional alloys, their use in high volume production sectors like automotive have until recently been very limited. This is especially true for the high performance composites. The reason for this is not only the cost of base materials compared to eg. composites used in civil engineering but also the high cost of development and of optimization of production process.

Liquid Composite Molding (LCM) is a group of high volume production methods for high performance composites. It has been largely automatized but some critical problems have yet to be addressed. In LCM the mold containing a dry fiber preform – the reinforcement, is impregnated with liquid resin, which after curing becomes a matrix. From the fluid mechanics point of view, this presents a problem of multiphase flow through dual-scale porous media. Multiphase because we are dealing with a resin that is undergoing a curing reaction and which is also interacting with air or other gas present in the mold. Dual-scale porosity is the result of morphology of textile preforms. When the base material is carbon or glass fiber, then a preform has a mesoscale geometry of fiber tows with millimeter size diameters. The tows are in turn composed of thousands of single fibers of diameters from several to ca. 15 µm. This means that we are dealing with two porous media, one with mesoscale porosity of spaces between fiber tows and layers of textile and the second related to microscale porosity between the individual fibers inside the tows.

The process of impregnation of porosity at both scales depends on different flow mechanisms. At mesoscale, the cavities are filled by viscous flow and at microscale by capillary flow. The dominating type of flow is related to the flow injection rate as has been shown by [5]. What follows is that a deficiency in flow of one type will lead to lacks of impregnation at one scale or the other and to formation of voids. This is of utter importance as it has been shown in [3] that the volume fraction, type and distribution of voids can seriously impact the mechanical properties of a final composite material.

In conclusion, the geometry of the fiber reinforcement plays a major role in definition of the actual properties of the final composite as compared to those that are theoretically possible. Additional factors, like fiber bed deformation while draping, i.e., arranging the textile preform in the mold and when compacting the preform after mold closure also strongly impact the meso- and microscale geometry of the reinforcement and thus will change the required production parameters. Large scale production of high performance composites requires means to simulate those processes. Our goal is to provide tools to

<footnote n="1">¹ CNRS UMR7337 Roberval, Centre de Recherche Royallieu, Sorbonne universités, Université de technologie de Compiègne, {anna.madra,piotr.breitkopf}@utc.fr</footnote>

<footnote n="2">² Chaire sur les Composites à Haute Performance (CCHP), École Polytechnique de Montréal, {anna.madra, trochu}@polymtl.ca</footnote>
reconstruct fiber geometry suitable for further modeling and simulation.

The problem of textile reconstruction has been already addressed by [7,1]. These reconstructions are based on theoretical approximations of fiber and textile shapes. This leads to an error that may be too substantial to plausibly model permeability and impregnation phenomena. On the other hand, geometry reconstruction based on X-ray microtomography like presented in [3] is very accurate but leads to formation of meshes of very high space complexity, inadmissible for mechanical nor fluid mechanics calculations. The solution to this problem could be found by combining the advantages of both approaches. Our idea was to extract the relevant information on the textile geometry from the big data sets obtained with X-ray microtomography and reconstruct the fiber tow geometry with dual kriging algorithms. This allows to further reconstruct the relevant geometry at different precision levels, depending on the application. The method permits to adjust not only the precision of the reconstruction but also the smoothing of the final curves and surfaces by the use of nugget effect. The approach presented here has also been semi-automatized to treat big data sets with only limited initial input from the user.

# 2 Dual kriging

## 2.1 Kriging of curves

Kriging is an interpolation and estimation method especially popular in geostatistics to assess eg. ore distribution or population characteristics. Dual kriging has not yet been used extensively in geometrical modeling, although it allows creating automatically continuous and derivable models of curves, surfaces and solids. This also extends to reconstruction of closed contours and solids of revolution which is realized by parameterization. In dual kriging formulation a given function $U(x) \in \Re$ may be represented as [8]

$$U(x) = a(x) + W(x) \quad (1)$$

where $a(x)$ is the general trend of the function and $W(x)$ represents the fluctuations. The general trend is the averaged behavior of a function described with the variance of the estimation error.

We assume that the general trend follows a linear function

$$a(x) = a_0 + a_1 \cdot x \quad (2)$$

Then modeling a one-dimensional function (simple curve) with kriging can be considered as a minimization problem of

$$\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_i(x) \cdot \lambda_j(x) \cdot K(|x_i - x_j|) - \sum_{i=1}^{N} \lambda_i(x) \cdot K(|x_i - x|) \quad (3)$$

with linear unbiased criteria

$$\begin{cases} \sum_{i=1}^{N} \lambda_i - 1 = 0 \\ \sum_{i=1}^{N} \lambda_i \cdot x_i - x = 0 \end{cases}$$

where $\lambda_i \in \Re^N$ are coefficients in the minimization problem of the expression (1) under unbiased conditions and $K$ is the generalized covariance. The solution of the (3) results with the one-dimensional kriging system which could be further transformed into dual kriging system (for a function with linear trend) [8]

$$\begin{bmatrix} K(|x_i - x_j|) & 1 & x_1 \\ & \vdots & \vdots \\ & 1 & x_N \\ 1 & \dots & 1 & 0 & 0 \\ x_1 & \dots & x_N & 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} b_1 \\ \vdots \\ b_N \\ a_0 \\ a_1 \end{bmatrix} = \begin{bmatrix} u_1 \\ \vdots \\ u_N \\ 0 \\ 0 \end{bmatrix} \quad (4)$$

where $x_i$ are the coordinates of the known values $u_i$ of the original function $U(x)$. The system is solved for values $a_i$ and $b_i$. They are then used as coefficients in the final dual kriging equation of an interpolated function $u(x)$

$$u(x) = a_0 + a_1 \cdot x + \sum_{i=1}^{N} b_i \cdot K(|x - x_i|) \quad (5)$$
Generalized covariance K is also called a shape function. In practical implementations shape func-
tions are chosen so as to give 0 values for the same coordinate (i.e. 0 on the diagonal of the matrix in 4).
The common examples are

- linear $K(h) = h$

- cubic $K(h) = h^3$

- logarithmic $K(h) = h^2 \ln(h)$

- sinusoidal $K(h) = \sin(h)$

Functions for general trends are also a subject of choice and are solved under adequate unbiased criteria.
For one-dimensional kriging this gives the following formulations :

- constant $a(x) = a_0$

- linear $a(x) = a_0 + a_1 x$

- square $a(x) = a_0 + a_1 x + a_2 x^2$

- trigonometric $a(x) = a_0 + a_1 \cos(\omega x) + a_2 \sin(\omega x)$

The above formulation works only for curves representable by functions in Cartesian plane. In case of
curbs and closed contours a parametric formulation is required. The original points $P_i$ describing the curve
can be assigned to a parameter $t_i$ with discrete values from the interval $[0, 1]$. And the dual kriging system
from eq. 4 is solved for a matrix with the number columns equal to the number of coordinates that are
parameterized. The interpolated function has a form analogous to that in eq. 5

$$
x(t) = a_0^x + a_1^x \cdot t + \sum_{i=1}^{N} b_i^x \cdot K(|t - t_i|) \quad (6)
$$

This formulation after suitable modifications can be also applied to model surfaces and three-dimensional
objects. Further modifications to the dual kriging system may be introduced to treat special cases, like
discontinuities by kriging with derivatives. This problem though has been deemed irrelevant in current
case of textile geometry reconstruction where smoothing is one of the goals and discontinuities are
avoided.

2.2 Kriging of surfaces

Extension of the dual kriging system to the second dimension requires different approaches depen-
ding on the case if the surface is described without or with parameterized functions. In the first case,
the Euclidean distance between x and y coordinates of a given point P can be used to solve the kriging
system under modified unbiased conditions. In the parametric case, the problem of surface modeling can
be divided into two subproblems. The first is the representation of a local cross-section (usually parame-
trized) solved by kriging along the profile no. 1 and the second is the interpolation of the curve along
which those cross-sections shall be placed – kriging along the profile no. 2. For the fiber tow case, this
would mean that at each point on the tow-length there is a different kriging system to solve pertaining to
given tow cross-section. There is also another kriging system to calculate the paths joining corresponding
points on the neighboring cross-sections.

First, let's assume that profile no. 1 is parameterized with $s \in [0, M]$ with shape function $K_1$ and constant general trend. We parameter profile no. 2 by $t \in [0, N]$ with shape function $K_2$ and linear general trend. Then for each $t_j$ there we formulate a dual kriging system like in 4 with kriging matrix $S$. Idem for each $s_i$ where we rename the kriging matrix $T$. Further we construct an auxiliary matrix $P_x$

$$
P_x = \begin{bmatrix}
x(s_1, t_1) & \dots & x(s_1, t_N) & 0 & 0 \\
\vdots & \ddots & \vdots & 0 & 0 \\
x(s_M, t_1) & \dots & x(s_M, t_N) & 0 & 0 \\
0 & \dots & 0 & 0 & 0
\end{bmatrix} \quad (7)
$$

and vectors $k_1(s)$ and $k_2(t)$

$$
\begin{align}
k_1(s) &= [K_1(|s - s_1|) \quad \dots \quad K_1(|s - s_M|) \quad 1] \notag \\
k_2(t) &= [K_2(|t - t_1|) \quad \dots \quad K_2(|t - t_N|) \quad 1 \ t] \tag{8}
\end{align}
$$

Then an expression for the x coordinate of the interpolated function is

$$
x(s,t) = k_1(s)^T \cdot S^{-1} \cdot P_x \cdot T^{-1} \cdot k_2(t) \quad (9)
$$
and idem for coordinates y and z. This approach has been employed to model fiber tows as their cross-section always requires parameterization. In the case of complex three-dimensional weaves, the parameterization along the tow length may also be required.

## 2.3 Nugget effect

An important property of an interpolation function obtained with dual kriging is that it will pass through all of the original points used to model given geometry. This can be adjusted by introducing a nugget effect $\sigma^2$. The value attributed to $\sigma^2$ is usually the variance in measurement of the original geometry. The nugget effect may be set constant for the whole data set or adjusted for each measurement point. It is introduced into dual kriging system by adding it to the diagonal (shape functions are 0 for $i = j$) of the matrix in 4

$$ \begin{bmatrix} \vdots & K(|x_i - x_j|) & \dots & 1 & x_1 \\ K(|x_i - x_j|) & \sigma_i^2 & \vdots & \vdots & \vdots \\ \vdots & \vdots & \ddots & 1 & x_N \\ 1 & \dots & 1 & 0 & 0 \\ x_1 & \dots & x_N & 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} b_1 \\ \vdots \\ b_N \\ a_0 \\ a_1 \end{bmatrix} = \begin{bmatrix} u_1 \\ \vdots \\ u_N \\ 0 \\ 0 \end{bmatrix} \quad (10) $$

In case of fiber geometry reconstructions, nugget effect can be used to control the final area of the fiber tow cross-section. The minimal area occupied by a fiber tow is the solution of the circle packing problem in the plane. Due to the method of production, glass and carbon fibers can be considered as uniform in size and shape. Thus, the highest density of packing in a plane is as given by Joseph Lagrange for hexagonal lattice

$$ \eta_h = \frac{\pi}{2\sqrt{3}} \approx 0.9069 \quad (11) $$

It follows that given a certain fiber volume, the final tow area has to be increased by a factor of $\approx 1.1027$.

Another use of nugget effect is in handling of the segmentation error of X-ray tomograms. This allows to individually adjust every vertex of the modeled geometry. The information on original fiber content and segmentation is obtained during data preprocessing step which is described in the next section.

## 3 Data preprocessing

The material used in this study was supplied by Ecole des Mines de Douai, France. A laminate of thermoplastic reinforced with 2/2 twill weave glass textile of Tex=266 g/km has been tested. The X-ray microtomography scans were performed in Laboratoire Mateis of INSA Lyon, France with the following acquisition parameters :

- Focus-to-Detector Distance (FDD) : 577.3 mm

- Focus-to-Object Distance (FOD) : 7.7 mm

- CCD camera effective area : 1920 x 1536 pixels

- magnification : 74.7 x

- voxel resolution : 1.7 x 1.7 x 1.7 $\mu$m$^3$

The sinograms acquired during X-ray scanning were then transformed by backprojection filtering algorithms to obtain final tomograms like the one shown in Fig. 1. This data has then been analyzed in ImageJ/Fiji [6] to yield the following principal characteristics of the material. For an example glass fiber tow

- total cross-sectional area of glass filaments : 99 127 $\mu$m$^2$

- average diameter of a single glass filament : 6.9229 $\mu$m

- total number of glass fiber filaments in a fiber tow : 658

Then according to the Lagrange criterium 11 for hexagonal lattice, the minimal area of a fiber tow is

$$ 99 127 \cdot 1.1027 = 109 303 \, \mu\text{m}^2 \quad (12) $$

This value will be further used to adjust the reconstruction parameters by nugget effect. After the initial measurements, the goal of data preprocessing is to extract the data describing the contour of the given
tow cross-section and prepare them to input into kriging algorithms. To achieve this the following image processing approach has been developed. The whole process is illustrated in Fig. 2. First the selected fiber tow is classified using Fast Random Forest algorithm. The results are shown in Fig. 2a. This image shows the probability of classification of a given pixel to the fiber tow class. The possible values span from 0 – not in the tow, to 255 – definitely in the tow. They will be further used to adjust the interpolation via nugget effect.

The probability image is then copied and its local thickness is being calculated with the result shown in Fig. 2b. This step enables to further automatically detect the outer contour of the tow by local automatic thresholding operation which results in filled cross-sections of the fiber tow presented in Fig. 2c. Finally, smaller areas which do not belong to the tow are discarded and the result is skeletonized to provide 1 pixel thick outline of a fiber tow in form of a binary image shown in Fig. 2d, which can be readily exported as a list of input vertexes for the kriging algorithm.

It has be mentioned in this place that all of the operations may introduce considerable error in contour definition. For this reason during extraction of vertex data, the probability of tow classification at these points is also added to enable further modifications to the contour during reconstruction by the use of nugget effect. This would mean that the actual error of the reconstruction will be directly tied to the quality of the classifiers trained with the Fast Random Forest algorithm. Given that the mesoscale morphology of the fiber tow can be a subject to some arbitrary assumptions and the same situation may occur for establishing the criteria for segmentation, this type of error is deemed acceptable for the further analyses.

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – An untreated tomogram image of glass fiber tow cross-section obtained from the X-ray microtomography. The intensity levels in the image indicate different coefficients of linear accentuation of different phases present in the material.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2 – Steps of the image processing approach developed to extract fiber tow contour data.</figcaption></figure>

## 4 Results

The final kriging library and program for parametric curve and surface reconstruction has been developed entirely in Python 2.7.6 using 64-bit deployment. First the approach has been tested for parametric kriging of a single fiber tow cross-section outline based on the data obtained from the curve extracted in Fig. 2d. The curve has been approximated using cubic shape function and constant general trend. The nugget effect has been set to 0. The results are shown in Fig. 3 for different discretizations of the *t* parameter. It may be observed that for the minimal amount of points *N* = 5 (Fig. 3a) where $N = \|[0, 1]\|$, the
curve is a convex quadrangle respecting the maximal dimensions of the original curve. With $N = 20$ (Fig. 3c) concave features begin to appear and for $N = 100$ (Fig. 3e) the error of reconstruction approaches 0. No significant improvements were noted for $N = 1000$ (Fig. 3f).

In the above reconstructions the error of segmentation has not been taken into account. For the same initial contour data (Fig. 2d) a parametric kriging system with cubic shape function and constant general trend has been used but with different strengths of the nugget effect. The results are shown in Fig. 4. The nugget effect is used to adjust the contour shape at each point individually based on the segmentation error. The impact of this adjustment on the overall shape is controlled by changing the strength of the nugget effect by multiplying its value by a constant. The greater the strength of the nugget the value, the higher is the degree of smoothing of the final curve. It may be observed that with strength 100 the reconstructed shape differs significantly from the original. To eliminate deviations from the original we calculate the reference between the original cross-section area which is taken for the nugget effect strength 1 and that of the given reconstructed contour. For all the examples shown here, only the curve with the highest degree of smoothing 100 was 6.5% smaller than the original area, whereas the rest of the contours deviated at most by 1.1%.

The final application of the presented method was parametric kriging of the whole fiber tow surface. The geometry has been reconstructed from 23 tomogram slices taken at different points on the tow length. The reconstruction used cubic shape function and constant general trend. The reconstruction has been performed for different strengths of nugget effect which was imposed arbitrarily. All reconstructions are shown in Fig. 5. As previously with curves, the final mesh is considerably smoother when the $\sigma^2$ is increased. To control the quality of surface reconstruction, the difference in total volume encompassed by the surface can be taken into account. As here the value of nugget effect has been chosen arbitrarily, this comparison has been deemed unjustified from the topological point of view. To compare the surface mesh obtained with kriging to the original, based on microtomography, both have been superposed in three-dimensional domain and rendered as shown in Fig. 6. The nugget effect was chosen high enough to enable maximum smoothing while still preserving important features like fiber tow thinning.

<figure><img src="image_5.png" /><figcaption>FIGURE 3 – Contour reconstruction with parametric kriging for different discretizations of $t$ parameter with $N = |[0, 1]|$. Nugget effect $\sigma^2 = 0$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 4 – Contour reconstruction with parametric kriging for $N = 250$ with different strengths of nugget effect $\sigma_i^2$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 5 – Tow surface reconstruction with parametric kriging for $N = 12$ and $M = 100$ (discretizations along profiles no. 1 and 2) with different values of constant nugget effect $\sigma^2$.</figcaption></figure>

# 5 Conclusion

An automated solution for geometry reconstruction from raw X-ray microtomography scans has been presented. It permits to either automatically (via cross-section area or volume control) or based on the user's decision adjust the geometry of the reconstructed shape creating mesh data. The accuracy of the final surface mesh is fully determined by the user and the final purpose of the reconstructed mesh, from the smoothed outlines suitable for textile compression analysis to meshes including complex microscale features interesting in studies on permeability and composite impregnation.
<figure><img src="image_8.png" /><figcaption>FIGURE 6 – Comparison of untreated mesh of fibers tows with smooth geometry envelopes obtained with dual kriging. Nugget effect $σ^2=0.01$, constant general trend, cubic shape function, parameter discretization : N = 12, M = 30.</figcaption></figure>

The analysis of algorithm performance has shown that parametric kriging of contours is a very interesting method that can be further coupled with other approaches to quickly generate modifiable meshes. The parametric surface generation was deemed too computationally expensive to find practical applications although it can be used for benchmarking purposes as it is already included in the current implementation. In any case, the final mesh is definitely more suitable for further numerical computations than the one obtained in previous works from the segmented tomography images.

## Acknowledgements

The authors would like to thank Jerôme Adrien and Prof. Eric Maire from Laboratoire Mateis at INSA Lyon, France for the help in performing the X-ray micromotomography scans. Further thanks go to the Centre d'Innovation de l'Université de Technologie de Compiègne, France for financing the X-ray microtomography testing. Thanks go also to the staff of Chaire sur les Composites à Haute Performance (CCHP) at École Polytechnique de Montréal, Canada for valuable insights into the LCM processes of composite manufacture and explanation of the permeability phenomena in fibrous preforms.

## Références

[1] L. P. Brown, X. Zeng, A. Long, and I. A. Jones. *Recent developments in the realistic geometric modelling of textile structures using TexGen*, International Conference on Digital Technologies for the Textile Industries, 2013.

[2] S. Lomov, D. Ivanov, I. Verpoest, M. Zako, T. Kurashiki, H. Nakai, and S. Hirosawa. *Meso-FE modelling of textile composites: Road map, data flow and algorithms*, Composites Science and Technology, Elsevier, 1870-1891, 2007.

[3] A. Madra, N. El Hajj and M. Benzeggagh. *X-ray microtomography applications for quantitative and qualitative analysis of porosity in woven glass fiber reinforced thermoplastic*, Composites Science and Technology, Elsevier, 50-58, 2014.

[4] J. Pazmino, V. Carvelli, and S. V. Lomov. *Micro-CT analysis of the internal deformed geometry of a non-crimp 3D orthogonal weave E-glass composite reinforcement*, Composites Part B: Engineering, Elsevier, 2013.

[5] E. Ruiz, V. Achim, S. Soukane, F. Trochu, and J. Breard. *Optimization of injection flow rate to minimize micro/macro-voids formation in resin transfer molded composites*, Composites Science and Technology, Elsevier, 475-486, 2006.

[6] J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, and S. Saalfeld. *Fiji: an open-source platform for biological-image analysis*, Nature Methods, Macmillan Publishers Limited, 676-682, 2012.

[7] I. Verpoest and S. Lomov. *Virtual textile composites software: Integration with micro-mechanical, permeability and structural analysis*, Composites Science and Technology, Elsevier, 2563-2574, 2005.

[8] F. Trochu. *Krigeage dual et applications en modélisation géométrique*, Notes du cours MEC6310, École Polytechnique de Montréal, 2003.