ANALYSIS OF DRY SPOT FORMATION IN RTM
THROUGH OPTIMIZATION OF

NON-NEWTONIAN RESIN PARAMETERS

A Thesis

Presented to the Faculty of the Graduate School

of Cornell University

in Partial Fulfillment of the Requirements for the Degree of

M.S.

by

Rahul Ghosh

August 2023



© 2023 Rahul Ghosh

ALL RIGHTS RESERVED



ABSTRACT

Resin Transfer Molding, or RTM, is one of the most popular methods in today’s

world for manufacturing composite materials. The capacity for creating flexi-

ble and complex geometries of high volume fraction at lower costs makes it a

very appealing manufacturing process for both the automotive and aerospace

industries. However, no process is without its own drawbacks. One of the ma-

jor issues faced in producing composites through this technique is the inclusion

of voids and air bubbles. Naturally, these act as points of failure or low strength

in the final created part. The first contribution of this thesis deals with identify-

ing void formation mechanisms through simulations of the RTM filling process

in commercial software (ANSYS FLUENT). This is followed by the computa-

tional exploration of resins that exhibit Non-Newtonian behavior and how they

affect void formation. Finally, we use Machine Learning Optimization tools to

determine the best set of resin rheology parameters from the conceivable sam-

ple space, and study the corresponding effects on fill. Through this innovative

approach, we are able to establish a simulation framework for RTM, a unique

fluid-based approach to tackle an industry-wide problem and a deeper under-

standing of how certain material properties can be cultivated to our advantage.



BIOGRAPHICAL SKETCH

Rahul Ghosh is a Master of Science candidate in Mechanical Engineering at the

Sibley School of Mechanical & Aerospace Engineering, Cornell University. He

received a Bachelor of Technology in Mechanical Engineering from the National

Institute of Technology, Durgapur, India, in the year 2020. His field of research

is at the intersection of Fluid Dynamics and Materials Engineering. He works

at the Mechanics for Material Design (MMD) Lab under his PI, Dr. Meredith

Silberstein.

iii



This document is dedicated to my mother, Nilanjana Ghosh, my father,

Anirban Ghosh and my sister, Sanjna Ghosh.

iv



ACKNOWLEDGEMENTS

This body of work would not have been possible without the guidance of my

committee, led by Dr. Meredith Silberstein and Dr. Mahdi Esmaily Moghadam.

I would also like to thank the Air Force Office of Scientific Research (AFOSR) for

their funding and support of this project. I would like to thank the members of

the Mechanics for Material Design (MMD) Lab, especially Dr. Robert Wagner,

for their continued efforts to help me improve all aspects of my research. Finally,

this thesis would not have been possible without the backing of my friends and

family, namely Karan Jha, Ashwin Sahasranaman, Kushan Mitra, Shinjini Roy,

Ashhad Ameer, Kiranmayi Yenduri, Jai Khanna and many others.

v



TABLE OF CONTENTS

Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi

1 Introduction 1
1.1 Polymer Matrix Composites . . . . . . . . . . . . . . . . . . . . . . 1
1.2 PMC Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Contributions of the Thesis . . . . . . . . . . . . . . . . . . . . . . 6

2 Dry Spot Formation and Simulation 10
2.1 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.1.1 Choice of Software . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.2 Governing Equations . . . . . . . . . . . . . . . . . . . . . . 12
2.1.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2 Dry Spot Formation . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1 Multi-scale problem . . . . . . . . . . . . . . . . . . . . . . 28
2.2.2 Causes behind Dry Spot formation . . . . . . . . . . . . . . 33
2.2.3 Replication in FLUENT Simulations . . . . . . . . . . . . . 39

3 Non-Newtonian flow through Porous Media 42
3.1 Non-Newtonian Fluids . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.1.1 Carreau Model . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2 Flow through Porous Media . . . . . . . . . . . . . . . . . . . . . . 46

3.2.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2.2 Porous Media Simulations . . . . . . . . . . . . . . . . . . . 49

4 Statistical Optimization of Parameters 57
4.1 Baseline Model setup . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.1.1 Parameter Sample Space . . . . . . . . . . . . . . . . . . . . 57
4.1.2 Shear Thickening . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.3 Model set up . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.2 Bayesian Optimization . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

5 Conclusions 77

A VOF Model at High Viscosity Ratio 79

B Time Step Convergence Study 80

C Pressure Gradient Comparison in Validation Simulations 81

vi



Bibliography 82

vii



CHAPTER 1

INTRODUCTION

1.1 Polymer Matrix Composites

Polymer reinforced composites are of great importance to the automotive and

aerospace industries. They are strong, light-weight and can conform to a wide

variety of geometries. PMC’s consist of many strong fibers that are surrounded

by a polymer matrix. While the fibers provide strength and stiffness to the

composites, the matrix secures the fibers in place and transfers loads between

them. [1]

Matrix materials vary between thermosets and thermoplastics. Thermosets are

vastly preferred for their superior mechanical properties as well as a higher

melting point. Epoxy resins are most commonly used for this purpose, along

with curing agents or hardeners to enhance curing properties. The epoxy

molecule itself can cross-link with a wide array of molecules to create resins

that display a wide range of viscosities and do not soften on heating.

As the world relies more on the applications of composite materials, there has

been increased demand for low cost, high performance and good quality man-

ufacturing techniques [2]. Additionally, the versatility of polymer composites

presents a large research domain when it comes to the enhancement of material

properties by incorporating established methods from the fields of fluid me-

chanics and chemical engineering.

1



Figure 1.1: A typical structure of a carbon-fiber-reinforced polymer, repro-
duced from [3]

1.2 PMC Processing

Current techniques for industrial manufacturing of composites include prepeg

molding, pultrusion, resin transfer molding, filament winding, compression

molding, and hand lay-up [1]. Resin Transfer Molding has gained popularity

for its ability to produce high-quality geometrically-flexible materials at a com-

paratively low cost. The popularity has been compounded by the introduction

of low viscosity resins that enhance flow properties.

The process involves placing aligned fiber preforms in a mold, which act as

the structural framework. The mold is closed and resin is injected into the mold

cavity, which displaces the air, causing it to escape from strategically placed

vents. There are many variations to this part of the process that often show in-

creased fiber content, such as Vacuum-Assisted RTM (use of vacuum at the exit

vent of the molding tool [4]) or Compression RTM (compression of mold to re-

duce filling time and void content [5]).

The resin is then allowed to cure and the part is removed from the mold. In

some cases, the mold is preheated to accelerate the curing process and reduce

2



the processing time. The parts may also be post-cured in an oven. Curing of

the resin directly affects the mechanical performance and dimensional stability

of the final part, as well as the desired quality. Prior work has been done in

simulating and studying the curing mechanism in RTM [6].

Many efforts have been made to streamline the process and negate the disad-

vantages, in order to successfully replace autoclave processing as the primary

method of choice. Some examples of this include the work done by Huang, Z.,

et al [7] in their investigation of resin flow behavior in a complex mold, through

which they were able to design an automotive wheel rim of carbon fiber pre-

form and epoxy resin. Similarly, the work done by Brocks, T., et al [8] in the

assembly of a laboratory scale RTM system for different carbon fiber fabrics dis-

played how appropriate parameter choices helped achieve the 2% void content

range that is characteristic of the process.

Beckwith, S.W., et al [9] were able to demonstrate the reduction of labor, skill

and volatile emissions in RTM as compared to hand layup processes. Gonza-

lez, R., et al [10] and Poodts, E., et al [11] presented a study on the RTM process

parameters and how they affect the quality of the final part. The impregnation

process itself was studied by Chan, A.W., [12] while Lee, J., et al [13] provided

more information on the mold filling and cure process.

When it comes to fluid properties of resins, the most typical ones used are New-

tonian in behavior. This is often done to eliminate the complications that arise

from Non-Newtonian flow, as we shall see in Chapter 3. As mentioned in Br-

uschke, M.V., et al [14], non-Newtonian effects are often neglected as these ef-

3



fects are not fully understood. Similarly, in Hourng, L. W., et al [15], the change

in resin viscosity during flow is either neglected or modeled as power law.

The market for Resin Transfer Molding has increased multifold in recent years,

in part due to the decreased volatile emissions, making it a sustainable and cost-

effective alternative to many processes. As compared to VARTM and Com-

pressed RTM, Resin Transfer Molding provides great design flexibility, bet-

ter surface finish and performance, large-scale production capabilities and the

framework for implementing various design controls through the use of trans-

ducers and sensors [2]. The aerospace industry has always preferred autoclave

Figure 1.2: A typical RTM process scheme diagram, reproduced from [16]

manufacturing when it comes to composite parts. As far as processes are con-

cerned, it is well-researched and can produce high quality and fiber volume

fraction, while greatly reducing the presence of voids. However, with increas-

ing material costs, the need for fast and reliable methods at lower running costs

has taken precedence. RTM acts as a useful alternative method, being less labor-

intensive and concerning for the environment as compared to wet lay-up and

4



autoclave methods [17]. It also requires much simpler implementations to cre-

ate the aforementioned process variants. However, as is common with all fab-

rication processes, voidage in the final laminate should be as low as possible,

and since RTM (generally) involves no further consolidation, the impregnation

process needs to be controlled with some care. It has also been established by

Judd, N.W.C, et al [18] that increased void content leads to detrimental mechan-

ical properties. They were able to show that even a 1% increase in void content

results in a 7% decrease in short beam shear properties, which makes the reduc-

tion of voids a major priority.

Table 1.1 takes a look at autoclaving versus RTM from a manufacturer’s perspec-

tive [17]. As one can see, autoclaving is preferred when higher product quality

is desired over mass production. This makes it particularly useful in specialized

aerospace fields and automotive racing environments as well. RTM is a low-cost

alternative, especially due to its great high volume production properties. Fur-

ther work also shows how RTM offers around 25% reduction in process time as

well as reduced exposure to the resin itself [17].

There are many existing methods of choice for simulation modeling of RTM. As

far as modeling flow through porous media is concerned, Dave, R., [19] was able

to provide a unified approach for modeling composite processing, through gov-

erning equations derived from Darcy’s Law. Many attempts have been made

to simulate resin flow through a porous media by developing computational

codes. For instance, Martin, G., et al [20] amended an existing 2D finite ele-

ment routine called POLYFLOW, for RTM systems. Coulter, J.P., et al [21] sim-

ulated 2D resin impregnation with the help of a finite difference code named

5



TGIMPG. By solving Darcy’s Law to obtain a pressure profile, Trochu, F., et

al [22] were able to develop a 3D finite element solver for RTM problems called

RTMFLOT. The use of finite volume models by Souza, J.A., et al [23] showed

great promise in multi-layer reinforced media. Multiple efforts have involved

RTM process simulations based on a finite element-modified control volume

