JOURNAL OFCOMPOSITE
M AT E R I A L S

ArticleStiffness and permeability multi-objective
optimization of carbon-fiber-reinforced
plastic mesostructures using
homogenization method

Journal of Composite Materials
2019, Vol. 53(13) 1865–1880! The Author(s) 2018
Article reuse guidelines:
sagepub.com/journals-permissions
DOI: 10.1177/0021998318812160
journals.sagepub.com/home/jcm

Ryosuke Matsuzaki1,2 , Tomohiro Ishikawa1,
Tomonaga Okabe3,4, Yutaka Oya3 and Shigeki Yashiro5

Abstract
This paper presents a stiffness and permeability multi-objective optimization method for carbon-fiber-reinforced plastic
mesostructures based on a homogenization method. To reduce the computational cost of dealing with multiple design
variables for complicated fiber mesostructures, we generate and extract effective design variables from optimization results
derived from a smaller number of design variables. We applied the proposed method to optimization of the in-plane andout-of-plane stiffness and permeability of non-crimp fabrics. The optimization results showed that the application ofeffective design variables enabled attainment of an improved elastic modulus and permeability. From investigations of theobtained optimized design, we clarified the trade-off relationship between the elastic modulus and permeability, andelucidated the effects of dimensions of non-crimp fabric mesostructures on the elastic modulus and permeability.

Keywords
Heterogeneous materials, homogenization, non-crimp fabric, permeability, porous media, shape optimization

Introduction
For fiber material used in the resin transfer molding(RTM) process of composites, optimization of themesostructure is important for improving the productivity while maintaining a constant level of quality. Studies
have been conducted on the elastic modulus and permeability coeﬃcient of fiber material mesostructures using
the homogenization method.1–4 However, to date, nostudies have been conducted that successfully realized
the optimization of multi-objective functions, such as
those related to the elastic modulus and permeability
coeﬃcient for fiber material mesostructures.
Although the location in which design variables are
set is an important factor for the optimization, thedetermination criteria remain unclear and are often
dependent on the experience of the designers. Studies
are therefore being conducted to achieve topology optimization5–9 that derives optimized shapes, even with noconsideration of determination criteria. In topology
optimization, there is no need to set design variables,
and optimized shapes can be derived by setting only theboundary conditions. However, there are no reported

examples of the application of topology optimization tofiber material mesostructures with restrictive conditions
for fabrication. This lack is attributed to the diﬃculty
associated with obtaining optimization solutions.
Topology optimization is suitable for developing thefundamental design during the early stages of the process. Nevertheless, shape optimization is considered tohave potential for implementation to practical design
applications. As a result, to derive optimized solutions,many researchers have performed shape optimizations
that incorporate several design variables that are
1

Department of Mechanical Engineering, Tokyo University of Science,Japan
2
Institute of Fluid Science, Tohoku University, Japan3
Department of Aerospace Engineering, Tohoku University, Japan
4
Department of Materials Science and Engineering, University ofWashington, USA
5
Department of Aeronautics and Astronautics, Kyushu University, Japan
Corresponding author:
Ryosuke Matsuzaki, Tokyo University of Science 2641 Yamazaki, Noda
Chiba, 278-8510 Japan.
Email: rmatsuza@rs.tus.ac.jp

1866
available in advance. To derive superior optimization
solutions in an eﬀective manner, design variables must
be reduced to eﬀective design variables.10–14 This isbecause the problematic aspect of increasing computational costs exists when the number of design variables
is increased. However, there are two primary problems
with previously used methods for optimization of fiber
material mesostructures.
First, in the extraction of eﬀective design variables, theoptimization results were analyzed by incorporating as
many design variables as possible into the optimization
model. Because models that incorporate a large numberof design variables have complex shapes, the mesh generation becomes a bottleneck during analyses owing to
the amount of time required. Second, the criteria for the
extraction of eﬀective design variables are arbitrary andnot based on any quantitative parameters.
In this paper, we therefore propose an optimization
method for the stiﬀness and permeability of carbonfiber-reinforced plastic (CFRP) mesostructures in aneﬃcient manner. The method generates eﬀective
design variables from optimization results that arederived from a smaller number of design variables.
Then, we examine the eﬀectiveness of the proposed
method by considering the optimization problems ofthe mesostructures of non-crimp fabric (NCF) material
arranged in molding during the RTM process. Finally,
we discuss the trade-oﬀ problems and dimension eﬀects
of multi-optimization of fiber material mesostructures.

Methods
Optimization outline
Figure 1 outlines the proposed method in comparison
with a conventional method. The proposed method

Journal of Composite Materials 53(13)involves the preparation of a simple fiber-material
mesostructure model that is expressed using theminimum number of design variables to minimize thecomputational time. We employ the homogenization
method to derive the elastic modulus and permeability
coeﬃcient of the fiber material mesostructure. To maximize the elastic and permeability coeﬃcients, we optimize the dimensions of the fiber-material mesostructureusing the non-dominated sorting genetic algorithm 2(NSGA-2),15 which is a multi-objective genetic algorithm (MOGA).
The optimization results are next used as sample
data. These data include new design variables, which
are dependent variables of the original design variables.
The latter variables were produced by dividing the optimization model into regions. Here, we verify the multicollinearity using a scatter plot matrix (SPM), whilevariables that have a high correlation coeﬃcient are
combined. After that, we generate the mathematicalformula for predicting the value of the objective functions. To this end, we use genetic programing (GP) toclarify the relationship between the design variables andthe objective functions, as well as to generate andextract eﬀective design variables for improving the optimization solutions. The mathematical formula we use
to predict the value of the objective functions is referred
to as the prediction model, as described below.
To determine the relationship between the design
variables and the objective functions, an evaluation
standard for the prediction model must be establishedwhen deriving a prediction model using the GP. We use
the GP to incorporate the concept of Akaike’s information criterion (AIC)16 and the tolerance (TOL)17
concept into the evaluation standard for deriving prediction models. The AIC indicates a good fit to thesample data and prevents overfitting throughFigure 1. Conceptual scheme of the multi-objective optimization of fiber-reinforced composites using the homogenization method.