(CV/FE) method to investigate resin flow patterns, such as the work done by

Laurenzi, S., et al [24].

All these methods can be experimentally verified by RTM setups, where the

resin flow front distance is measured and compared with the numerical so-

lution results. Naturally, RTM has its own disadvantages, some of which in-

clude higher tooling cost, greater tooling construction, complicated molds and

reinforcement loading for complex parts. However, the primary reason behind

preferential treatment for other processes stems from its trial and error nature,

Racetracking, Void Formation and high cycle time. This body of work provides

multiple ways to overcome these challenges and help establish RTM as the pri-

mary industrial process.

1.3 Contributions of the Thesis

Around 50% of polymer matrix composites are generated for the aerospace in-

dustry. The combination of strength and durability makes them great building

materials for airplanes, spacecrafts, jetplanes, etc. They have been widely used

in various commercial sectors as well, Airbus being one of the notable benefi-

ciaries. Automotive manufacturers have greatly increased their usage of com-

posites as well, in the form of dashboards, car covers, engine covers, car floors

6



Property RTM Autoclave

Production
output

Moderate
(102–105/annum) Low (101–102/annum)

Void
content Low (< 2%) Low (< 2%)

Labour
intensity Low–moderate High

Achievable
Vf Moderate (30–65%) High (50–70%)

Typical
applica-
tions

Low to medium volume,
automotive parts, non-
structural components
(e.g. automotive spoil-
ers), structural items (e.g.
propeller blades, missile
boxes)

Aerospace industry, For-
mula 1 automotive, sport-
ing goods

Table 1.1: A comparison of RTM and autoclaving processes from a manu-
facturer’s perspective, reproduced from [17]

and even seats. Construction remains an important domain as well, whether it

is storage tanks, various house components and even different types of boats.

When it comes to mass production, cost and sustainability are major factors.

RTM assures that both constraints are well satisfied, while providing a compar-

atively high quality of product. When it comes to strength and durability, the

identification of weak points and failure mechanisms become extremely impor-

tant, as well as the consequent reduction in the same. Improving these aspects

will go a long way in boosting the economy, by enhancing the utility of such

materials.

As we shall see in further chapters, the resin injection into the mold pre-form is

essentially governed by Darcy’s flow through a porous media of given isotropic

or anisotropic permeability. In itself, the flow through porous media is a topic of

7



extensive research, both theoretical and experimental. By setting up a validated

simulation framework, we are able to analyze flow patterns and the effects of

various parameters at much reduced costs. If the variation of fill with process

parameters was to be replicated experimentally, it would lead to a large waste

of resources, increased testing costs and would also vastly limit our ability to

study resins with properties different from those commercially available. The

only costs incurred in our case would relate to computing resources, which has

been improved greatly in the last few decades.

In Chapter 2, we discuss the various causes behind void formation. Some of

them include the presence of obstacles in the mold and regions of increased per-

meability that provide less resistance to resin flow. All of these lead to uneven

filling, which is the major challenge to tackle. Conventionally, RTM utilizes

Newtonian resins, or fluids which have the same viscosity irrespective of the

strain rates faced. Our novel approach is to use non-Newtonian resins, or resins

whose viscosities change with strain rate, to see if they are able to overcome the

hurdles causing void formation. In itself, non-Newtonian flow through Porous

Media is an entire research field that is still being understood, with great bene-

fits to oil recovery and chemical engineering. Our work benefits from prevailing

research and poses additional questions for further accomplishments.

Once we have set up our simulation framework and calibrated it to study the

effect of non-Newtonian resins, a wide range of parameters need to be sam-

pled to be able to narrow down the regime which will give us the best fill. This

is where a multi-modal statistical optimization like Bayesian Optimization is

employed, for its ability to traverse through a large sample space and deduce

8



where the best performance is expected based on certain output parameters. It

uses a Gaussian approach to ensure that all regimes are both explored and ex-

ploited in the right balance through an informed search algorithm. Overall, our

work combines materials engineering, fluid mechanics and machine learning

to provide a platform for better understanding and execution of a widely-used

industrial manufacturing process.

9



CHAPTER 2

DRY SPOT FORMATION AND SIMULATION

2.1 Simulation Setup

2.1.1 Choice of Software

As discussed in the previous chapter, a lot of research has gone into develop-

ing numerical computation codes to be able to simulate the position of the resin

flow front interface during fill. The advantage of such methods include the

ability to predict fill time, fill behavior and the effect of changing process pa-

rameters, with relative ease. As the computational codes can easily be validated

by experimental setups, they can be used to further research in the same field.

Many commerical softwares have also been developed to emulate the indus-

trial process, such as Moldex3D, PAM-RTM and Autodesk’s Moldflow. These

have built-in RTM modules and in some cases, incorporate cure evolution and

modeling. Our software of choice is ANSYS FLUENT, which does not have a

built-in RTM module, nor does it incorporate cure evolution models.

The reason behind our choice is that the advantages of using FLUENT out-

weigh the disadvantages. ANSYS FLUENT is mainly a commercial software

for fluid dynamics simulation that is widely used across industries. Resin flow

into a mold is itself a subset of fluid dynamics called multi-phase flow, where

the two fluid phases involved are the injected resin and the escaping air. As far

as cost effectiveness and utility goes, ANSYS is available on a student license

10

https://www.moldex3d.com/
https://www.esi.com.au/software/pamrtm/
https://www.autodesk.com/products/moldflow/overview
https://www.ansys.com/products/fluids/ansys-fluent


that only restricts the maximum number of elements and certain paralleliza-

tion capabilities. FLUENT is incorporated into ANSYS Workbench. As it is a

widely used industrial tool, there is a plethora of guiding information and ex-

pert tutorials on the internet, as well as a wealth of shared knowledge amongst

users. FLUENT comes with the strong capability of parallelization over multi-

ple cores and GPUs, drastically reducing processing time for simulations. The

software also comes coupled with CFD-Post, a post-processing software that

allows us to analyze flow patterns, flow velocity, pressure variation and vari-

ous other state variables at all points and elements. Additionally, most CAD

models generated from third-party software can be directly imported into AN-

SYS Workbench, allowing us to study a wide variety of geometries. Finally,

when it comes to other available software, there is a huge restriction in terms

of command line or Python usage that can invoke the use of the software. A

Python interface is particularly important during the optimization part of our

research, since we rely on built-in Python modules to run statistical algorithms,

for which we need an interface with the software. ANSYS FLUENT has incor-

porated a particularly useful module known as PyFluent that allows us to run

Python code that directly interacts with the Text-User-Interface (TUI) of FLU-

ENT to execute commands. All of this helps us expand the domain of research

by making it readily available for future use and analysis.

There are associated constraints in using ANSYS FLUENT. Certain restrictions

in viscosity and process parameter choices are discussed in Chapter 4. There is

a large amount of fine-tuning involved to avoid convergence issues. With larger

time durations, the software tends to slow down due to accumulation of data.

There is also a lack of appropriate rheology models for shear thickening fluids.

11

https://cfd2012.com/cfd-post.html
https://fluent.docs.pyansys.com/version/stable/


As a commercial software, very little can be done to change the computational

code to suit our various needs. There is a lack of a cure modeling infrastructure

as well, which is an integral part of the RTM process. However, our assump-

tion is that the cure does not proceed significantly before the fill of the mold is

achieved and causes negligible changes to viscosity before fill.

There is an established history of FLUENT usage in RTM simulations as well.

Porto, J.S., et al [16] [25] displayed the use of FLUENT in both conventional

RTM and Light RTM (resin is injected into an empty channel with no porous

medium) processes, as well as the use of FLUENT to simulate the manufactur-

ing of complex parts like propellers for naval propulsion. The work done by

Ribeiro, G.G., et al [26] proved to be a framework for FLUENT usage, as they

were able to incorporate user-defined functions to the software and compared

the numerical solutions with known analytical solutions for RTM processes.

2.1.2 Governing Equations

When it comes to any numerical solution, it is very important to establish the

underlying governing equations which are being solved. Without a proper un-

derstanding of the governing equations, it is difficult to derive any meaning

from the corresponding graphs and results. At its heart, the injection of a sin-

gle phase fluid into a cavity is governed by the fundamental equations of fluid

mechanics- conservation of mass, conservation of momentum and conservation

of energy. The numerical solution of these equations form the basis for the re-

search field of fluid mechanics and heat transfer, and has been studied for cen-

turies. The equations are essentially coupled equations of velocity, pressure and

12



temperature, which are the main state variables when it comes to fluid flow. De-

pending on the flow constraints, changes in density, viscosity, body forces and

surface forces might also play a role.

Two important distinctions are necessary for our process. Since we are dealing

with a multi-phase problem (two fluids co-exist, resin and air), volume frac-

tion of each fluid in the domain is another important parameter and both the

density and the velocity are representative of this. At the same time, since we

are assuming that the curing of the resin occurs after fill, temperature remains

constant throughout our resin fill, discarding the need to solve the conserva-

tion of energy equation. The equations for continuity and momentum are often

grouped together and known as the Navier-Stokes equation for fluid flow. The

differential form of the equations boil down to [27]:

∂ρ

∂t
+
∂(ρui)
∂xi

= 0 (2.1)

∂(ρui)
∂t
+
∂(ρuiu j)
∂x j

= −
∂p
∂xi
+
∂τi j

∂x j
+ ρ fi (2.2)

where ρ represents density, t represents time, x represents the spatial domain,

p is the pressure and fi represents accelerations due to body forces. The sub-

indices of i and j represent the respective variables in the corresponding di-

mension (i.e. 1,2,3 or x,y,z). For instance, in a 3D domain, the conservation of

momentum would give us three separate equations. τ represents the deviatoric

stress tensor, which is given by:

τi j = µ

(
∂ui

∂x j
+
∂u j

∂xi

)
(2.3)

13



where µ represents the Newtonian viscosity. Now, for flow through porous me-

dia, Darcy’s Law takes effect, which is given by:

v⃗ = −
K⃗
µ
∇P (2.4)

Here, v⃗ is the fluid velocity vector, K⃗ is the permeability tensor of the porous me-

dia and ∇P is the pressure gradient in the respective directions. For anisotropic

porous media, the various components of the permeability tensor in the respec-

tive directions play a major role [28]. Permeability is an intrinsic property of all

materials and refers to how connected pore spaces are to one another. Higher

permeability signifies lesser resistance to fluid transmission. According to FLU-

ENT Documentation, when it comes to incorporating Darcy’s Law in our mo-

mentum equations, a source term F⃗ is added to the right-hand side of Equation

2.2, which is used to model the resistance imposed to the flow by the fibrous

reinforcement (porous media):

F⃗ = ∇P = −
µ

K⃗
v⃗ (2.5)

Porosity is another material property of importance in porous media. It is the

measure of the vacant or empty space in the material. For fibrous media, it is

directly related to the fiber volume fraction as [28]:

ϕ = 1 − v f (2.6)

where ϕ is the porosity and v f is the fiber volume fraction. The permeability

of any material directly relates to the porosity, as shown by Carman, P.C. [29].

14

https://www.afs.enea.it/project/neptunius/docs/fluent/html/ug/main_pre.htm
https://www.afs.enea.it/project/neptunius/docs/fluent/html/ug/main_pre.htm


The velocity vector represented by Darcy’s law is the apparent velocity or Darcy

velocity, which is the product of the actual fluid velocity v⃗ f and ϕ, the porosity

of the fiber reinforcement:

v⃗ = v⃗D = ϕv⃗ f (2.7)

All these equations hold true for single phase fluid flow through porous media.

When it comes to multi-phase flow, another state variable becomes incredibly

important, which is the volume fraction of the primary phase in the cell domain

f. The volume fraction of the secondary phase in the cell domain automatically

becomes 1-f. In all numerical solutions, it is very important that the volume

fraction of each phase is conserved. The continuity of the advective transport of

volume fraction is represented by the following differential equation [16]:

∂ f
∂t
+
∂( f ui)
∂xi

= 0 (2.8)

Hence, this equation is solved in parallel with Equation 2.1 and modified Equa-

tion 2.2. The phase properties such as density and viscosity are represented by

phase-averaged equations, which account for the respective volume fractions of

each phase in the mixture. For a resin-air system like ours, the phase-averaged

equations for density and viscosity are given by [25]:

ρ = fρresin + (1 − f )ρair

µ = fµresin + (1 − f )µair (2.9)

Overall, these set of equations govern the multi-phase flow of resin and air

15



in porous media, while neglecting temperature effect on resin flow. These equa-

tions are further simplified by assuming low Reynolds number of flow and in-

compressible, isotropic Newtonian fluids. Surface tension and air pressure ef-

fects on resin front are also neglected, which is discussed in detail when we talk

about void formation later in this Chapter.

Boundary Conditions

As mentioned earlier, the continuity of the phase volume fraction is extremely

important in numerical solutions. For resin injection into a fibrous reinforce-

ment, the accurate modeling of the flow front interface is crucial to the respec-

tive numerical simulations. As the flow front enters the cell domain, Equation

2.8 is obeyed. At the same time, since it is a free surface, the free-boundary

condition applies [30]:

σ · n⃗ = 0 (2.10)

where n⃗ is the unit normal vector to the corresponding surface and σ stands for

the stress tensor. The fibers are assumed to be rigid, with perfectly fixed location

and the air influence is neglected. The mold walls are assumed to be stationary

during resin injection and the fluids obey no-slip conditions, given by [30]:

v⃗ = 0 (2.11)

At the resin injection gates, either a prescribed velocity of injection or an applied

pressure gradient starts the filling of the mold, given by some v⃗ = v⃗0 or p = p0.

16



2.1.3 Methodology

Volume of Fluid method

We have already defined the governing equations for the flow front interface

of the primary phase in multiphase flow. It is imperative that our numerical

method is able to track the movement of the interface accurately, which helps

us validate our simulations and predict the locations of void formation. The

Volume of Fluid (VOF) method is one such numerical technique, developed by

Hirt and Nichols in 1981 [31].

It follows a Eulerian approach that can track and locate the free surface in a

two-phase flow. The main principle follows the use of a phase fraction α that

varies between 0 and 1. A value of unity signifies that the location is entirely

occupied by the primary phase, and a value of zero indicates that the location is

completely devoid of it. An average value of αwould signify fractional occupa-

tion of the primary phase in the location. If any location has an average value

between unity and zero, it means that an interface exists in that region. The sum

of the fractional volumes for all phases is equal to unity [32].

The VOF method consists of three ingredients: a scheme to locate the surface,

an algorithm to track the surface as a sharp interface moving through a compu-

tational grid, and a means of applying boundary conditions at the surface. The

advection equation (Equation 2.8) governs the fraction of volume of the primary

phase as the interface moves in time and space. The advection equation makes

the interface discontinuous and therefore, reconstruction of the calculated in-

terface is required to obtain a smooth curve for curvature, as well as normal

17



estimations at each interface location. The interface is typically smeared over a

few cells owing to the change in volume fraction as it goes from one phase to

the other.

The issues of divergence in solution are encountered when there is a sudden

jump or deviation in the fluid properties. For instance, miscible fluids would

create substantial problems when it comes to interface tracking. Accurate rela-

tions between fluid properties at the interface are absolutely vital for modeling.

At the same time, large differences in viscosities of the two fluids cause similar

problems in numerical prediction, as shown in [33]. In Appendix A, we have

been able to show the same discrepancies by choosing rather large viscosities

for our resin as opposed to air, and then compared the same simulation with a

secondary phase of comparable viscosity. The failure of the former sums up the

main drawbacks of using the VOF solver.

Finite Volume Method

We discussed various numerical methods in Chapter 1 that have shown promis-

ing results for simulating RTM. ANSYS FLUENT is a commercial software that

utilizes the Finite Volume Method (FVM) to solve problems of fluid flow and

heat transfer. As compared to the Finite Difference Method, FVM is applicable

for non-uniform geometries and is conservative in nature. As the conservation

laws themselves act as the governing equations in fluid mechanics, it makes

sense to utilize an approach that uses the integral form of the conservation equa-

18



Figure 2.1: A reprsentative domain showcasing the Volume of Fluid (VOF)
approach for two immiscible fluids, reproduced from [32]

tion as the starting point [27].

Essentially, the domain is split into a finite number of control volumes, and

the conservation laws are applied on each control volume in the integral form.

The volume and surface integrals are suitably approximated using numerical

laws and quadrature. Boundary conditions are then implemented to obtain an

algebraic set of equations, which could be linear or non-linear in nature. Vari-

ous methods are prevalent for solving these set of equations, such as iterative

methods, which involve linearizing around guess values, solving iteratively and

updating guess values till we attain convergence (or fall below a specified tol-

erance threshold). The flow-chart shown in Figure 2.2 signifies the process for

transient multi-phase flows, such as ours. Enthalpy and Species Conservations

do not apply in our simulations.

19



The transient formulation of the solver is fully implicit, and hence, no sta-

bility criterion has to be met. However, due to the complexity of the multiphase

coupling, the time step should be limited to ensure sufficient computational ac-

curacy [34]. In Appendix B, the optimal time step has been obtained empirically

as 0.1 s, through time-step convergence studies that show the trade-off between

solution accuracy and simulation run-time. For additional stability and accu-

racy, we choose the Bounded Second Order Implicit transient formulation.

The flow is laminar for all cases, with sufficiently low Reynolds numbers. The

porous media is configured by altering Cell Zone Conditions to Porous Me-

dia Formulation, utilizing Superficial Velocity. The porosity is registered ac-

cordingly, and the permeability is input as viscous resistance in each direction,

which is simply the inverse of absolute permeability in the respective direction.

For Pressure-Velocity Coupling, we choose SIMPLE (Semi-Implicit Method for

Pressure Linked Equations), which utilizes pressure correction and shows good

results for transient simulations. Spatial Discretization for the Gradient is Least

Squares Cell Based, the staggered PRESTO for Pressure, Second Order Up-

wind for Momentum and Compressive for Volume Fraction. As shown in [34],

around 50 iterations are sufficient to decrease the normalized residuals to a

value below the convergence limit. We have chosen a lower convergence thresh-

old of 10−6 for all residuals to obtain accurate solutions.

20



Figure 2.2: A flow-chart of the multi-phase solution procedure in ANSYS
Fluent, reproduced from [34]

21



2.1.4 Validation

Contour Verification

A large portion of the work in this section was contributed by Shreyasvi Gowda

(B.S., Mechanical Engineering, Cornell University ’21). For our initial process

and simulation verification, we take a look at the work done by Isoldi, L. A.,

et al [35] in the three-dimensional numerical modeling of RTM processes using

FLUENT. Their verification included the use of a mold cavity in the shape of a

rectangular box of dimensions 0.4 × 0.3 × 0.1 m with inlet and outlet nozzles of

8 mm diameter. A resin of density 916 kg/m3 and Newtonian viscosity 0.07115

Pa s is injected into a fibrous reinforcement of permeability K = 3.89 × 10−9m2

and porosity ϕ = 0.88. An injection pressure of 70 kPa is applied at the inlet.

We are able to compare the respective fill contours in FLUENT and PAM-RTM

Figure 2.3: FLUENT Fill contours at t = 81 s and t = 181 s respectively,
image on left is from [35] and image on right is our simulation

22



(which has a built-in RTM module) at different time durations as a means of

software verification. As shown in Figure 2.3, we are able to match the FLUENT

contours at those times. Figure 2.1 shows how a cell completely filled with resin

is signified in red, a cell completely devoid of resin (or filled with air) is signified

in blue, and the interface is displayed by the respective colors in between the

spectrum. While contour tracking reinforces our confidence in the use of the

VOF method in FLUENT to track the flow-front interface, we need to validate

the same with experimental data to fully validate our simulations.

Experimental Validation

For experimental validation, we take a look at the work done by Oliveira, I.R.,

et al [36] [37]. They used particle-filled resins for the manufacturing of the com-

posite part, which presents various advantages such as cost reduction, better

mechanical properties or improvements to flame resistance. However, the in-

crease in viscosity upon adding particles contributes to a slower flow of resin,

which is often responsible for defects such as void formation, poor saturation of

fiber tows and longer production cycle time.

The experiments were performed by injecting orthophthalic polyester resin

filled with calcium carbonate particles (CaCO3) of various concentrations (40%

and 30%) into a fiberglass mat. These were performed in a stainless steel RTM

mold with cavity dimensions of 320 × 150 × 3.6 mm with the inlet and outlet

points (vents) shown in Figure 2.4. There is a camera positioned above the mold

for monitoring the flow-front advance with the aid of a timer. The rectangular

region of the mold is specially designed to facilitate the formation of a linear

23



Figure 2.4: The experimental mold set up and preform insertion, repro-
duced from [37]

profile in resin flow front advancement.

For the 40% CaCO3 filled resin of density ρ = 1430 kg/m3, the viscosity is given

as µ = 2113 cP and the corresponding permeability and porosity is found to be

K = 1 × 10−10m2 and ϕ = 0.7 respectively. It is worth noting that the permeability

is found through experimental measurements of 1D rectilinear flow.

For the 30% CaCO3 filled resin of density ρ = 1380 kg/m3, the viscosity is given

as µ = 1352.53 cP and the corresponding permeability and porosity is found to

be K = 5.5 × 10−10m2 and ϕ = 0.76.

The injection pressure is ramped from 0 to 21200 Pa for around 100 seconds,

at which point it remains constant for the rest of the fill. For the simplicity