Matsuzaki et al.

1867

minimization. Furthermore, we can use the TOL toexpress the relationship between the design variables
and the objective functions with a valid prediction
model because multi-collinearity, which renders as
uncertain the codes that are included in the prediction
model, is prevented. Terms that are included in the prediction model having large-valued coeﬃcients may be
extracted as eﬀective design variables.

The following equilibrium equation, strain–
displacement expression, and constitutive equation
derived from the mechanics of elasticity are used toderive the elastic modulus of the macrostructure elicitedby the mesostructure

Calculation of elastic and permeability coefficient for
the fiber-material mesostructure

1 @u!i @u!jþ"!kl ¼
2 @xj @xi

In the case of fiber-reinforced composites, Young’s
modulus and the shear modulus can be predictedusing Halpin–Tsai equations,18 while the permeability
coeﬃcient can be predicted using the Gebart method.19
However, these methods cannot predict the macroscopic properties of CFRP fabrics that have arbitrary
mesoscale structures. The homogenization method20
must therefore be used to connect macro- andmesoscales.
The homogenization method defines the macrostructure (region !) having the same size as the inhomogeneous overall structure of the target considered for the
analysis, as well as the periodical unit structure thathas micro-inhomogeneity, as the mesostructure (region
Y). The ratio of the representative length of the macrostructure l! to the representative length of mesostructurelY is also defined, as shown in equation (1) (see Figure 2).
lY
l!

!¼

ð1Þ

Furthermore, the coordinate of the macrostructureis defined by x, while the coordinate of the mesostructure is defined by y. The relationship between x and y isdescribed below using scale ratio k
y¼

x
!

ð2Þ

Homogeniza!onLocaliza!on

lW

lY
Y

W

x

y

Figure 2. Multiscale simulation of the elastic body using thehomogenization method.

@"ij!
¼0
@xi

ð3Þ
!

"ij! ¼ D!ijkl "!kl

ð4Þ
ð5Þ

where the physical quantities of the overall structure
depend on the mesostructure. Accordingly, the superscript k is added. Furthermore, u!i is the displacementof the overall structure. The principle of virtual work forthe overall structure shown below is derived by applying
to equation (3) the boundary conditions for the overall
structure as well as the Gauss divergence theorem
ZZ!
! @#ui"ij
d! ¼
ti #u!i d"
!
@x
!
j
!
"t

ð6Þ

where "!t is the mechanical boundary surface for theoverall structure, while ti is the surface force per unit
area.
Let us consider that the displacement u!i for the overall structure, which depends on the mesostructure, canbe expressed using an asymptotic expression, as shownby the following equation
u!i ¼ u0i ðxÞ þ !u1i ðx, yÞ þ !2 u2i ðx, yÞ þ % % %ﬃ u0i ðxÞ þ !u1i ðx, yÞð7Þ

where u0i is the displacement of the macrostructure
and u1i ðx, yÞ is the displacement disturbance of themesostructure. Equation (7) expresses the displacementof the overall structure, together with the mesostructure
displacement disturbance and macrostructure displacement. The mesostructure unit is assumed to be repeatedinfinitely using the homogenization method. In otherwords, we assign the extreme boundary of k! 0. The
partial diﬀerential of the overall structure then becomesequation (8), as shown below, by using the combinationof equations (2) and (7)@u!i@u0 @u1¼ i þ i@xj @xj @yj

ð8Þ

We apply the averaging approximation theorem
expressed by equation (9) to separate the principle of1868

Journal of Composite Materials 53(13)virtual work for the overall structure, as expressed byequation (6), into the integration domains of themacrostructure and mesostructurelim

Z!!0 !!

!

d! ¼

Z

Z
1
ðx, yÞdYd!
! jYj Yð9Þ

The integrand (x,y) represents an inhomogeneous
physical quantity and a material constant. For equation(9) to be valid, a periodicity must exist for Y.
We substitute equations (3), (5), (7), and (8) into
equation (6), and we apply the averaging approximation theorem to derive the macrostructural formula ofequation (10) and the mesostructural formula of equation (11).
! 0
"
Z
Z
Z
@uk @u1k @#u0i
1
Dijkl
þ
dYd! ¼
ti #u0i d"@xl @yl @xj! jYj Y"t
! 0"
Z
@uk @u1k @#u1i
1
Dijkl
þ
dY ¼ 0jYj Y@xl @yl @yj

ð10Þ

ð11Þ

We discretize equations (10) and (11) using the finiteelement method (FEM) based on the hypothesis thatthe mesostructural displacement disturbance u1i is proportional to the macrostructural strain, represented byequation (12), to derive the homogenized macrostructural elastic stiﬀness constant of equation (13)u1i ðx, yÞ ¼ $kl
i ðx, yÞDH
ijmn ¼

@u0k
@xl

Z

#
$
1
Dijmn ' Dijkl B1l% $mn
k% dY
jYj Y

ð12Þ

v!i ¼ 0

@v!i
¼0
@xi

ð14Þ

@p!
@2 v!iþ &!
¼0@xi@xj @xj

ð15Þ

on

"

ð16Þ

Here, when the extreme boundary of k! 0 is
assigned, &k is defined as&! ¼ ! 2 &ð17Þ

The above assumes that, when the extreme boundaryof k! 0 is assigned, the overall flow of the liquid
becomes slower with respect to k2. Let us consider
that the flow velocity and pressure for the overall structure can be expressed using an asymptotic expression,as was the case in the previous section
v!i ðx, yÞ ¼ v0i ðx, yÞ þ !v1i ðx, yÞ þ !2 v2i ðx, yÞ þ % % %ﬃ v0i ðx, yÞ þ !v1i ðx, yÞ

p! ðx, yÞ ¼ p0 ðx, yÞ þ !p1 ðx, yÞ þ !2 p2 ðx, yÞ þ % % %ﬃ p0 ðx, yÞ þ !p1 ðx, yÞ

ð18Þ
ð19Þ

where v ji and pi are periodic in Y. Based on the relationship in equation (2), we can express the chain ruleof diﬀerentiation as follows
%
%
%
@ %%
1 @ %%@ %%¼
þ
@xi %!! ! @yi %Yf @xi %!f

ð20Þ

@v0i
¼0
@yi

ð21Þ

ð13Þ

where superscript H signifies the homogenized macrostructural model. $kl
i ðx, yÞ is referred to as the characteristic displacement function, with periodicity in Y.
Subscript a in equation (13) represents the constituent
node number of an element.
Similar to the derivation of the elastic modulus, weconsider the Stokes flow of Newtonian fluid, and we
consider the solid phase to be undeformed to simplify
the derivation of the permeability coeﬃcient. There are
two governing equations, namely the equation of continuity and the balance equilibrium of the liquid phase,and we use the following equations'

where p represents the pressure, vi denotes the flow
velocity, and & represents the viscosity. Because thephysical quantities for the overall structure depend on
the mesostructure, we add superscript k. The boundary
condition for the solid and liquid boundaries is thenonslip conditionfEquations (14) and (15) must be valid for therespective orders of k. Equations (18) and (20) are
therefore substituted into equation (14) and rearranged
for the respective orders of k.

@v0i @v1i
þ¼0@xi @yiðorder of !'1 Þðorder of !0 Þ

ð22Þ

Rearranging equation (15) in a similar manner
results in the following
'
'

@p0
¼0
@yi

ðorder of !'1 Þ@p1 @p0@ @v0i'þ&¼0@yj @yj@yi @xi#$
order of !0

ð23Þ
ð24Þ

Matsuzaki et al.

1869

The following equation must be valid in order forequation (22) to have v1 as the solution on themesoscale:Z@v0i
dy ¼ 0
Y @xi

ð25Þ

Vi is the flow velocity on the macroscale and isdefined as
Vi ¼

Z
1
v0 dy
Y Y i

ð26Þ

Based on equations (25) and (26), it then becomes
evident that the equation of continuity on the macroscale must be valid, as shown by the following equation
@Vi
¼0
@xi

ð27Þ

From equation (23), we see that p0 is not a functionof y.
Let us assume that the hypothesis in equations (28)and (29) is valid based on the use of the characteristic
function21–23'j @p0v0i ¼ ' i
& @xj
p1 ¼ 'Aj# $@p0
þ p# 1 xj
@xj

ð28Þ
ð29Þ

where 'ji and Aj are characteristic functions withperiodicity in Y. In the same way, p denoted with abar is independent of yj. If we substitute equation (28)
into equation (21), and equations (28) and (29) aresubstituted into equation (24), we can derive thefollowing3,21,22,24,25'

permeability coeﬃcient tensor Kij is derived by taking
the volume average of 'ji in unit cells, as described
below
Kij ¼

Z
1
'j dY
Y Y i

ð33Þ

The validation of the homogenized permeability calculation is described in Appendix 2.

Genetic programming for the prediction model
For a good fit of the homogeneous elastic modulus andthe permeability coeﬃcient, we use the GP to optimize
the prediction model, which indicates the relationshipbetween the objective functions and design variables
with the added dependent variables. We derive thedependent variables by dividing the model into regions
and expressing the regions with combinations of design
variables. Furthermore, we normalize the respective
variables using equation (34) to assign a value within
the range of zero to one. This is because we compare
the impact of the design variables in terms of the coefficients when the eﬀective design variables are extracted
from the prediction modelx~ i ¼

ðxi ' xmin Þðxmax ' xmin Þ

ð34Þ

In the above equation, the variable denoted with atilde (() represents a normalized variable.
When deriving a prediction model using the GP, we
implement the AIC to a degree of fitness fi. We establish the AIC using the unbiased estimator of the
expected average logarithmic likelihood with reference
to the Kullback–Liebler (K-L) divergence, which indicates the quality of the prediction model. The AIC for
the multiple regression models is described below
because prediction models derived by the GP are multiple regression models

@'ji
¼0
@yi

ð30Þ

@Aj
@ @'jiþ¼ '#ij@yi @yk @yk

! !
"
"
Se
fi ¼ n ln 2(þ 1 þ 2ð m þ 2Þn

ð31Þ

where n represents the number of samples, Se denotes
the sum-of-squares error, and m represents the numberof explanatory variables. The first term on the rightside of equation (35) signifies the size of the error forthe prediction model and the target of prediction, while
the second term signifies the complexity of the equation
for the prediction model. We can derive the intended
prediction model by minimizing fi.
Furthermore, there exists a problem of multicollinearity upon generating a prediction model.
Multi-collinearity is a problem that occurs when a

Similarly, for the boundary conditions, we substituteequation (28) into equation (16) to derive the following'ji ¼ 0

on

"

ð32Þ

Equations (30) to (32) comprise the governing equations on the mesoscale. By discretizing equations (30)and (31) using the FEM and GLS method,26 we canderive 'ji (see Appendices 1 and 2). The homogeneous

ð35Þ

1870

Journal of Composite Materials 53(13)number of variables with a high degree of correlation
exist within a prediction model. The task performed
involves the visualization of the correlation coeﬃcient
between variables, which include newly defined dependent variables. This is achieved using the SPM shown in
Figure 3 as well as classifying those variables that havea correlation coeﬃcient exceeding 0.7 as the same variables, thereby reducing them into one.
To prevent multi-collinearity of terms in the prediction models generated with the GP, we use the TOL,
which is one of the indices for determining multicollinearity. The equation used to derive the TOL isshown below.
TOLi ¼ 1 ' R2ið36Þ

In equation (36), Ri represents the determinationcoeﬃcient, with the number of degrees-of-freedom
being adjusted when the regression model is set up inaccordance with variable i. As described in Figure 4,
to prevent multicollinearity, we derive a prediction
model with a TOL for all variables equal to 0.1 orhigher by repeating the TOL evaluation of all variables
and removing the worst TOL variables. In this manner,
we derive the variables that significantly aﬀect theobjective functions, and they are used as references
for the design variables for re-optimization.

Application to non-crimp fabricModeling
We applied to NCF the extraction method of the eﬀective design variables proposed in this paper. We fabricated an NCF model with a mesostructure, as shown in
Figure 5, and performed a multi-objective optimization
using NSGA-2. The fiber bundle is modeled by orthotropic elastic materials, whereas the resin is an isotropic
elastic material. The properties of these materials
are listed in Table 1. The objective functions consistedof the four coeﬃcients, namely the elastic modulus andpermeability coeﬃcient in the in-plane and out-of-plane
directions. The four objective functions are optimized

Figure 5. Design variables of the unit cell of the non-crimp
fabric operating at multi-objective genetic algorithm (MOGA0.
Figure 3. Scatter plot matrix showing the value of the correlation coefficient.

Table 1. Material properties of fiber bundle and resin.
Fiber bundlea(Vf 70%)Make predic!on model by GP

Young’s modulus
Shear modulus

Evaluate TOLof predic!on model variables

Poisson’s ratio
Remove worst TOL variables

Worst TOL > 0.1

N

Y
END

Figure 4. Scheme used to eliminate the multi-correlation in theprediction model.

Resin

DensityYoung’s modulus
Poisson’s ratio
Density

E1 (GPa)
E2, E3 (GPa)G12, G13 (GPa)G23 (GPa))12, )13)23* (g/cm3)
E (GPa))* (g/cm3)