of our simulation, we have applied a constant pressure difference of 21200 Pa

across the inlet and outlet throughout our simulation. With these parameters

in mind, we are able to set up a 2D simulation of a 320 × 150 mm rectangular

section and compare the various contours with the respective experimental im-

ages in Oliveira, I.R., et al [36]. We are also able to compare our contours with

those of the 2D PAM-RTM simulations in the same body of work. As we can

see in Figures 2.5 and 2.6, the flow front advancement in our contours closely

24



Figure 2.5: Our FLUENT simulation contours based on the specified pro-
cess parameters in [36].

Figure 2.6: (a) The corresponding PAM-RTM simulation contours and (b)
The experimental flow front position for the 40% CaCO3 filled
resin at different times, as reproduced from [36].

matches those in the experimental images in terms of distance travelled. When

it comes to the shape of the flow front, it can be attributed to the 2D nature of

our simulations. As mentioned in Oliveira, I.R., et al [36], the region close to the

injection port assumes a ring or radial shape form in the numerical simulation,

as compared to the experiment in which the injection is performed through the

bottom of the mold. The flow front would become rectilinear when the pres-

sure gradient becomes linear. We are able to confirm the same by comparing

the pressure gradient profiles in Appendix C at a particular time duration.

25



In the case of the 30% CaCO3 filled resin, we can compare our flow front ad-

vancement with experimentally recorded data [37]. We have chosen the fur-

thermost iso-surface where resin volume fraction is its average value (0.5) as

a marker for the flow front position in our simulations. As we can see in the

contour comparisons in Figures 2.7 and 2.8, our flow front advancement closely

matches experimental images in distance travelled and matches the expected

contour shape in the 2D PAM-RTM simulations.

As we can see in Figure 2.9 (b), the FLUENT simulation flow front advance-

ment closely matches the flow front position from both experimental data and

PAM-RTM simulations. The maximum error between FLUENT data points and

experimental data (after initial 100 seconds) is found to be under 10%. The error

can be easily explained by the 2D nature of our simulations and the difference in

our injection pressure profile. As corroborated by Oliveira, I.R., et al [37], there

is also an agglomeration of CaCO3 particles at the beginning of fill, which makes

flow front progression difficult. This is evident from how the simulation flow

front travels faster than the experimental flow front for both PAM-RTM and

FLUENT simulations in the beginning. Hence, the FLUENT model is able to

represent RTM and the associated physics really well, and can help us simulate

resin fill and predict fluid flow profiles in the mold.

26



Figure 2.7: Our FLUENT simulation contours based on the specified pro-
cess parameters in [37]

Figure 2.8: (a) The corresponding PAM-RTM simulation contours and (b)
The experimental flow front position for the 30% CaCO3 filled
resin at different times, as reproduced from [37].

27



Figure 2.9: (a) The 30% CaCO3 filled Resin front position vs. time, repro-
duced from [37]; and (b) The extracted data points from [37]
compared with our FLUENT simulation data points

2.2 Dry Spot Formation

2.2.1 Multi-scale problem

One of the key aspects of RTM involves the passage of resin through a fibrous

reinforcement or preform, which is considered as a porous medium. Depending

on the type of preform used, they can be composed of either natural or artifi-

cial fibers, which could be arranged in tows and yarns. Each preform consists

of interlaced yarns in different patterns, and each yarn consists of thousands

of solid filaments. These are formed by weaving, knitting, braiding, crochet-

ing, knotting or pressing fibers together [2]. Their respective arrangements,

alignments and composition can create a wide array of fibrous reinforcements,

each of which would behave differently with regards to resin injection. Long,

A.C., [38] defined a hierarchy for fiber arrangement, which included fibers at

the microscopic scale, repeated yarns at the mesoscopic scale and fibers on the

macroscopic scale. Each of these are characterized by respective length scales,

28



which makes the filling of resin a multi-scale problem.

Void formation essentially becomes a concern from a multi-scale mechanics

standpoint. Air bubbles could get trapped inside fiber bundles, in between in-

terlaced fiber tows and even in various larger areas of the mold. These are re-

spectively categorized as micro scale (in fiber tows), meso scale (between fiber

tows) and macro scale voids. In Table 2.1, we see the different length scales of

voids commonly found in composites. These are especially pronounced when

resin flow itself is non-uniform. It is important to note that the inter-bundle and

intra-bundle resin velocity determines whether meso or micro scale voids will

form [39] [40].

Observation
scale

Size of com-
putational
domain

Flow characteristics Medium

Macro scale >0.1 m Darcy’s flow Porous
medium

Meso scale 1 - 10 mm Stokes’ flow, Darcy’s flow
with capillary pressure

Open gap,
Porous
medium

Micro scale 10um to
1mm

Stokes’ flow with surface
tension

Solid fila-
ments and
tiny pore

Table 2.1: The different scales associated with a fibrous reinforcement, and
the corresponding flow characteristics, reproduced from [41]

Macroscale voids are often referred to as dry spots, or locations in the mold

where resin has not impregnated. At their length scale, they are visible to the

naked eye and hence, pose a major threat to mass production. The motion of

resin flow front on this scale is governed by Darcy’s flow through porous me-

dia (Equation 2.4). It is worth noting that in most cases, the permeability of the

preform is anisotropic. Variation in permeability in different directions will lead

29



to non-uniform flow of resin, and ultimately becomes one of the major causes

behind dry spot formation. Higher fiber volume fraction leads to lower perme-

ability and higher resistance to resin flow. It is also worth noting that surface

tension is neglected when it comes to modeling at this length scale, since sur-

face tension forces are inversely proportional to the characteristic length. This

is corroborated by Park, C.H., et al [42].

Mesoscale voids are formed in between yarns, where open gaps are encoun-

tered. In these gaps, flow is represented by Stokes equation and viscous force is

dominant. However, when it comes to the interface of resin and the air inside

the yarns, Darcy’s law applies once again, with special emphasis on capillary

pressures. Mesoscale voids are caused by the differences between a potentially

high value of intra-bundle resin velocity as compared to the inter-bundle resin

velocity. In some literature, mesoscale voids are referred to as macro voids or

macro pores, which leads to a considerable bit of confusion when discussing the

scale of this problem.

Microscale voids are formed inside fiber yarns, in between solid filaments. They

are caused by the differences between a potentially high value of inter-bundle

resin velocity as compared to the intra-bundle resin velocity. At such a small

length scale, surface tension plays a major role. The flow is modeled through

Stokes equations with a special emphasis on surface tension.

The formation of all three voids at their various length scales has been shown

by Park, C.H., et al [42], as can be seen in Figure 2.10. The formation of meso and

micro scale voids is dependent on the differences between viscous and capillary

flow [43]. If capillary forces between fibers dominate, meso-voids form. If vis-

30



Figure 2.10: The multi-scale void defects, shown as : (a) dry spot (race-
tracking effect), (b) spherical meso-voids between fiber tows,
and (c) cylindrical micro-voids within fiber tows, all repro-
duced from [42]

cous flow dominates, micro-voids form. Hence, this has been aptly quantified

by a dimensionless velocity known as capillary number (Ca), which represents

the ratio of viscous force to capillary force, as modified by Patel, N., et al [44] :

Ca∗ =
µV
γcosθ

(2.12)

where µ, V and γ are the resin viscosity, the macroscopic resin velocity, and the

resin surface tension, respectively, and θ is the resin–fiber contact angle. As is

31



Figure 2.11: (a) A schematic of dual-scale void formation; and (b) The op-
timum capillary number for achieving minimum voidage, all
reproduced from [43]

evident in Figure 2.11, small capillary numbers lead to meso-voids, and large

ones cause micro-voids. Hence, minimization of both these voids can be accom-

plished by selecting an optimal capillary number.

As we can clearly see, the meso and micro voids show more complicated

physics in their formation and evolution, compared to the dry spots. A large

amount of research has been performed in this area, by creating multi-scale me-

chanics models as well as creating a dual-scale preform framework for the study

of these voids. Furthermore, because of the associated length scale, it would re-

quire advanced testing and analysis to study such void formation.

In comparison, Macro voids or dry spots are visible to the naked eye. Their elim-

ination becomes increasingly important in terms of mechanical performance

[18]. A good amount of research has been conducted in getting rid of dry spots

through automated injection strategies and optimized process parameters. Our

research focuses on eliminating or reducing the problem from a materials per-

spective.

32



2.2.2 Causes behind Dry Spot formation

The work done by Park, C.H., et al [42] specified the various reasons behind dry

spot formation. In macroscopic flow, the main cause behind dry spot forma-

tion is mechanical air entrapment, as shown in Figure 2.12. These air bubbles

naturally get trapped when the resin flow inside the mold encounters certain

hurdles that create non-uniform progression. It is also worth noting that in the

real-life RTM processes, the resins undergo gelling as they fill the mold. This

prevents certain dry spots from being completely impregnated with resin. Es-

pecially if resin starts flowing out of the outlet, which is the conventional indus-

trial indicator for a filled part. There are several reasons behind this mechanical

entrapment, such as premature gelling of the resin, air entrapment due to the

merging of multiple flow fronts, irregular preform permeability and the race

tracking effect. The premature gelling of resin depends extensively on the ma-

terial properties of resins and their respective gel times, and is outside the scope

of this body of work. It has been shown by Laurenzi and Marchetti [45] that pre-

mature gelling can be avoided by optimizing the injection velocity so that fill is

fast enough to avoid gelling but slow enough to ensure proper impregnation

of the fibers. The other causes of void formation are discussed in the following

sections in detail.

Merging Flow Fronts

One of the key aspects of RTM involve the creation of complex and flexible

geometries. These complex geometries often requires the presence of various

33



Figure 2.12: An experimental image of dry spots formed in a composite
part, reproduced from [46]

inserts, ribs, cores, etc. during the fill process, and can even lead to heterogene-

ity in the preform lay-up. These can cause the resin flow front to split or branch

apart. When the flow fronts merge again in the absence of an outlet vent, air

gets entrapped in that region and leads to the formation of dry spots. Liu, B.,

et al [47] showed an example of dry spot formation with a flat panel consisting

of two open windows in the path of resin flow. These windows essentially stop

air from displacing and causes the resin flow front to branch (Figure 2.13). The

branching of the resin flow front leads to the formation of a dry spot.

To overcome this defect, many approaches are possible. Their work showcased

the advantages of relocating gates and vents to better control the branching out

of the resin flow front. Similarly, a sensor approach has also been trialed, in

which a secondary gate opens when a sensor is triggered by the presence of

resin in that region (shown in Figure 2.14).

34



Figure 2.13: Dry spots due to
merging flow fronts, from [47]

Figure 2.14: Use of sensors and gates
to avoid dry spots, from [47]

Irregular Permeability

As discussed earlier, the permeability in fibrous preforms is often anisotropic,