176
32.8
8.55
5.46
0.24
0.19
1.6
3.50
0.37
1.0

a
Subscript 1 denotes the fiber bundle direction, whereas subscripts 2 and3 denote the transverse direction.

Matsuzaki et al.
simultaneously by obtaining the Pareto front, which isthe non-inferior solution group. Note that thePareto front cannot be obtained by the single objectivefunction comprising all parameters with some weights.
There were three design variables, as shown in Figure 5,
namely the fiber-bundling width (x1), fiber-bundle
thickness (x2), and fiber-bundle interval in the crossdirection within a unit cell (x3). The unit cell repeats
periodically in three directions in this study. It should
be noted that directions 1 and 2 (i.e. in-plane) areequivalent for the unit cell.

Generation and extraction of effective
design variables
The NCF model used in this study classifies the regionsof resins from the upper surface, as shown in Figure 6.
We then newly define the proportion of the volume ofresin in each region with respect to the volume of theunit cell as dependent variables, Cr, Cba, Cbb, Hoa, Hob,and Hoc.
To reduce the sum of the original design variables
normalized by the length of the unit cell, x~ 1 , x~ 2 , and x~ 3 ,and the six above-mentioned newly defined dependent
variables, we obtain their correlation coeﬃcients using
the SPM. Based on the results, x1, Cr, Hob, and Hoc areclosely correlated and reduced to x1. Meanwhile, x3 andHoa are reduced to x3, which results in the formulationof a prediction model with five variables, namely x10 ,x20 , x30 , Cba, and Cbb.
We employ the cross-validation method to verify thequality of the fit for the prediction model.27 This is performed by separating all the data into the training data
used for the generation of the prediction model and thetest data that are handled as true values. The leave-oneout cross validation involves the use of respective models
included in the Pareto solution as test data in sequence.
The results obtained from conducting a leave-oneout cross validation on respective objective functions

1871
are shown in Figure 7. The relationship between theestimated values and true values indicates that the determination coeﬃcient is equal to one in Figure 7(a) and(b), while it is 0.94 in Figure 7(c) and (d). In both cases,
the values exceed 0.9, and the model is therefore considered a good-quality prediction model. We derive thefollowing equations after performing the GP.
E1 ¼ 7:0 ) 103 ) x~ 1 ' 8:6 ) 102 ) C~ bb ' 4:8 ) 102 ) C~ ba&'þ 4:0 ) 102 ) x~ 2 ) x~ 1 þ C~ bb þ 1:8 ) 102 ) x~ 2 3 ) C~ baþ 1:0 ) 104

ð37Þ

&
'
2
E3 ¼ 2:5 ) 103 ) x~ 1 ' 1:7 ) 103 ) x~ 1 2 ) x~ 3 ) C~ bbþ 3:8 ) 102 ) x~ 2 þ 2:8 ) 101 ) x~ 2 2 ) ðx~ 1 þ x~ 3 Þþ 5:6 ) 103

ð38Þ

K1 ¼ '2:1 ) 10'3 ) ðx~ 1 Þ þ 2:7 ) 10'3

ð39Þ

K3 ¼ '9:4 ) 10'3 ) ðx~ 1 Þ þ 9:6 ) 10'3

ð40Þ

Equations (37) to (40) are the prediction equations
for the in-plane and out-of-plane elastic modulus andpermeability coeﬃcients. An overall assessment ofequations (37) to (40) reveals that x~ 1 is included in allof the prediction models, and the coeﬃcients are larger,
significantly aﬀecting the objective functions. A carefulassessment of x~ 1 also reveals that the relationshipbetween the elastic modulus and permeability coeﬃcients is a trade-oﬀ relationship.
Because x1 is a variable that is highly correlated with
Cr, Hob, and Hoc, we derive a design knowledge indicating that, in addition to the existing design variables,
increasing the design variables with respect to the

Figure 6. Defined dependent variables obtained by dividing the optimization model into regions.

1872

Journal of Composite Materials 53(13)

16000

12000

test points

Estimated value [MPa](b)Estimated value [MPa]

(a)

12000
16000
True value [MPa]

0.002

0.001
0.001

7000
test points
6000
6000

7000
8000
True value [MPa]

(d)

0.003

test points
0.002
True value [m 2]0.003

Estimated value [m 2]

Estimated value [m 2]

(c)

8000

0.01

0.005
test points
0
0

0.005
0.01
True value [m 2]Figure 7. Validation of the genetic programing (GP) prediction model derived using the proposed method. (a) In-plane Young’s
modulus, (b) Out-of-plane Young’s modulus, (c) In-plane permeability and (d) Out-of-plane permeability.

Figure 8. Design variables of the unit cell of the non-crimp fabric applied to the effective design variables operating at the multiobjective genetic algorithm (MOGA).

regions of Cr, Hob, and Hoc may realize an improved
fiber material mesostructure.

Re-optimization
Reflecting the design knowledge of the NCF model
derived in the preceding section, we increase thenumber of design variables in the regions Cr, Hob, andHoc, as the fiber-bundle width (x1), fiber-bundle thickness (x2), fiber-bundle width interval (x3), change rate of
the fiber-bundle width (x4), and fiber-bundle oﬀset (x5),as shown in Figure 8. It should be noted that the fiberbundle width change rate and fiber-bundle oﬀset can becontrolled by stitching fiber bundles together28 in actual
applications. The cross-sectional areas of the fiber bundles and fiber-volume fraction of the fiber bundle areassumed to be fixed. Note that the fiber volume fractionin the unit cell is not constant.
Figure 9 shows the Pareto solutions of in-plane andout-of-plane properties obtained by optimization
before and after application of the eﬀective design variables. The area marked by the dashed line indicates that
the Pareto solutions improve upon the application ofeﬀective design variables. From Figure 9, we confirm
that the solutions of stiﬀness and permeability multiobjective optimization improve in the overall Pareto
front by applying eﬀective design variables.

Discussion
The optimal shapes that maximize each objective function are shown in Figure 10. The maximum values forMatsuzaki et al.

1873

(a)

(b)

10 )
2

3
2

Especially
improved area

1

BeforeA"er8
10 12 14 16 18
In-plane Young's modulus [GPa]Especially
improved area

Maximize

5

Before
A"er
0
0

〜

0
0

10-3)
10

Maximize〜

In-plane Permeability [m ](Out-of-plane Permeability [m2]

(

-3

5
6
7
8
9
Out- of-plane Young's modulus [GPa]Figure 9. Comparison of the Pareto front of each objective function before and after application of the effective design variables.

Figure 10. Comparison of the maximum values of each objective function before and after application of the effective design
variables. (a) In-plane Young’s modulus, (b) Out-of-plane Young’s modulus, (c) In-plane permeability and (d) Out-of-plane permeability.

the respective objective functions and the optimal
dimensions for each objective function are shown in
Table 2. From Figure 10 and Table 2, we observethat the change rate of the fiber-bundle width (x4) exhibits an interdependent relationship between the elastic
modulus and the permeability coeﬃcients; the maximum modulus has maximum x4, while the maximum
permeability has minimum x4. This implies that thefiber-bundle width (x4) is attributed to the trade-oﬀ
for the elastic modulus and the permeability coeﬃcients. This is because, as x4 increases, the transverse
fiber bundles spread in the in-plane direction, which
increases the elastic modulus. On the other hand, theflow path between the fiber bundles decreases, which
results in a low permeability.
Moreover, by comparing Figure 10(c) and (d), thefiber-bundle oﬀset (x5) exhibits an interdependent

relationship between the in-plane and out-of-plane permeability coeﬃcients; the maximum in-plane permeability has maximum x5, while the maximum out-ofplane permeability has minimum x5. This implies thatthe fiber-bundle oﬀset (x5) is attributed to the trade-oﬀ
for in-plane and out-of-plane permeability coeﬃcients.
This is because, as x5 increases, the main (i.e. largest)flow path between the fiber bundles located in the adjacent unit cell increases and the in-plane permeability
increases, while the projected flow path in the out-ofplane direction decreases. Therefore, the out-of-plane
permeability decreases.
As a detailed investigation of the optimized dimensions of NCF, we employed optimized solutions before
and after application of eﬀective design variables as an
example in Figure 11. Comparing models (a) and (b)in Figure 11, although model (b) has a 2.5% smaller

1874

Journal of Composite Materials 53(13)Table 2. Design variables and objective functions of Pareto solutions having the maximum value of each objective function (numbers
marked with yellow are maximum objective function values).
Design variables

In-plane Young’s modulus
Out-of-plane Young’s modulus
In-plane permeability
Out-of-plane permeability(a)Objective functions

x1

x2

x3

x4

x5

E1(GPa)E3(GPa)K1(m2)

K3(m2)

0.60
0.60
0.30
0.30

0.10
0.10
0.10
0.10

0.30
0.30
0.30
0.30

1.03
1.03
1.00
1.00

0.03
0.03
0.03
0.00

17.9
17.9
10.2
10.2

8.59
8.59
6.01
6.00

7.70 ) 10'4
7.70 ) 10'4
3.18 ) 10'3
3.07 ) 10'3

1.24 ) 10'3
1.24 ) 10'3
1.15 ) 10'2
1.16 ) 10'2(b)

Figure 11. In-plane and out-of-plane permeability inside the unit cell of the non-crimp fabric obtained by homogenization method.
(a) Before application of the effective design variables and (b) After application of the effective design variables.

Matsuzaki et al.

1875

Table 3. Correlation coefficient between fiber bundle widthand each objective function.
Fiber-bundlewidth
at the end ofunit cell (x1)In-plane Young’s modulus
Out-of-plane
Young’s modulus
In-plane permeability
Out-of-plane permeability

Fiber-bundle
width atcrimp(x1 ) x4)

0.98
0.98

0.98
0.99

'0.99
'0.97

'0.98
'0.96

fiber-bundle width at the end of unit cell (x1) and a
1.1% lower fiber volume fraction in the unit cell (Vf),model (b) has a 2.3% and 4.4% higher in-plane andout-of-plane permeability, respectively, while maintaining the elastic modulus. This is because the fiber-bundle
width at crimp (x1)x4) has a strong influence on
the elastic modulus (Table 3), and it is almost thesame value between the two models.
In addition, model (b) has a higher permeability thanthat of model (a) because it has a 2.5% smaller fiberbundle width at the end of unit cell (x1), which adversely
aﬀects the permeability (Table 3). Moreover, the smaller
x1 in model (b) means the expansion of the high permeability region. From those facts, by adding the parameter of change rate of the fiber-bundle width (x4 > 1), thefiber-bundle width at the end of unit cell (x1) can besuppressed while maintaining the fiber-bundle width at
the crimp (x1)x4), which enables attainment of the highmodulus and permeability.