and varies with direction of flow. The varying permeability impacts the flow

of the resin, as it will impregnate areas of high permeability faster while the

resin in areas of lower permeability would be considerably slowed in move-

ment. The irregularity of permeability is often a problem in parts consisting of

areas of high density, which have been incorporated to increase the strength of

the final part in that particular region. Liu, B., et al [47] modeled an analogous

2D square panel, which includes an area of 100 times less permeability than

the surrounding regions. The variation of permeability causes the flow front to

seemingly wrap around this area as it is more difficult to penetrate, leaving a

dry spot behind (shown in Figure 2.15). Poodts, E., et al [11] displayed the same

phenomena in the fill of a con-rod beam consisting of an anisotropic fibrous

35



Figure 2.15: Dry Spots due to
Irregular Perm., from [47]

Figure 2.16: Resin volume fraction in areas of
varying permeability, from [11]

preform where regions of high and low permeability have been manually intro-

duced (shown in Figure 2.16). In these situations, optimized location of inlets

and outlets greatly reduce the presence of voids (such as starting fill from the

center of the less permeable zone in Figure 2.15). A more prevalent display of

the same issue is a phenomena known as Race-tracking, and immense research

efforts have been dedicated to reduce their effect, as we see in the next section.

Race-tracking

We have discussed how appropriate gating and venting solutions can help

greatly reduce dry spot formation in RTM. This has been widely shown in re-

search, such as the work done by Kang, M.K., et al [48] in displaying how the

use of multiple sequentially controlled injection gates can avoid merging flow

fronts and air entrapment. However, pre-form variability and permeability is

often outside our control space, and have been shown to cause the formation

of dry spots. Even with a perfectly isotropic fibrous preform, one of the big-

36



ger issues faced in RTM that is difficult to solve is known as the Race-tracking

effect. This arises between the edges of the preform and the walls of the mold it-

self [2]. Even with perfect compaction, there is always a small gap that remains,

which is known as the race-tracking channel. Essentially, open channel flow of

resin takes place in this gap, which is a few orders of magnitude faster than the

resin flow through the porous media. This may lead to incomplete saturation,

forming dry spots in the composite. While fibrous reinforcements and gaps in

the mold might vary greatly in the industry, the problem of race-tracking is one

which is universal to RTM. Therefore the optimization of the race-tracking phe-

nomenon is of high importance for the manufacturing process of any composite

material. The formation of a dry spot due to race-tracking is shown in Figure

2.17.

Figure 2.17: A representative experimental image of a dry spot formed
due to race tracking in an E-glass/epoxy composite, repro-
duced from [49]

37



We have already established how in this gap area the resin usually meets

lower flow resistance than the fibrous reinforcement, which signifies that the

permeability along a racetracking channel is higher. It is worth understand-

ing how the thickness of such a gap would affect the permeability of the race-

tracking region. Hammami, A., et al [50] were able to solve the Navier-Stokes

equation for an open channel in terms of Darcy’s Law to provide us with an

equivalent average permeability for the race-tracking channel:

Keq =
h2

12
(2.13)

where h is the clearance between the edge of the fibrous reinforcement and the

walls of the mold. It is clearly seen that a clearance as low as 1 mm would

give us an equivalent average permeability of 8.3 × 10−8m2, which is consider-

ably higher than the average permeabilities of commerically available preforms

(10−9 − 10−11m2). Hammami, A., et al [50] showcased the effect of a 3mm clear-

ance in Figure 2.18.

Naturally, the question arises as to what this gap thickness h depends on and

how it can be reduced. Koutsonas, S. [2] was able to show how h varied with

the height of the mold HMold and compaction pressure P for a 90◦ curved plate:

h = HMold − α(P−1)β (2.14)

where α and β are material constants.

38



Figure 2.18: A showcase of racetracking effect through (a) a simulation
and (b) an experiment of fill with a clearance of 3 mm between
the mold wall and preform, all reproduced from [50]

2.2.3 Replication in FLUENT Simulations

In order to create test geometries for the analysis of dry spots, we first need to

validate whether our simulations are able to capture dry spot formation accu-

rately. For this, we rely on the images and schematics shown in the literature

reviewed in the previous section. Most of the voids shown in literature are ei-

ther micro or meso voids. Dry spots are shown rather sparingly, as the industry

recommendation is centered around proper injection strategies involving mul-

tiple gates and sensors. Even the schematics and images shown above provide

very little information that is essential for our simulation set up (such as pre-

form permeability, porosity, resin viscosity, injection pressure, geometry speci-

fications etc.).

39



We have been able to simulate the following dry spots, shown in Figure 2.19,

based on a trial-and-error method to determine the process parameters that

have been used in the schematics and experiments of Liu, B., et al [47] and

Hamidi, A., et al [49], which gives us approximately similar results in terms of

contour and shape. To replicate the dry spots, shown in Figure 2.19 (a) and (c), a

permeability of the order of 1× 10−9 is chosen, while the central region in Figure

2.19 (b) has a lower permeability by two orders of magnitude. The resin chosen

is Newtonian of µ = 0.07115Pas. For the creation of Figure 2.19 (c), we define

a race-tracking channel of arbitrary width (shown by the solid black lines) and

utilize a Newtonian resin of viscosity µ = 0.37Pas. Both the Race-tracking and

Merging Flow Front void formation cases have utilized a low pressure differ-

ence of 10,000 Pa across the inlet and outlet domain length of 1cm (along Y-axis)

to showcase the formation of the respective voids. This is particularly of signif-

icance in Chapter 3, as we shall see later.

40



Figure 2.19: Dry spots formed in simulations by : (a) Merging flow fronts
(based on Fig.2.13) (b) Irregular Preform Permeability (based
on Fig.2.15) (c) Race-tracking effect (based on Fig.2.17)

41



CHAPTER 3

NON-NEWTONIAN FLOW THROUGH POROUS MEDIA

In all the RTM simulations we have come across, the resins used are Newto-

nian in nature, i.e. their viscosities remain constant at all shear rates. This is

also specified in the assumptions related to the governing equations (Equation

2.3). Our novel approach to the problem of void formation is one from a ma-

terials perspective. We take advantage of the shear thinning or shear thicken-

ing behavior of certain fluids, known as non-Newtonian behavior. Essentially,

the viscosities of these proposed resins would change with the strain rates they

face. In the cases of dry spots, we have seen the formation of areas of varying

resin flow-front velocity, either due to irregular permeability or race-tracking. It

would be advantageous to have a less viscous resin in low resin velocity regions,

or vice-versa, to ensure that resin flow is uniform throughout the mold. Non-

Newtonian resins can be utilized to achieve the above-mentioned fluid charac-

teristics, however, a comparison of the degree of success needs to be established.

In itself, non-Newtonian flow through porous media is a research area that is ex-

tremely complex, as we shall see in this chapter.

3.1 Non-Newtonian Fluids

A Newtonian fluid is one which obeys Newton’s law of viscosity, which states

that the fluid shear stress τ is directly proportional to the rate of shear γ̇ [51]:

τ = µ
∂γ

∂t
= µγ̇ (3.1)

The Newtonian viscosity µ is independent of the shear rate γ̇. However, that is

42



not the case for non-Newtonian fluids, for which the shear stress τ does not fol-

low the linear law. Non-newtonian viscosities are known as apparent viscosity

η and are functions of the shear rate γ̇. This leads to a modified Equation 3.1:

τ = η(γ̇)γ̇ (3.2)

The most commonly used model for non-Newtonian viscosity behavior is given

by the power law:

η = Kγ̇n−1 (3.3)

Where K is known as the consistency index and n is the power law index. A

value of n < 1 shows shear thinning behavior and n > 1 shows shear thicken-

ing behavior. Time independent non-Newtonian fluids can be categorized as

purely viscous (shear thinning, shear thickening and Newtonian) and viscoplas-

tic (Bingham fluid). Time dependent fluids can be classified as rheopectic and

thixotropic fluids. The various classifications are shown in Figure 3.1.

Figure 3.1: Various types of non-Newtonian fluids, both time-independent
(left) and time-dependent (right), all reproduced from [51]

43



3.1.1 Carreau Model

Multiple fluid models exist that accurately capture non-Newtonian behavior

over a large range of shear rates. These are mainly phenomenological mod-

els that are fitted to rheometer data. They include Non-Newtonian Power Law

model, Cross Model, Herschel Buckley model, Carreau model, to name a few.

Out of all these, the Carreau model is found to be a popular model for cap-

turing realistic viscosity at very low and high shear rates [52]. According to

Hackley and Ferraris [53], the Carreau model is able to accurately showcase

shear-thinning flow with asymptotic viscosities at zero and infinite shear rates,

and with no yield stress, as shown in Figure 3.2.

Figure 3.2: The rheological behavior of a shear-thinning fluid specified by
the Carreau model, reproduced from FLUENT Documentation

The Carreau model is able to successfully negate the concerns associated

with the power law model (governed by Equation 3.3) with regards to viscosity

at very low or high shear rates, by incorporating a Newtonian plateau in each

of those regions. Hence, four parameters define the Carreau model for a shear

thinning fluid- the zero shear Newtonian viscosity η0, the infinite shear Newto-

nian viscosity η∞, the critical shear rate γ̇c at which shear thinning begins and the

slope of the power law region that is proportional to (n − 1). The corresponding

44

https://www.afs.enea.it/project/neptunius/docs/fluent/html/ug/node297.htm#sec-viscosity-newt-power


equation for apparent viscosity η is given by [54]:

η = η∞ + (η0 − η∞)(1 + (λγ̇)2)
n−1

2 (3.4)

where λ is known as the characteristic time constant and is simply the inverse of

critical shear rate γ̇c, at which non-Newtonian behavior begins. We will use the

Carreau model to simulate non-Newtonian flow as it accurately represents the

real-world behavior of many shear thinning fluids [55]. It has also been shown

to explain various polymeric resin properties, as established by Kol, R. et al [56].

Figure 3.3: a) The viscosity curves at different temperatures of an ideally
mixed resin, reproduced from [57]; and b) Our least squares
curve fitting of a Carreau Model to the available experimental
data from [57]

We are able to show good agreement of a fitted Carreau model with the ex-

perimental data from the work of Meier, R., et al [57]. They have studied an

epoxy resin system (EPIKOTETM Resin MGS® RIMR135 and EPIKURETM Cur-

ing Agent MGS® RIMH1366) that is frequently used in RTM for the manufac-

ture of wind turbine blades, and have obtained data from experimental tests in

the rheometer at various temperatures. We are able to fit a Carreau model with

some of the more detailed experimental data through a least-squares curve fit-

ting method, shown in Figure 3.3. It is worth mentioning that Meier, R., [58]

45



utilized a modified version of the Carreau model known as the Bird-Carreau-

Yasuda model [59], which consists of an additional transition parameter a, and

established the model parameters as η0 = 53.37Pas, η∞ = 0.37Pas, λ = 55.76s,

n = 0.0268 and a = 8.11. We shall these parameters as points of reference

throughout our work.

3.2 Flow through Porous Media

Our research objective is to optimize fill by varying the parameters of the above

mentioned Carreau model. However, a big question arises with regards to the

scale of the simulation. We are not explicitly modeling the porous media geom-

etry, but we are characterizing the porous media space by a given permeability

and porosity. The question that arises is with regards to the modeling. If we

want to approximate the fibers as a porous media space captured by parame-

ters (permeability and porosity), naturally the question arises as to whether we

can use the same homogenized equation used for Newtonian fluids.

As we recall, Darcy’s Law (Equation 2.4) consisted of a viscosity term µ that

was considered Newtonian and constant for all processes, irrespective of veloc-

ity v⃗, pressure gradient ∇P and permeability of the preform K⃗. This is expected

to change as the apparent viscosity for non-Newtonian fluids depends on the

shear rate (Equation 3.2), and the shear rate itself is a non-linear function of the

pressure gradient, as shown by Fayed, H.E., et al [52]. Even though our op-

timization deals with the same permeability throughout, it is natural that the

shear rates faced by the fluid would also change with changes in permeability,

as it is a measure of resistance to flow. This has been shown by Meier, R., [58],

46



who verifies that changing fiber volume fraction (i.e. changing permeability)

affects the maximum shear rates encountered. Overall, this becomes a problem

worth analyzing through prevalent research and our own simulation models

before the optimization of parameters.

3.2.1 Literature Review

Rheometer-measured steady state viscosity often does not entirely encapsulate

the behavior of polymeric fluids or various flow conditions in such complex

tortuous geometry. Skauge, A., et al [60] were able to show various phenom-

ena such as mechanical degradation, extensional viscosity and the effect of

contraction-expansion channels that result in significantly different rheological

behavior compared to bulk flow.

The work done by Sochi, T. [61] showcased the complex nature of flow in porous

media and how it involves many approximations. After careful analysis, they

were able to conclude that no analytical fluid mechanics solution exists for the

problem at hand, and is unlikely to appear in the near future. However, they

were able to highlight the various methods and models used in practice to pre-

dict non-Newtonian flow behavior in porous media, along with their corre-

sponding qualities and drawbacks. Some of these include Continuum models,

Capillary Bundle models, Numerical methods and Pore-network modeling.

Eberhard, U. et al [55] confirmed that the relation between Darcy velocity and

pressure drop cannot be described by the linear function mentioned in Equation

2.4 for Newtonian fluids. Instead, they suggested the use of a bulk equation

47



approach so that flow velocity is linear in the pressure drop. This is possible

through the incorporation of an effective viscosity µe f f across the length L of the

domain of concern, as shown below:

v =
K
µe f f

∆P
L

(3.5)

The models suggested by Cannella, W.J. [62] and Hirasaki and Pope [63] to cal-

culate effective viscosity are empirical and semi-empirical in nature. They use

a capillary bundle model to estimate an effective shear rate γ̇e f f by comparing

the non-Newtonian power law (Equation 3.3) with the equations of Newtonian

Poiseuille flow, which gives us the following equation:

γ̇e f f =
1
√
ζ

(
3n + 1

4n

) n
n − 1 4v√

8Kϕ
(3.6)

where n comes from the non-Newtonian model, ζ is the tortuosity and γ̇e f f

varies linearly with Darcy velocity v. The effective viscosity is found by re-

placing the effective shear rate (Equation 3.6) in the respective non-Newtonian

constitutive model. Hirasaki and Pope [63] incorporated the semi-empirical

tortuosity of flow ζ, while Cannella, W.J. [62] argued the use of an empirical

correction factor (C-factor) instead of (1/
√
ζ), which is found to vary over sev-

eral orders of magnitude. Eberhard, U. et al [55] pointed out the flaws in this

approach when using a more complicated non-Newtonian model such as the

Carreau model, but both their non-empirical approach and experimental data

showed similar results to Hirasaki and Pope’s method [63].

It is clearly understandable how varied and complicated non-Newtonian flow

48



through porous media is, with various empirical models, numerical methods

and the distinct lack of a reliable analytical solution. The use of Equation 3.6

is able to approximately answer a lot of the complications, and Sorbie, K.S. [64]

was also able to establish through 2D connected pore network modeling that the

average shear rate correlates linearly with the flow rate. This is not withstand-

ing some of the related issues with a capillary bundle approach, which include

neglecting complex features such as tortuosity, pore size distribution, intercon-

nectivity between pores, anisotropy, variation in cross-sectional area and exten-

sional viscosity [60].

3.2.2 Porous Media Simulations

Simulation Test Case

As an initial test case to study the non-Newtonian flow through a porous media,

we create two representative 2D unit cells of the same fiber area fraction v f =

0.342 and dimensions 1 × 1 cm, based on a capillary bundle approach shown

in Figure 3.4. The diameters of the grains or circular representations of fibrous

media are 0.02 cm for the first test case and three times that (0.06 cm) for the

second test case, which we will refer to as Model A: 1x Fiber and Model B: 3x

Fiber test case moving forward. Fiber area fraction is calculated by dividing the

combined area of all the fibers by the total area of the 2D unit cell, which gives

us v f = 0.342 for both Model A and Model B. A Newtonian resin of viscosity µ =

0.37Pas is injected through the left edge (inlet) and exits through the right edge

(outlet) of each unit cell in multi-phase flow, across which pressure differences

are varied from 0 to 700 kPa.

49



Figure 3.4: Our porous media representative unit cells for simulation, la-
belled as: a) Model A: 1x Fiber and b) Model B: 3x Fiber

We are able to estimate the permeability of Model A and B (described above)

as K = 9.86 × 10−10m2 and K = 5.70 × 10−9m2 respectively, by computing the av-

erage outlet velocity at various pressure gradients and applying Darcy’s Law

(Equation 2.4) to the respective plots, shown in Figure 3.5 (a). It is worth noting

that Model B is more permeable, despite having the same fiber area fraction as

Model A, due to larger free space. The permeability calculations are based on

1D rectilinear flow. We are also able to validate these calculations by varying the

Newtonian viscosity µ, which gives us the same permeability each time. We use

a fictitious Carreau Model that shows shear thinning below our current Newto-

nian threshold (µ = 0.37Pas) to model a non-Newtonian fluid. The parameters

are as follows- η0 = 0.37Pas, η∞ = 0.0037Pas, λ = 1 × 105s and n = 0.8. The

values of λ and n are chosen to ensure we have a large shear thinning regime, as

shown in Figure 3.5 (b), so that the resin always undergoes shear thinning at all

our chosen pressure gradients to fully demonstrate non-Newtonian effects. The

range of shear rates are shown in Figure 3.5 (b) for both test cases, confirming

50



Figure 3.5: a) The comparison of Newtonian flow between Model A: 1x
Fiber and Model B: 3x Fiber; and b) The fictitious Carreau
model created to analyze non-Newtonian flow

that we are in the shear-thinning regime throughout our simulations. We are

Figure 3.6: A plot comparing the average outlet velocity vs. the pressure
gradient for non-Newtonian flow in Model A (shown in black)
and Model B (shown in red)

able to compare the outlet velocity vs. pressure gradient profile for both cases

in Figure 3.6. The plot confirms the fact that velocity is no longer linear with

pressure gradient, as is characteristic for non-Newtonian flow through porous

media. We are also able to show that for the same fiber area fraction, differ-

51



ent permeabilities lead to different shear rates and consequently, different flow

variables.

Isotropic Fibers

Now that we have established non-Newtonian behavior in our test simulations,

we model a system more closely related to the fiber bundle models encountered

in RTM to use for our optimization. According to Abraham, D., et al [17], the

achievable fiber volume fraction for RTM lies between 0.3 to 0.65. Contrary to

our simulation test case, the fibers in an actual preform would be less ordered.

The aim was to develop a more realistic isotropic fiber bundle model, where the

density of fibers is the same throughout the representative unit cell. Through

the use of a Poisson point process, we were able to develop stochastic distri-

butions of 2D fiber cross-sections. Fiber volume fraction is modulated by the

number density of the fibers, which is constrained by certain factors such as the

maximum number of fibers that can be packed into our unit cell and still be

meshed accurately for simulations. For this reason, we limited ourselves to a

maximum fiber volume fraction of 0.5, up to which reasonable mesh element

size could be incorporated. We also slightly reduced our domain of measure-

ment along the horizontal axis to avoid convergence errors caused by skewed

elements near the inlet/outlet of flow. A representative unit cell of v f = 0.392 is

shown in Figure 3.7 (a). We generated 7 unit cells in the above-mentioned fiber

volume fraction range- v f = 0.295, 0.327, 0.362, 0.392, 0.421, 0.455, 0.483- to use

for our study and optimization. We were also able to show similar flow results

upon increasing the number density of fibers while keeping the fiber volume

fraction the same, as shown in 3.7 (b), which confirms the accuracy of our fiber

52



Figure 3.7: a) A representative fiber bundle unit cell of generated v f =
0.392; and b) A plot of the average outlet velocity vs. the pres-
sure gradient for v f = 0.425 fiber bundle, compared at different
fiber number densities

generation method. We then proceeded to normalize our domain with regards

to real-life RTM processes. For this, we chose a carbon fiber preform (commonly

used in aerospace applications), where the fibers have an average diameter of

10 µm. The pressure gradients were normalized accordingly, with the max pres-

sure gradient being limited at 2 × 107Pa/m. This was specified by Bodaghi, M.,

et al [65], who mentioned that higher injection pressures may cause preform de-

formation and difficulties in mold filling.

Using the same process as in the test cases, we were able to evaluate the per-

meabilities by utilizing a Newtonian resin of viscosity µ = 0.37 and plotting

the average outlet velocity against the pressure gradient, as shown in Figure 3.8

(a). We then use the same Carreau model specified in Figure 3.5 to compute

the non-Newtonian velocities at different pressure gradients, as shown in Fig-

ure 3.8 (b). We are able to see the same non-linear relationship between velocity

53



and pressure gradient for all our fiber bundles.

Figure 3.8: The flow variable plots of all fiber bundles, including: a) New-
tonian velocity vs. pressure gradient; and b) non-Newtonian
velocity vs. pressure gradient

Carreau Model Correction

For choice of porous media for optimization, we go with our established fiber

bundle of v f = 0.392 (shown in Figure 3.7 (a)) as it lies around the average value

of our chosen range for porous media simulations (v f = 0.3−0.5). It has a poros-

ity of 0.608 (Equation 2.6) and a corresponding permeability of 1.98 × 10−9m2