Conclusions
We proposed an optimization method for the stiﬀness
and permeability of CFRP mesostructures by using ahomogenization method. The method generates eﬀective
design variables from optimization results derived froma smaller number of design variables. Furthermore, therelationships between the objective functions and design
variables are formulated into mathematical expressions
that are based on statistical parameters. We examinedthe eﬀectiveness of the proposed method by considering
the optimization problems of the mesostructures of
NCF. In this application, we newly derived the eﬀective
design variables as the change rate of the fiber-bundle
width and fiber-bundle oﬀset in the NCF model in addition to the initial variables. By reflecting the design
knowledge, we conducted another optimization.
The results confirmed that the fiber material structure attained higher elastic and permeability coeﬃcientscompared to the case where the design knowledge was
not applied. In addition, from the investigation of optimization results, we also clarified the trade-oﬀ

relationship between the in-plane and out-of-plane elastic modulus and permeability, as well as the eﬀects ofdimensions of NCF mesostructures on the elastic
modulus and permeability. In this study, we set thein-plane and out-of-plane stiﬀness and permeability
along the axis as objective functions. Consideration ofthese properties along an arbitrary direction as objective functions may expand the performance space;hence, this will be investigated in future work.
Acknowledgements
The authors acknowledge the vitally important encouragement
and support made through the University of WashingtonTohoku University: Academic Open Space (UW-TU:AOS).

Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest withrespect to the research, authorship, and/or publication of thisarticle.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of thisarticle: The Cross-ministerial Strategic Innovation Promotion
Program (Funding agency: JST, Japan).

ORCID iD
Ryosuke Matsuzaki http://orcid.org/0000-0001-8237-1638
Shigeki Yashiro http://orcid.org/0000-0001-5266-2039

References
1. Kouznetsova VG, Geers MGD and Brekelmans WAM.
Multi-scale second-order computational homogenizationof multi-phase materials: a nested finite element solution
strategy. Comput Meth Appl Mech Eng 2004; 193:
5525–5550.
2. Takano N, Zako M, Okazaki T, et al. Microstructurebased evaluation of the influence of woven architecture
on permeability by asymptotic homogenization theory.
Compos Sci Technol 2002; 62: 1347–1356.
3. Song YS and Youn JR. Asymptotic expansion homogenization of permeability tensor for plain woven fabrics.
Compos A Appl Sci Manufact 2006; 37: 2080–2087.
4. Terada K, Ito T and Kikuchi N. Characterization of themechanical behaviors of solid-fluid mixture by the homogenization method. Comput Meth Appl Mech Eng 1998;
153: 223–257.
5. Bendsøe MP and Kikuchi N. Generating optimal topologies in structural design using a homogenization method.
Comput Meth Appl Mech Eng 1988; 71: 197–224.
6. Kobayashi M, Nishiwaki S, Izui K, et al. An innovative
design method for compliant mechanisms combining
structural optimisations and designer creativity. J Eng
Design 2009; 20: 125–154.
7. Yamada T, Nishiwaki S, Izui K, et al. A structural optimization method incorporating level set boundary

1876
expressions based on the concept of the phase field
method. J Environ Eng 2011; 6: 567–578.
8. Wang MY, Wang X and Guo D. A level set method forstructural topology optimization. Comput Meth Appl
Mech Eng 2003; 192: 227–246.
9. Guest JK and Prévost JH. Optimizing multifunctional
materials: design of microstructures for maximized stiffness and fluid permeability. Int J Solid Struct 2006; 43:7028–7047.
10. Tatsukawa T, Nonomura T, Oyama A, et al. Aerodynamic
design exploration for reusable launch vehicle using genetic
algorithm with Navier Stokes solver. Transact Jpn Soc
Aeronaut Space Sci Aerosp Technol Jpn 2012; 10:
Pe_55–Pe_63.
11. Jeong S, Murayama M and Yamamoto K. Efficient optimization design method using Kriging model. J Aircraft
2005; 42: 413–420.
12. Kosugi Y, Oyama A, Fujii K, et al. Multidisciplinary andmulti-objective design exploration methodology for conceptual design of a hybrid rocket. In: Infotech@Aerospace conferences, St. Louis, Missouri, US, 2011,
AIAA. DOI: 10.2514/6.2011-1634.
13. Tatsukawa T, Oyama A and Fujii K. Extraction of design
information from pareto-optimal solutions using genetic
programming: a first report. In: International workshop on
future of CFD and aerospace sciences, Kobe, Japan, 2012.
14. Obayashi S, Jeong S and Chiba K. Multi-objective design
exploration for aerodynamic configurations. In: 35th
AIAA fluid dynamics conference and exhibit, American
Institute of Aeronautics and Astronautics, 2005, AIAA.
DOI: 10.2514/6.2005-4666.
15. Deb K, Agrawal S, Pratap A, et al. A fast elitist nondominated sorting genetic algorithm for multi-objective
optimization: NSGA-II. In: Schoenauer M, et al. (eds)
Parallel Problem Solving from Nature PPSN VI. PPSN
2000. Lecture Notes in Computer Science. Vol 1917.
Berlin, Heidelberg: Springer, 2000, pp.849–858. DOI:
10.1007/3-540-45356-3_83.
16. Akaike H. A new look at the statistical model identification. IEEE Transact Automat Control 1974; 19: 716–723.
17. Berk KN. Tolerance and condition in regression computations. J Am Stat Assoc 1977; 72: 863–866.
18. Affdl JCH and Kardos JL. The Halpin-Tsai equations: areview. Polym Eng Sci 1976; 16: 344–352.
19. Gebart BR. Permeability of unidirectional reinforcements
for RTM. J Compos Mater 1992; 26: 1100–1133.
20. Sanchez-Palencia E. Non-homogeneous media and vibration theory. Lecture Notes in Physics, 127. New York:
Springer Verlag, 1978.
21. Mei C. Method of homogenization applied to dispersion in
porous media. Transp Porous Med 1992; 9: 261–274.
22. Mei CC, Auriault J-L and Ng C-O. Some applications of
the homogenization theory. Adv Appl Mech 1996; 32:
277–348.
23. Ulrich H. Homogenization and porous media. New York:
Springer-Verlag Inc., 1997.
24. Nedanov PB and Advani SG. Numerical computation of
the fiber preform permeability tensor by the homogenization method. Polym Compos 2002; 23: 758–770.