(computed from Figure 3.8 and confirmed with the Kozeny-Carman equation

[29]). Tortuosity of this porous media is found to be ζ = 1.13 after testing with

fluids of various viscosities and using the following equation [66]:

ζ =
< u >
< ux >

(3.7)

where < u > is the average magnitude of the velocity over the entire volume

54



and < ux > is the volumetric average of its component along the macroscopic

flow direction (X-axis).

As we had seen in Appendix C, the pressure gradient in our simulation will

Figure 3.9: Effective viscosity vs. Effective Shear rate for a v f = 0.392 fiber
bundle, compared to the Carreau model used in the simulation

reduce across the medium from the initial set value between the inlet and the

outlet. In other words, a porous media cell at the inlet will face the maximum

pressure gradient prescribed in our simulation model, and the rest of the cells

will all face lower pressure gradients. Our chosen pressure gradient for our

void formation simulations (Figure 2.19) is 10 kPa across a domain of length 1

cm, giving us a safe pressure gradient threshold of 1 × 106Pa/m [65]. The effect

of this is shown in Figure 3.9, where we trial various pressure differences rang-

ing up to 10 kPa for our given Carreau model (Figure 3.5 (b)). The simulations

through our chosen porous media give us the average outlet velocity at various

pressure gradients, and we can calculate the effective viscosity at each pressure

gradient using Equation 3.5. We then plot the effective viscosity against the

effective shear rate (from Equation 3.6), which depends on our choice of fiber

bundle. We then compare our plot with the Carreau model we had chosen for

55



our simulation (Figure 3.5). As we can see in Figure 3.9, the resulting plot is lat-

erally shifted to the left of our Carreau model, as confirmed by both Eberhard,

U., et al [55] and Skauge, A., et al [60]. Our plan for optimization will be to run

through our entire sample space of Carreau model parameters, and incorporate

a shift to the right of the best-performing Carreau model to account for the ef-

fect of non-Newtonian flow through porous media. There are associated limits

to the approach of using a representative capillary bundle model to account for

flow through porous media, as we have discussed in Section 3.2.1.

56



CHAPTER 4

STATISTICAL OPTIMIZATION OF PARAMETERS

4.1 Baseline Model setup

Now that we have established our simulation model and carefully explored

non-Newtonian flow through porous media, our final objective is to optimize

the parameters of the Carreau model to show the best overall fill. For this, we

create a base-line model for one of the void formation cases mentioned in Chap-

ter 2. We use the other models to test the effect of our optimized results. We

define our parameter sample space based on currently available data on var-

ious resins, polymers and related fluids of interest. We then select one of the

fiber bundle models from our study of non-Newtonian flow through porous

media, and use the corresponding material properties that we are able to de-

rive. Finally, we run an informed-search based optimization over our entire

sample space based on certain target values. Once we have settled on our best

set of parameters, we can test their effect on the other models of void formation,

and apply the chosen correction factor for porous media flow. In this way, we

are able to successfully define a non-Newtonian fluid model that gives us the

best fill in a known case of void formation.

4.1.1 Parameter Sample Space

Before setting up our model, we need to define a realistic parameter space to op-

erate in. We are able to define this by exploring the viscosities of resins used in

both RTM and non-RTM related processes as well as polymers of various types,

57



especially those that are non-Newtonian in nature. Different solvents have also

been studied to establish a variety of Carreau model parameters. Our aim is to

explore resins that are not commercially available as well, in case they are able

to show better fill properties.

Some of the explored materials include- High Density Polyethylene (HDPE)

Resins; Semipolymerised 52% Linseed Oil-glycerol-phthalate Polyester Resin

In Different Solvents; Bio-based Phenolic Resin; Polylactide (PLA) in

Methyl Methacrylate (MMA); Polylactic Acid; Partially Hydrolyzed Polyacry-

lamide (PAMH); Hydrolyzed Polyacrylamide Solution (HPAM); Poly (Lactic

Acid)/Poly (Ethylene Oxide)/Carbon Nanotubes Biosensor; Carboxymethyl-

cellulose (CMC); Xanthan Gum (XG); Isotactic PP Borflow HL504FB (76k),

HL508FB (64k) and HL512FB (56k); Nano-particle Filled Polymer Nanoparti-

cle Composites; 5 Wt% Polyisobytylene (PIB) In 95 Wt% Tetradecane (C14). On

the basis of the resulting viscosity parameters, we are able to set up our sample

space, shown in Table 4.1. It is worth noting that for a shear thinning case, η0

will always be greater than η∞. Values of n = 1 signify Newtonian behavior and

values of n > 1 signify shear thickening, as we shall discuss in the next section.

Carreau Parameter Range of established values

λ 10−4 − 105

n 10−4 − 2
η0 10−2 − 102

η∞ 10−3 − 102

Table 4.1: The Carreau model parameter sample space, obtained from an
exploration of various materials

58



4.1.2 Shear Thickening

Purely shear thickening fluids are not very common in practice. Dilatant behav-

ior is mostly seen in colloids in suspensions, and very rarely in polymers. The

limited data reported so far suggests that flow behaviour may be represented

by the power-law model (Equation 3.3) where n is taken to be greater than one.

It is not yet possible to ascertain whether these materials also display limiting

apparent viscosities, such as the Newtonian plateau shown at high shear rates

in Figure 3.2. However, it is a sample space worth exploring to see if void for-

mation is reduced. With that in mind, we take a look at some simulation ar-

eas that have incorporated non-Newtonian models that undergo shear thicken-

ing [67–71]. The value of n is found to vary from 1 to 2 with no infinite plateau,

justifying our sample space in Table 4.1. We can see a sample Carreau model

for shear thickening with all parameters staying the same as those in Figure 3.5,

except n which is chosen as 1.5, shown in Figure 4.1.

Figure 4.1: A fictitious shear thickening Carreau model with parameters
η0 = 0.37Pas, η∞ = 0.0037Pas, λ = 1 × 105s, n = 1.5

59



4.1.3 Model set up

We choose the simulation corresponding to the Racetracking effect (Figure 2.19

(c)) as our baseline model, shown in Figure 4.2 (a). The target function of our

optimization is the Volume-averaged resin volume fraction, which we will re-

fer to as Resin vf, and is given by the equations in the FLUENT Documentation.

This acts as a good parameter to analyze resin volume fraction from a software

perspective, as it computes the volume-weighted resin volume fraction of each

cell in the domain and divides it be the total volume of each cell. We can com-

pute the estimated void content by simply calculating (1 - Resin vf). Our max-

imum time duration for flow is taken as 10 minutes, while our simulation is

terminated as soon as the outlet reaches 95% fill. This acts as a metric to con-

firm that resin will be exiting the mold, as compared to the industrial process.

We will be performing a separate optimization for Shear Thickening utilizing

a user-defined function in ANSYS FLUENT. Finally, we shall compare our best

results with our other baseline model for void formation, based on an altered

version of Figure 2.19 (a), which is shown in Figure 4.2 (b).

4.2 Bayesian Optimization

Finally, we take a quick look at the basis of our optimization technique. The

problem at hand is one which is known as a black box function in Machine

Learning. This is because there is no discernible pattern or trend with regards

to our target function (Resin vf). It is very possible that a highly viscous fluid

would show better fill than a less viscous fluid. This also explain why we are

unable to use a gradient based approach, as the behavior of the target function is

60

https://www.afs.enea.it/project/neptunius/docs/fluent/html/th/node418.htm#sec-compute-volint


Figure 4.2: a) Our racetracking base-line model at Newtonian µ = 0.01,
Resin vf = 0.9785; and b) our merging flow front base-line
model at Newtonian µ = 0.37, Resin vf = 0.8127

unknown to us. Furthermore, from a statistical perspective, our problem is sus-

pected to be a multi-modal function. In essence, a multi-modal function would

signify the presence of multiple local maxima and minima across the domain,

while our research objective is to find the global maxima for resin volume frac-

tion.

There are multiple methods that deal well with multi-modal functions, as de-

fined by Roman, I., et al [72]. These include Evolutionary algorithms, Simulated

annealing, Differential evolution and Bayesian Optimization. Out of all of these

methods, Bayesian Optimization is preferred as it is an informed-search algo-

rithm, which utilizes prior information to make an educated guess as to where

the maxima would be. The method is based on a probabilistic approach that

uses a surrogate model to approximate the target function as well as an acquisi-

tion function to balance exploration and exploitation of the sample space.

61



Of course, there are various drawbacks to the use of Bayesian Optimization.

This includes computational costs for high-dimensional problems, sensitivity

to choice of hyperparameters, limited scalability and a lack of theoretical guar-

antees for any results generated. However, a probability based approach is our

best bet for exploring a sample space of several orders of magnitude for a target

function that is unknown in behavior.

Search Algorithm

We define how the surrogate model and acquisition function work together to

explore the parameter space. In Figure 4.3 (a) we see a hypothetical target func-

tion. The optimizer does not know what the target function looks like. The sur-

rogate model is formed based on sampled points and is represented by Gaussian

Processes that return functions fitted to the sampled data points. The process re-

turns several functions, each of which have probabilities attached to them. The

Gaussian Process is a probability distribution that can be updated with new in-

formation through inherently probabilistic Bayesian theories, shown in Figure

4.3 (b). After each iteration, the optimizer continues to look at the current sur-

rogate model, learning more about areas of interest by sampling, and updating

the model accordingly. After a certain number of iterations, the optimizer ar-

rives at a surrogate model that accurately captures the target function, and is

able to define a global maxima. The number of iterations depend on the size of

the sample space and the behavior of the target function. The model updates on

the basis of Bayes’ Theorem, which is given by [73]:

p(S core|Hyperparameters) =
p(Hyperparameters|S core) ∗ p(S core)

p(Hyperparameters)
(4.1)

62



Figure 4.3: A demonstration of the surrogate model estimation based on a
Bayesian probability approach, reproduced from The Beauty of
Bayesian Optimization [73] and Hyperparameter Optimization [74]

The surrogate model — represented as a probability distribution known as the

prior — is updated with an acquisition function. This function is responsible

for driving the sampling of new points to test, in an exploration and exploita-

tion trade-off. Exploration seeks to sample in locations where the uncertainty is

high. This ensures that no major region of the space is left unexplored. Exploita-

tion seeks to sample where the surrogate model predicts a good target, looking

for the maxima. Common acquisition functions include expected improvement,

maximum probability of improvement, upper confidence bound, etc. These are

explained in detail in Washington University in St. Louis’ course on Bayesian

Optimization Acquisition Functions [75]. We choose to go with the upper con-

fidence bound function, which is typically described in terms of maximizing a

target function. It follows the given equation:

aUCB(x, κ) = µ(x) + κσ(x) (4.2)

where µ(x) is the mean of the target function f (x) that promotes exploitation,

while σ(x) is the standard deviation that promotes exploration. Here, κ is a

63