Journal of Composite Materials 53(13)25. Wang JG, Leung CF and Chow YK. Numerical solutions
for flow in porous media. Int J Numer Anal Meth
Geomech 2003; 27: 565–583.
26. Sert C. Lecture Notes for ME 582. Finite element analysis
in thermofluids. In: Cüneyt Sert (ed.) Incompressible flow
solutions. Department of Mechanical Engineering Middle
East Technical University, 2016.
27. Rogers WH and Wagner TJ. A finite sample distributionfree performance bound for local discrimination rules.
Ann Stat 1978; 6: 506–514.
28. Loendersloot R. The structure-permeability relationof textile reinforcements. PhD Thesis, University ofTwente, Enschede, The Netherlands, 2006. Available at:https://research.utwente.nl/en/publications/the-structurepermeability-relation-of-textile-reinforcements.
29. Franca LP and Frey SL. Stabilized finite element methods: II. The incompressible Navier-Stokes equations.
Comput Meth Appl Mech Eng 1992; 99: 209–233.
30. Berdichevsky AL and Cai Z. Preform permeability predictions by self-consistent method and finite element
simulation. Polym Compos 1993; 14: 132–143.
31. Verleye B, Klitz M, Crose R, et al. Predicting
the permeability of textile reinforcements via a hybridnavier-stokes/brinkman solver. In: Proceedings of the
8th international conference on flow processes in composite
materials (FPCM8), Douai, France, 2006.

Appendix 1
Calculation process of characteristic function 'ji
The detailed calculation process of characteristic function 'ji from equations (30) and (31) is described below.
First, equation (30) is incorporated into equation (31)to give'@2 'jl
@2 'ji
@Aj
'
þ
¼ #ij
@yk @yk @yl @yi @yi

ð41Þ

Equations (30) and (41) are transformed into theweak form using the Galerkin method.
Z@'j'wc i d! ¼ 0
ð42Þ
@yi
!
(
Z "
Z
Z
j
@w @'ji @w @'l @wþ'
Aj d! ¼
w#ij d! þ wti d"@yl @yi @yi! @yk @yk!
"
ð43Þ
where ti is given as follows
ti ¼

@'j
@'ji
nk þ l nl ' Aj ni
@yk
@yi

ð44Þ

In equation (43), the area integral term becomes zero
because the weighting function is set to be zero on theboundary. Then, the characteristic functions 'ji and Aj

Matsuzaki et al.

1877

are discretized based on FEM as follows
NEN'
X

'ji ¼

Nm 'jim

m¼1

Aj ¼

NEN'
X )Z

NENA
X

Mm Ajmð45Þ

m¼1

where 'jim is a characteristic function of the flow velocity at node; Ajm is the characteristic function of thepressure at the node; Nm and Mm represent shape functions; and NENj and NENA denote the respective
numbers of 'jim and Ajm nodes over an element.
It should be noted that the weight function is thesame as the shape function as follows
w ¼ Nnwc ¼ Mn

ð46Þ

where n represents the constituent node number of anelement. Equations (45) and (46) are substituted into
equation (42) to derive
NEN'
X )Z
m¼1

!

'Mn

(
@Nm
d! 'jim ¼ 0
@yi

ð47Þ

Equations (45) and (46) are substituted into
Equation (43) to derive
NEN'
X )Z

!

" (
NEN'
X )Z !@Nn @Nm " ( j@Nn @Nmjd! 'im þ
d! 'lm@yl @yi! @yk @yk
!
m¼1
m¼1
(
Z
NENA
X )Z
@Nn
þ
'
Mm d! Ajm ¼
Nn #ij d!
@yi
!
!
m¼1

ð48Þ

Equation (48) can be expressed as equation (49) in
expanded form
NEN'
X )Z

!
" (
@Nn @Nm @Nn @Nm @Nn @Nm
2þþ
d! '11m
@y
@y
@y
@y
@y
@y
1
1
2
2
3
3
!
m¼1
(
NEN'
X )Z @Nn @Nm
þd! '12m! @y2 @y1
m¼1
)
(
Z
NEN'
X
@Nn @Nmþd! '13m! @y3 @y1
m¼1
)
(
Z
NENA
X
@Nn
þ
'
Mm d! A1m
@y1
!
m¼1
Z
¼
Nn d! ði ¼ 1, j ¼ 1Þ!

(@Nn @Nmd! '11m
@y
@y
1
2
!
m¼1
" (
NEN'
X )Z !@Nn @Nm@Nn @Nm @Nn @Nm
þþ2
þ
d! '12m
@y
@y
@y
@y
@y
@y
1
1
2
2
3
3
!
m¼1
(
NEN'
X )Z @Nn @Nm
þd! '13m! @y3 @y2
m¼1
)
(
Z
NENA
X
@Nn
þ
'
Mm d! A1m ¼ 0 ði ¼ 2, j ¼ 1Þ
@y2
!
m¼1
..
.

ð49Þ

Equations (47) and (49) can be expressed in thematrix form
2 * 11 +*K +6 K21
6 * 31 +4 K* 41 +
K* 12 +*K22 +*K32 +*K42 +
K* 13 +*K23 +*K33 +*K43 +
K

2n o3* 14 + 32 , j - 3
F1j'*K24 + 6 , 1j - 7 6 n o 77' 7 6 F2j 7
7*K34 + 76
6 , 2j - 7 ¼ 6
6, - 7
5*K44 + 4 , '3 - 5 4 F3j 5
, 4K
Aj
F

ð50Þ

Kij (i ¼ 1, 2, 3, and j ¼ 1, 2, 3, 4) in equation (50)corresponds to the left side of equation (49) and isexpressed as"
Z !
@Nn @Nm @Nn @Nm @Nn @Nm2þþd!
@y1 @y1@y2 @y2@y3 @y3!
Z@Nn @Nm¼
d!
K12
nm@y @y1
Z! 2@Nn @NmK13d!
nm ¼@y @y1
Z! 3
@Nn
K14
'
Mm d!
nm ¼
@y1
!
Z
@Nn @Nm
K21¼
d! ¼ K12nmmn! @y1 @y2"
Z !
@Nn @Nm@Nn @Nm @Nn @Nm
22
Knm ¼
þ2
þ
d!
@y1 @y1@y2 @y2@y3 @y3Z!
@Nn @NmK23d!
nm ¼! @y3 @y2
Z
@Nn
K24
¼
'
Mm d!
nm
@y2
Z!
@Nn @Nm
K31
d! ¼ K13nm ¼
mn! @y1 @y3
Z
@Nn @Nm
K32
d! ¼ K23
nm ¼
mn! @y2 @y3"
Z !
@Nn @Nm @Nn @Nm@Nn @Nm
K33
¼
þ
þ
2
d!
nm
@y2 @y2@y3 @y3! @y1 @y1
K11
nm ¼

1878

Journal of Composite Materials 53(13)

K34
nm ¼

Z

@Nn
'
Mm d!
@y3
!

ð51Þ

Kij (i ¼ 4 and j ¼ 1, 2, 3, 4) in equation (50)corresponds to the left side of equation (47) and isexpressed as

K43
nm ¼

Z

!

'Mn

Z

Fij (i ¼ 1, 2, 3 and j ¼ 1, 2, 3) in equation (50)corresponds to the right side of equation (49) and isexpressed as
F1j
n ¼
F2j
n ¼
F3j
n ¼

Z
Z

Nn #1j d!
Nn #2j d!

ð53Þ

Nn #3j d!

Because the right side in equation (47) equals zero
F4n ¼ 0ð54ÞK44
nm ¼ 0 in equation (52) means that the global stiﬀness matrix has a zero term on the diagonal, which leads
to instability in the iterative solver of the linear equation.
To avoid this instability, the pressure approximationmust be at least one order lower than the velocity
approximation. This is known as the inf-sup condition.29
To meet the inf-sup condition, the mixed interpolation
function must be used, which results in complex code
and a high computational cost. Therefore, we use the
GLS method to stabilize equation (50) instead of using
the mixed interpolation function.
In GLS, the following stabilization term is added toequations (42) and (43)Z#$+LðwÞ % R 'ji , Aj d!
ð55Þ
!

where L and R are as follows

1 @wc & @2 wiLðwi Þ ¼ 'þ
* @yi * @yj @yj

ð56Þ

!

1 @wc
+
*
@yi
!

ð52Þ

Z

1 le 2
3 4)

ð58Þ

where le is the measure of the element length, and ) is
the kinetic viscosity. The stabilization term that is
added to equation (42) is@Nm
d! ¼ K34
mn
@y3

K44
nm ¼ 0ð57Þ

+ is the stabilization parameter and is usually taken
to be26

+¼

Z

@Nm
'Mn
d! ¼ K14
mn
@y
1
Z!
@Nm
K42
'Mn
d! ¼ K24
nm ¼
mn
@y2
!
K41
nm ¼

#
$
1 @Aj & @2 'ji
þþ #ij
R 'ji , Aj ¼ '* @yi * @yk @yk!
& @2 'ji
1 @Ajþ
' #ij d!
'* @yl @yl * @yi

"

ð59Þ

The stabilization term that is added to equation (43) isZ!
!
"& @2 wi& @2 'ji
1 @Aj
+ '
þ
' #ij d! ð60Þ'
* @yl @yl
* @ym @ym * @yi!

Here, when the second derivative terms in equations(59) and (60) are assumed to be zero, only equation (59)remains, and it is expressed as!
"@wc @Aj @wc'#ij d!
@yi @yi
@yi
!

Z

+

ð61Þ

For simplicity of calculation, * is taken to be one.
Therefore, [K44] and {F4} in equation (50) are as follows
K44
nm ¼
F4n ¼

Z

!
"
@Mn @Mm'+d!
@yk @yk
!

Z

!

'+Mn

!

"
@#ij
d! ¼ 0
@yi

ð62Þ

ð63Þ

Equation (62) indicates that the zero diagonal termof the global system in equation (50) is removed, andthe characteristic function 'ji can be obtained without
the instability problem.

Appendix 2
Validation of homogenized permeability calculation
We validated the homogenized permeability calculationcompared with well-known permeability models:Gebart19 and Berdichevsky.30 Both models used here
assume the use of square-arrangement unidirectional

Matsuzaki et al.

1879

fiber-reinforced material and are expressed as

K== ¼

Gebart model
16
K? ¼ pﬃﬃﬃ
9( 2!rﬃﬃﬃﬃﬃﬃﬃﬃ
"5=2
(
'1
R2
4Vf

8R2 ð1 ' Vf Þ3
57
V2f

ð65Þ

Berdichevsky model
ð64Þ
R2
1 1 ' Vf 2
K? ¼
ln '
Vf 1 þ Vf 2
8Vf

!

#
$#
$
R2
1
K== ¼
ln 2 ' 3 ' Vf 1 ' Vf
8Vf
Vf

Figure 12. Micro-model of a square arrangement unidirectional
fiber-reinforced material used to calculate the homogenized permeability coefficient.

(a)

ð66Þ

!

ð67Þ

where R is the fiber radius, and K? and K== are thepermeability in the transverse and parallel flow to thefiber, respectively. The finite models in the case ofVf ¼ 60% are shown in Figure 12. For simplicity, theresin viscosity is set to one.

(b)

(c)

Figure 13. Validation results of homogenized permeability calculation for square arrangement of unidirectional fiber-reinforced
material. (a) Characteristic flow-velocity profile in transverse flow to the fiber ('11 ) and (b) parallel flow to the fiber ('33 ), where thefiber-volume fraction is 20%. (c) Relationship between fiber-volume fraction and permeability obtained from the Gebart model,
Berdichevsky model, and homogenization method.

1880
The calculation results of the characteristic flow velocity profile are shown in Figure 13(a) and (b) in
the case of Vf ¼ 20%. The relationship between
the fiber-volume fraction in the unit cell and the calculated homogenized permeability is shown in
Figure 13(c). In the transverse flow, the calculation

Journal of Composite Materials 53(13)results agree with the Gebart model. In the parallel
flow, the calculation results agree with theBerdichevsky model for low Vf, whereas they agreewith the Gebart model for high Vf. This tendency is
consistent with that of Verleye et al.31