trade-off parameter controlled by the user, the value of which decides which

one of the two, exploitation or exploration, takes precedence. Generally, a value

of κ = 2 signifies a balance of both.

Hence, our approach involves using the built-in Bayesian optimizer in Python

libraries to explore and exploit the same space. An example of Bayesian Opti-

mization in action can be found on GitHub [76]. Since we are dealing with a

large sample space that cannot be approximated by sampling a few points, we

let the optimizer decide where to go, based on Equation 4.1, instead of sampling

an initial set of points. We also vary the value of κ, as seen in the results section,

to first explore the entire space and then exploit regions of potential global max-

ima, in order to avoid the optimizer getting trapped in a local maxima.

4.2.1 Results

Shear Thickening

With regards to our Shear Thickening optimization, we restrict the maximum

possible viscosity of our Carreau model to η = 1Pas, as shown in Figure 4.4.

There are two reasons behind this. We create an artificial Newtonian plateau at

infinite shear rate to prevent the viscosity ratio from becoming too high, which

leads to the failure of the VOF model (Shown in Appendix A). We also found our

initial shear thickening simulations at higher viscosity than η = 1Pas to show an

increased fill and simulation run time. We vary the trade-off parameter κ over

the entire sample space, as shown in Table 4.2. We start with a balanced value

of κ = 2, which proves to be insufficient for exploring the large parameter space.

64



Figure 4.4: Our altered Carreau model for shear thickening, based on the
parameters of Fig. 4.1, with a generated Newtonian plateau at
η = 1Pas shown by the green line

We then ramp up exploration of the sample space by using values of κ = 10,

κ = 10, 000 and κ = 100, 000. We also break down the range for λ into two halves

of separate optimization, as it varies over quite a few orders of magnitude and

we want to ensure that the entire space is well explored. Once we have explored

our space sufficiently and loaded the results into our optimizer, we begin a more

narrow exploitation through the use of κ = 0.001 and an exploration of concen-

trated ranges around which best results are seen.

After 1000 iterations, we are able to show confidence in our results by showcas-

ing the variation in the target function (Resin vf) as we move through the iter-

ations, which shows initial perturbations followed by circling around the max-

imum values. The last 150 iterations are only exploring the Newtonian space,

which explains the relatively similar results. We have also shown the corre-

sponding fill times, which showcase how much they vary as we move through

the iterations. It is worth noting that the fill times do not reach very high val-

ues as our Newtonian plateau at infinite shear rate is very low (η = 1Pas). Our

65



Figure 4.5: a) All our shear thickening optimization results and corre-
sponding fill times vs. number of iterations; and b) our best
results (where Resin vf > 0.9999) vs. corresponding fill times

Figure 4.6: Based on our shear thickening optimization, we have the corre-
sponding resin volume fraction contours of: a) our best overall
result and b) our best non-Newtonian result; the parameters of
which are given in Table 4.3

66



κ Value Parameter Range Iterations

2 Full Range 100
10 Full Range 25
10,000 Full Range 25
100,000 Full Range 25
10,000 10 < λ < 1 × 105 50
10 1 × 10−4 < λ < 10 50
10,000 1 × 10−4 < λ < 10 50
0.001 Full Range 100
2 Concentrated around max 325
0.001 Concentrated around max 100
2 Only Newtonian viscosities 150

Table 4.2: The variation of our optimization parameters for the shear thick-
ening optimization of 1000 iterations

results show that for the given geometry (Figure 4.2 (a)) Newtonian fluids fill

better in comparison to shear thickening non-Newtonian fluids. The best result

and corresponding contour is shown in Figure 4.6 (a). While most of the best

results are Newtonian in nature, the best Shear Thickening result is shown in

Figure 4.6 (b). However, the value of n, shown in Table 4.3, for the shear thick-

ening fluid is marginally greater than 1, signifying that it practically behaves as

a Newtonian fluid. This can also be seen in Figure 4.7, where the value of the

slope n causes the Carreau model to be extremely gentle.

Shear Thinning

We vary the value of κ for our Shear Thinning Optimization, between κ = 0.001

to κ = 100, 000, to carefully explore and exploit the full sample space (from Table

4.1), the extent of which is shown in Table 4.4. After an initial set of explorations,

it was observed that higher values of λ (or, lower values of critical shear rate γ̇c)

lead to better fill. This clearly signifies a quick onset of shear thinning, which

67



Shear Thickening
Optimization results Rheology Parameters Fill output

Overall best result µ = 0.54 Pa s
Resin vf = 0.99999508

Fill Time = 1.3 s

Overall best non-
Newtonian result

η0 = 0.430615 Pa s

η∞ = 0.077964 Pa s

λ = 9.39285 s

n = 1.007428

Resin vf = 0.9999911

Fill Time = 1.2 s

Table 4.3: Our shear thickening optimization results, the corresponding
Carreau model parameters and fill output variables

Figure 4.7: A comparison of the corresponding Carreau models of the re-
sults shown in Table 4.3 and Fig. 4.6

68



κ Value Parameter Range Iterations

2 Full Range 25
100,000 1 < λ < 1 × 105 25
10,000 10 < λ < 1 × 105 100
10,000 1 × 10−4 < λ < 10 103
0.001 10 < λ < 1 × 105 25
10,000 10 < λ < 1 × 105, η0 < 10 100
10,000 1 × 10−4 < λ < 10, η0 < 10 100
0.001 η0 < 10 200
0.001 Concentrated around max 200

Table 4.4: The variation of our optimization parameters for the shear thin-
ning optimization of 878 iterations

Figure 4.8: a) All our shear thinning optimization results and correspond-
ing fill times vs. number of iterations; and b) our best results
(where Resin vf > 0.999) vs. corresponding fill times

is beneficial in terms of fill time as well. Furthermore, we were able to narrow

down our region of best performance to η0 values of less than 10 Pas. One

particular limitation of our exploration in spaces of larger viscosity than this

was the increased amount of time it took to run one set of simulations. This

was due to software complications of ANSYS FLUENT, which causes it to slow

down after large simulation time durations due to storage of larger amounts of

69



prior data. However, our overall best result from both optimizations was found

to be at η0 values of less than 10 Pas. We see how the optimizer moved through

the iteration space in Figure 4.8 (a). We also see the variation in terms of best

results (Resin vf > 0.999) in Figure 4.8 (b), which confirms our hypothesis that

comparatively good fill can be achieved by both low viscosity and high viscosity

shear thinning fluids.

Figure 4.9: The resin volume fraction contours of our: a) overall best result
from shear thinning optimization; and b) the corresponding in-
finite Newtonian plateau viscosity of our best result

The contour of the best result in shown in Figure 4.9 (a). We also record the

parameters of the Top 3 results from the Shear Thinning optimization, shown in

Table 4.5, for testing their effect on our other geometry (Figure 4.2 (b)). We see

that the Newtonian result from our Shear Thickening optimization (Table 4.3)

performs better than the 2nd and 3rd best result from the Shear Thinning opti-

mization, but not as well as our overall best result. The Carreau models for our

shear thinning results are shown in Figure 4.10, along with the corresponding

shear rate ranges over the entire fill. It is interesting to note that the shear rate

70



Shear Thinning
Optimization results Rheology Parameters Fill output

Overall best result

η0 = 9.086963 Pa s

η∞ = 0.409465 Pa s

λ = 98586.53 s

n = 0.574426

Resin vf = 0.9999956

Fill Time = 1.2 s

Overall 2nd
best result

η0 = 1.912514

η∞ = 0.1590934

λ = 42961.45

n = 0.649057

Resin vf = 0.9999938

Fill Time = 1 s

Overall 3rd
best result

η0 = 8.171539 Pa s

η∞ = 0.167079 Pa s

λ = 99324.41 s

n = 0.403321

Resin vf = 0.9999937

Fill Time = 1 s

Table 4.5: Our top 3 shear thinning optimization results, the correspond-
ing Carreau model parameters and fill output variables.

71



Figure 4.10: Shear rate range in the Carreau model of our: a) overall best
result; b) 2nd best result; and c) 3rd best result

range in all three results lies very close to the Infinite Newtonian plateau. This

raises the question of whether an entirely Newtonian fluid of the same viscosity

as that of the infinite Newtonian plateau would give us the same results. We

tested this by running a simulation at Newtonian viscosity µ = 0.409465, cor-

responding to the overall best result from Table 4.5, and found the Resin vf to

be only 0.99973 with a Fill Time of 1.1 s (shown in Figure 4.9 (b)). This shows

us that a certain amount of shear thinning can definitely improve fill as com-

pared to a complete Newtonian fill. Finally, we take our best result and correct

it according to our Porous Media lateral shift, given by Equation 3.6. Since we

have shown that a lateral shift to the left is expected in porous media, if our best

72



Figure 4.11: a) Our corrected Carreau model of our overall best results
with the incorporated lateral shift; and b) the shifted viscosi-
ties checked with the parameters of the corrected Carreau
model

results are shifted to the right, they will represent the rheometer-based Carreau

fluid model that would give us the best fill. An added complication is the de-

pendence of effective shear rate on the Darcy velocity v in Equation 3.6. It is is

difficult to accurately compute the shift to the right as we are unaware of the

magnitude of velocities generated by our corrected Carreau model. Our fix for

this involved computing the average lateral shift by incorporating our best re-

sult (Table 4.5) to our v f = 0.392 fiber bundle model (Figure 3.7), and correcting

our Carreau model by shifting it by that amount to the right. We find the aver-

age lateral shift to be around 25.6% of the shear rates. Incorporating the same

lateral shift to the right, we get our corrected Carreau model, the parameters for

which are given in Table 4.6 and is shown in Figure 4.11 (a). The parameters

are obtained by shifting all shear rates to the right by 25.6% and curve fitting

the Carreau viscosities at these adjusted shear rates. The curve fit shows good

agreement with the shifted results, shown in Figure 4.11 (b). We can test the

validity of our corrected Carreau model by passing a Carreau model fluid with

73



corrected parameters from Table 4.6 through our v f = 0.392 fiber bundle simula-

tion. We then plot the effective viscosity (computed using Equation 3.5) against

the effective shear rate (from Equation 3.6) from this simulation and match the

results with our uncorrected Carreau model (Figure 4.10 (a)). As we can see in

Figure 4.12, there is good agreement in results, confirming the validity of our

shift.

Parameter Value
η0 9.086963
η∞ 0.409465
λ 78490
n 0.574426

Table 4.6: Corrected Carreau Model for Porous media flow

Figure 4.12: The validation of our corrected Carreau model by comparing
computed simulation results to the original Carreau model of
best overall results

Sensitivity of results

We have already established in earlier sections how changing the pressure gra-

dient would have an effect on our results. It is worth also analyzing how our

74



results translate to other process parameters, such as change in geometry or rea-

son behind void formation, and change in pe