Inverse modelling of compression moulding of sheet moulding compound using CFD

(1)

2007:168 CIV

M A S T E R ' S T H E S I S

Inverse modelling of compression moulding of sheet moulding compound using CFD

Sofia Ebermark

Luleå University of Technology MSc Programmes in Engineering

Engineering Physics

Department of Applied Physics and Mechanical Engineering Division of Fluid Mechanics

2007:168 CIV - ISSN: 1402-1617 - ISRN: LTU-EX--07/168--SE

(2)

(3)

iii

Preface

This work has been carried out at the Division of Fluid Mechanics at Lule˚a University of Technology in Sweden as part of the Research Trainee Pro- gramme 2006/2007. This thesis constitutes the final part of the Master of Science Programme in Engineering Physics. This work has been part of the FYS-project (F¨orb¨attrade ytor i SMC-produkter) that is lead by Scania and sponsored by VINNOVA and the participating industries.

First, I would like to thank my examiner and supervisor Professor Staffan Lundstr¨om for guidance and support throughout the work and also my as- sistant supervisor Ph.D. Daniel Marjavaara for his valuable help during this work. I would also like to thank the people involved in the FYS-project.

Finally, I would like to thank the Division of Fluid Mechanics for always making me feel part of the group and for a great atmosphere and the Re- search Trainee group of 2006/2007 for making this a great year.

Lule˚a, June 2007.

Sofia Ebermark

(4)

(5)

v

Abstract

Sheet Moulding Compound (SMC) is primary used in compression moulding which is a common manufacturing process for composite materials. When manufacturing SMC, voids are likely to appear on the surfaces of the parts causing large defects after painting. The automotive industry often use SMC in their applications implying that the parts often must have a class A appearance. When defects are discovered on the surfaces of the parts they need to be rejected or repaired which is a costly process. The formation of voids is strongly related to the spatial distribution of the pressure. It is therefore interesting to reveal the pressure distribution during moulding. If the pressure distribution is known critical areas can be localized. The purpose of the work is to investigate if an inverse modelling approach together with Computational Fluid Dynamics (CFD) can be used to predict the pressure distribution during compression moulding of SMC.

The geometry in focus is a circular plate with a radius of 100 mm placed between two heated mould halves. To be able to find the pressure distribution the viscosity as a function of time and spatial coordinate needs to be known since it relates the deformation rate to stresses and hence the pressure. The viscosity of the SMC is known to be dependent on many parameters such as temperature, shear rate, degree of curing, amount of fill-material and fibre orientation. Here, the viscosity is assumed to be a function of temperature only since it is difficult to derive an explicit expression for the viscosity as a function of all the parameters mentioned above. A way forward is instead to do measurements of the pressure during moulding in a simple geometry and at the same time do numerical simulations that are fitted, by adjusting parameters in the viscosity model, to the experimental results by inverse modelling.

The results from the experiments show a rapid increase in pressure initially.

After the pressure has reached its peak value it decreases almost linearly with time. The same behaviour is captured by the viscosity model used in the numerical simulations. By an optimization procedure it is possible to find the parameters in the viscosity model that gives the best fit between the experiments and the simulations. The optimization is performed at one spatial coordinate and the resulting pressure is similar to the one obtained in the experiments. The pressure however differs when compared at an other spatial coordinate. Comparisons between visualizations made previously of the in-mould flow and the simulations show however similar flow behaviour.

In other words, the numerical model serve as a good start but needs further improvements in order to capture the true mould flow behaviour.

(6)

(7)

Chapter 1 Introduction

1.1 Background

Compression moulding is a common manufacturing process for composite materials, for example Sheet Moulding Compound (SMC). During this process voids are likely to appear on the surface as flaws and great efforts must be spent on minimizing the number of them. This is especially important in the automotive industries since their applications often must have a class A appearance. Surface voids can be transported out of the SMC or dissolve into it and these mechanisms are strongly related to the spatial distribution of the pressure. Hence, if the pressure distribution is known, critical areas where voids appear can be discovered.

1.2 Aim and Objective

The purpose of the work is to investigate if inverse modelling together with Computational Fluid Dynamics (CFD) can be used to predict the pressure distribution during compression moulding of SMC. The aim with this work is to implement a numerical model of the compression moulding process and to find a viscosity model that describes the pressure response of the material.

1.3 Previous work

Traditionally the formation of voids were investigated with experiments but since the capacity of computers in the recent years has rapidly increased, full 3D simulations of the mould filling process is possible to perform to get a deeper insight into the physical phenomena governing the process. One issue to be solved however is the modelling of the material i.e. the SMC-material since its dependence on several parameters is somewhat complex. Inverse modelling may be an alternative approach. In previous work done by P.T.

Odenberger [1] experiments and numerical simulations of the compression 1

(10)

moulding process together with inverse modelling in a simple geometry were established. The viscosity of the SMC was allowed to vary as a function of time but was constant throughout the domain at each time.

1.4 Description of project

In this project the same arrangement as in previous work [1] will be used.

A numerical model of the compression of a circular disc between two plates during compression moulding will be implemented. In contrast to previous work the viscosity will change as a function of temperature and the temperature in its turn will vary as a function of time and spatial coordinate.

Hence, the equations governing the heat transfer need to be solved. By inverse modelling the parameters in the viscosity model will be adjusted so that the pressure in the numerical simulations agrees with the experiments.

1.5 Outline of this thesis

First, a short introduction to composites manufacturing is presented to give some insight of the process. Then, the inverse modelling procedure will be explained further and also some theory about how the governing equations are solved numerically. After that, the experimental and numerical setup will be presented followed by the results from the experiments and the inverse modelling. Finally the results will be discussed and some conclusions drawn and also some suggestions for future work will be given.

(11)

Chapter 2 Theory

2.1 Fibre Reinforced Polymer Composites

Fibre reinforced polymer composites are because of their high strength to weight ratio commonly used in the automotive industry to reduce weight and hence fuel consumption. A fibre reinforced polymer composite is a polymer matrix reinforced with fibres. These fibres are usually made from glass, carbon or aramid.

2.2 Compression Moulding

Compression moulding is a high-volume manufacturing method for polymer composites and is suitable for manufacturing complex and high-strength parts. In compression moulding the moulding material is placed in an open heated mould cavity, see figure 2.1. The mould is then closed and a pressure applied to force the material to fill the cavity. Curing takes place and when the mould is opened the part is removed.

2.3 Sheet Moulding Compound (SMC)

Sheet moulding compound (SMC) or sheet moulding composite is a fibre- reinforced polymer material of typically 30% glass fibres primarily used in compression moulding. It is composed basically of four ingredients that is thermosetting resin, fibres, fillers and additives [2]. The automotive industry has the highest level of consumption of SMC. For example some of the exterior parts on the Scania truck in figure 2.2 are made from SMC.

2.3.1 Defects in SMC

The parts made from SMC used in automotive applications such as hood panels must have high quality surfaces. When defects are discovered on

3

(12)

Figure 2.1: Schematic figure of the compression moulding process.

Figure 2.2: Scania truck with exterior parts made from SMC.

(13)

2.4. INVERSE MODELLING 5 the surfaces of the parts they need to be rejected or repaired which is a costly process. The defects originate from small voids on the surfaces of the parts. When the parts are painted, the air entrapped in these voids expands resulting in large bubbles which finally collapse and cause large defects on the surfaces.

2.4 Inverse Modelling

A system can be described by a set of model parameters. Some physical relation should exist between the model parameters and some measurable parameters. By inverse modelling the model parameters are tuned to best fit the measurable parameters obtained in experiments [3]. Inverse modelling may be a helpful tool when dealing with complex materials whose dependence on several parameters is difficult to model.

2.5 Surrogate based optimization

Instead of performing time-consuming simulations for every parameter com- bination in a numerical model, a surrogate based optimization approach can be used. The surrogate model is an approximate model that is intended to mimic the behaviour of the simulation model as closely as possible. The exact function of the simulation code is not assumed to be known (or even understood), solely the input-output behaviour is important. The code thus behaves like a black-box. The aim with the surrogate modelling is to deter- mine a relationship between of a set of design variables x = (x₁, x₂, ..., x_N) and their calculated response y for a set of data points.

The true response of a system is given by:

y = f (x), (2.1)

while the response of the surrogate model is given by:

ˆ

y = g(x), (2.2)

such that:

y = ˆy + ε, (2.3)

where ε is the modelling error [4].

Surrogate modelling can be used to capture the global behaviour of the system but it can also be used to find an optimum. The surrogate based optimization (SBO) process can be divided into four steps;

1. Experimental Design

(14)

Figure 2.3: The surrogate model is a least square approximation as illus- trated here.

2. Numerical simulations at selected locations 3. Construction of surrogate model

4. Model validation

In the first step, the experimental design, the set of data points is chosen and these can be randomly distributed or form a pattern. In the second step, numerical simulations are performed for the chosen set of data points and then a surrogate model is fitted to the calculated data. Finally, the model is used predict the response for given values of the design variables in order to validate it [5].

2.5.1 Response Surface Methodology

One of the most common surrogate models is the polynomial response sur- faces (PRS) model in which polynomials are fitted to the calculated data.

An example of a second-order response surface model in two variables is:

y = β₀+ β₁x₁+ β₂x₂+ β₃x²₁+ β₄x²₂+ β₅x₁x₂+ ε. (2.4) This is a multiple regression model and the parameters β_i are the regression coeffients and the independent variables x_i are the regressor variables. More generally this regression model can be written as:

y_i= β₀+ Xk j=1

β_jx_ij + Xk j=1

β_jjx²_j+^X Xk j<i=2

β_ijx_ix_j+ ε_i (2.5) or in matrix notation:

y = Xβ + ε

(15)

2.6. COMPUTATIONAL METHODS 7 The method of least squares selects the β’s so that the sum of squares of the errors ε is minimized. The least squares estimators of the regression coefficients are denoted b₀, . . . ,b_k.

b =^¡X⁰X^¢⁻¹X⁰y (2.6)

and the polynomial response surface model can be written as [6]:

ˆ

y(x) = x^Tb (2.7)

2.6 Computational Methods

2.6.1 Governing equations

Fluid motion is governed by the conservation of mass, momentum and energy [7]. These conservation laws can be expressed in several forms e.g.

differential and integral form. Here they will be presented in differential form.

The differential form of the mass conservation equation is

∂ρ

∂t +∂ (ρu_i)

∂x_i = 0 (2.8)

where u_i((i = 1, 2, 3) or (i = x, y, z)) are the components of the velocity vector in the direction of the coordinates x_i and ρ is the density of the fluid and t is the time. This equation is known as the continuity equation.

By applying Newton’s second law of motion to an infinitesimal fluid element, the differential form of the conservation of momentum can be derived yielding the Cauchy’s equation of motion

ρDu_i

Dt = ρg_i+∂τ_ij

∂x_j (2.9)

where ρg_i and ^∂τ_∂x^ij

j is the i-component of the body force and the surface force per unit volume respectively. The stress tensor τ_ij is a constitutive equation relating the stress to the deformation and for a Newtonian fluid and it can be written as

τ_ij = − µ

p +2 3µ∇ · u

¶

δ_ij+ 2µe_ij (2.10) in which p denotes the static pressure, µ the dynamic viscosity, δ_ij the Kronecker delta and e_ij the strain rate tensor [8],

(16)

e_ij = 1 2

Ã∂u_i

∂x_j +∂u_j

∂x_i

!

. (2.11)

A substitution of the constitutive equation (eq. 2.10) into the Cauchy’s equation (eq. 2.9) results in the equation of motion for a Newtonian fluid that can be expressed in the following way

ρDu_i

Dt = −∂p

∂x_i + ρg_i+ ∂

∂x_j

·

2µe_ij− 2

3µ (∇ · u) δ_ij

¸

(2.12) being the general form of the Navier-Stokes equation. The viscosity µ in this equation can be a function of the temperature.

The conservation of energy is provided by the first law of thermodynamics.

ρD Dt

µ e +1

2u²_i

¶

= ρg_iu_i+ ∂

∂x_i(τ_iju_i) − ∂q_i

∂x_i (2.13)

in which the internal energy is represented by e and the heat flux per unit area by q_i. For a perfect gas e = C_vT , where C_v is the specific heat at constant volume. If the mechanical energy equation is subtracted from (eq.

2.13) the thermal energy equation is obtained,

ρDe

Dt = −∇ · q − p (∇ · u) + φ. (2.14)

2.7 Computational Fluid Dynamics (CFD)

The set of partial differential equations (PDEs) describing fluid motion, re- ferred to as the Navier-Stokes equations, can only be solved analytically for a few special cases. Although known solutions are useful in helping to under- stand fluid flow, they can rarely be used directly in engineering analysis or design. Traditionally simplifications of the equations and experiments were used to investigate fluid flow. For many types of flow these methods are not sufficient and an alternative or complementary method is needed. With the help of computers the partial differential equations of fluid mechanics can be solved numerically by an iterative procedure. This is known as Compu- tational Fluid Dynamics (CFD).

In order to solve the differential equations numerically they need to be approximated by a system of algebraic equations using a discretization method.

The algebraic equations are solved at discrete locations in space and time in the computational domain [8].

(17)

2.7. COMPUTATIONAL FLUID DYNAMICS (CFD) 9 2.7.1 The Finite Volume Method

The Finite Volume Method (FVM) is used as a discretization method in most CFD codes such as ANSYS-CFX. In this method, the computational domain is subdivided by a grid into a finite number of control volumes (CVs). Val- ues of the variables are calculated at the computational nodes being located at the centre of the control volumes. An integral form of the conservation equations are applied to each control volume and the integrals are approximated using variable values at the control volume surfaces [8]. Variable values at locations other than the computational nodes are approximated in terms of the nodal values using a difference scheme. In order to preserve the accuracy, the order of the integral approximation must be at least the same order as the difference scheme. Many difference schemes are based on series expansion approximations such as the Taylor series. The more terms of the expansion used in the difference scheme, the more accurate the approximation will be. More terms in the expansion will however increase the computational load. The order of the largest term in the truncated part of the series expansion determines the order of the scheme.

2.7.2 Errors and accuracy

Numerical solutions are approximations and they will always contain errors arising from different parts of the process of obtaining the numerical solution.

Three kinds of systematic errors are however always included in a numerical solution [8]:

• Modelling errors

• Discretization errors

• Iteration errors

Modelling errors

The modelling errors are defined as the difference between the actual flow and the exact solution of the mathematical model of the problem. These errors originate from assumptions made when deriving the transport equations for the variables. Turbulent flows, multiphase flows, combustion and other complicated flows involve physical phenomena that are not perfectly described by current scientific theories and therefore modelling errors for these types of flow may become very large. Simplifying the geometry and/or simplifying the boundary conditions may also introduce modelling errors. The modelling errors can only be estimated by comparing numerical solutions, in which the discretization and iteration errors are negligible, with accurate experimental data or data from more accurate models.

(18)

Discretization errors

The discretization errors are defined as the difference between the exact solution of the conservation equations and the exact solution of the corresponding discretized equations. These errors should become small as the grid spacing goes to zero. A way to estimate the discretization error is to compare solutions from grids with different number of cells. If three grids are available an estimate of the discretization error ε^d_h on a grid with grid spacing h can be obtained with:

ε^d_h≈ φ_h− φ_2h

2^p− 1 (2.15)

where φ denotes the numerical solution and p is an exponent defined by:

p = log^³^φ_φ^2h^−φ^4h

h−φ2h

´

log 2 (2.16)

An approximation of the exact solution can then be obtained by adding the error estimate ε^d_h to the solution on the finest grid φ_h. This is method is known as Richardson extrapolation and is valid only when the convergence is monotonic which can be expected only on sufficient fine grids.

Iteration errors

The system of algebraic equations obtained in the discretization process is usually non-linear and therefore needs to be linearized in order to be solved. This linear system is then solved by iterative methods in which a starting guess is successively improved. The iteration errors are defined as the difference between the iterative and the exact solutions of the system of algebraic equations and it decreases as the number of iterations increases.

The residual is a measure of iterative convergence as it relates to whether the equations have been solved.

2.7.3 Multiphase Flow

In multiphase flow more than one phase is present. Each phase may have its own flow field or all phases may share a common flow field. In multiphase flow the fluids are mixed on a macroscopic scale in contrast to multicompo- nent flow where they are mixed on a microscopic scale. There exist a number of models for these types of flows. Two distinct multiphase flow models are the Eulerian-Eulerian multiphase model and the Lagrangian Particel Track- ing multiphase model.

The flow can be either homogeneous or inhomogeneous. In the first, all the fluids share the same velocity fields and other relevant fields such as

(19)

2.7. COMPUTATIONAL FLUID DYNAMICS (CFD) 11 temperature, turbulence and pressure. In the second case separate velocity fields and other relevant fields exist for each fluid. The pressure field is however shared by all fluids. The fluids interact via interphase transfer terms [9].

2.7.4 Fluid Structure Interaction (FSI)

In Fluid Structure Interaction (FSI) simulations the solution fields in fluid and solid domains are coupled. The coupling between the solution fields can be uni-directional, which means that one solution field strongly affects the other fields but is not being affected by the others [9].

2.7.5 Transient simulation

The time dependence of the flow characteristics is either steady state or transient. In transient simulations the flow characteristics change with time and real time information is required to determine the time intervals at which the flow field is calculated [9].

Courant number

The dimensionless Courant number is one of the key parameters in com- putational fluid dynamics and can be used to determine the time step in a transient simulation. The Courant number is defined as the ratio of the time step ∆t to the characteristic convection time, which is the time required for a disturbance to be convected a typical length of a computational cell, ∆x [8]:

c = u∆t

∆x (2.17)

When diffusion is negligible the criterion for stability is:

c < 1 (2.18)

or

∆t < ∆x

u (2.19)

i.e. the Courant number should be smaller than unity.

(20)

(21)

Chapter 3 Methods

3.1 Experimental Setup

To give input to the inverse modelling process measurements are performed during pressing in a Fjellman 310 ton press. Charges are compressed in a plate to plate mould held at a constant temperature of 150˚C. The pressure is measured with two Kistler pressure transducers mounted in the upper mould at two spatial locations that is, at the centre and 37.5 mm away it, P0 and P1 respectively. The mould closing speed is set to 2 mm/s corresponding to real production conditions. The height of the upper mould as a function of time is also recorded. The SMC-charge has a circular diameter of 100 mm and consists of 6 layers of SMC (20 wt% fibre glass content) which gives a total height of 15 mm.

Figure 3.1: Schematic figure of the experimental setup with the position of the pressure transducers. The area inside the dotted line is included in the numerical model.

13

(22)

3.2 Numerical Setup

3.2.1 Geometry

The geometry in focus is a circular plate of radius 200 mm, which symbol- izes the region between two mould halves and it includes both the SMC- charge and the air surrounding it. In order to save computational time axi-symmetric conditions are applied in the θ-direction. The dimensions of the SMC-charge are taken from experiments and correspond to real production conditions. That is, a height of 15 mm and a radius of 50 mm. The geometry is created in ANSYS Workbench 10.0 and in order to facilitate the grid generation process the centre of the plate is not included in the geometry. Instead, it has an inner radius of 1.5mm.

Figure 3.2: The geometry used in the simulations. It has a radius of 200 mm and a height of 15 mm.

3.2.2 Grid generation

The grid is generated in ANSYS ICEM CFD 10.0 using structured hexahe- dral elements. The elements are uniformly distributed in order to preserve a good mesh quality during the entire simulation process. Three grids are created for the grid convergence study consisting of 2k, 4.2k and 8.5k elements respectively. Due to the geometrically simple domain almost regular elements was achieved guaranteeing an almost 90 degree angle between element sides. The aspect ratio which is a measure of the relative length between the sides of the element had a maximum value of 15.7 for the 4.2k mesh.

(23)

3.2. NUMERICAL SETUP 15

Figure 3.3: Part of the mesh at the outlet boundary. The elements are uniformly distributed and there is only one element in the theta-direction.

3.2.3 Physics definition

To obtain the resulting three-dimensional flow field from the fluid structure interaction, a two-phase, transient simulation is performed with a high res- olution advection scheme. The flow field is assumed to be laminar due to the relatively high viscosities. The commercial CFD code CFX-10.0 from ANSYS [9] is used for obtaining the solution. To solve the fluid structure interaction, the mesh displacement equations are activated and the mesh stiffness value is set to 1 to ensure that the mesh displacements are homoge- neously diffused throughout the mesh. Since this is a transient simulation a time step is set and this is done in an adaptive manner in order to prevent the root mean square value of the courant number to exceed 1 anywhere in the the domain.

3.2.4 Domain models

An Eulerian-Eulerian homogeneous multiphase flow model is used in which the fluids share the same velocity field. The two phases present are SMC at a variety of temperatures and air at 25˚C and they are both considered continuous. Values for the SMC-material parameters are presented in table 3.1.

All the parameters except the viscosity are assumed, for simplicity, not to be a function of temperature, which is not true for some of them. Some of the material parameters are actually unique for each SMC-recipe but they serve as an acceptable approximation.

(24)

Table 3.1: Parameter values for the SMC-material Molar mass 7000 g·mol⁻¹

Density 1.9 g·cm⁻³

Specific Heat Capacity

1500 J·kg⁻¹·K⁻¹ Thermal

Conductivity

0.8 W·m⁻¹·K⁻¹ Thermal Ex-

pansivity

1.4·10⁻⁵ K⁻¹

The SMC-charge is placed in the domain by setting its volume fraction to 1 and the air to 0 for r < 50 mm and the other way around for

50 mm < r < 200 mm. The free surface model is used since the two phases are assumed to be separated by a distinct surface. As a heat transfer model, the thermal energy model is used since kinetic effects are assumed to be negligible.

This is a fluid structure interaction (FSI) simulation where a solid wall is displacing a fluid. In this case the upper mould wall is displacing the SMC-charge. The mesh will therefore be deformed according to a preset displacement.

3.2.5 Boundary conditions

Mesh motion boundary conditions are specified on all boundaries in the domain. They are applied on all nodes on the boundary region unlike other boundary conditions, which are set on mesh faces. The upper wall is set to move downwards with a velocity of 2 mm/s, the lower wall is specified as stationary and the other boundaries has an unspecified mesh displacement.

Outlet boundary conditions: At the outlet an opening boundary condition is set with atmospheric pressure and flow direction normal to the boundary.

Wall boundary conditions: The mould surfaces are assumed to be held at constant temperature, 150˚C and therefore the upper and lower wall has isothermal boundary conditions with no slip. The wall at the inner radius is a non slip adiabatic wall. On the two remaining boundaries, symmetry boundary conditions are applied.

(25)

3.3. INVERSE MODELLING AND RHEOLOGY 17

3.3 Inverse Modelling and Rheology

The aim with this thesis is to predict the pressure in the material during moulding. To be able to do that the viscosity as a function of time and spatial coordinate needs to be known since it relates the deformation rate to stresses and hence the pressure [10]. The relation between the deformation rate and stresses is non-linear for the SMC-material, which is the definition of a non-Newtonian fluid. This non-linear relationship will however not be considered here. The viscosity of the SMC is known to be dependent on many parameters such as temperature, shear rate, degree of curing, amount of fill-material and fibre orientation. It is not possible to derive an explicit expression for the viscosity as a function of all these parameters. A way forward is instead to do measurements of the pressure during moulding in a simple geometry and at the same time do numerical simulations that are fitted to the experimental results by inverse modelling.

3.3.1 Viscosity model

In this work the viscosity is assumed to be a function of the temperature in the material according to equation 3.1 below,

µ = Ae^{−B(T −298)} (3.1)

where A represents the value of the viscosity at room temperature and B represent the rate of decrease with temperature. Initial values of A and B are taken from literature [11]. These variables are unique for each SMC- recipe but the chosen values serve as a good start and will be tuned by the inverse modelling process to best fit the SMC used in the experiments.

3.3.2 Optimization

The regression model chosen for this problem is as follows:

y = β₀+ β₁A + β₂B + β₃A²+ β₄B²+ β₅AB + ε (3.2) where A and B are the parameters in the viscosity model. The dependent variable, y, in the regression model is defined as the sum of squares of the difference between the measured and calculated value of the pressure,

y = X32 i=1

(p_exp− p_sim)², (3.3) where p_exp is the pressure obtained in the experiments and p_sim is the corresponding pressure from the simulations (see also figure 3.4). The pressure differences are calculated for every 0.1 second for the first 3 seconds at the

(26)

Figure 3.4: Schematic figure of the objective function y which is defined as the sum of squares of the difference between the measured and calculated value of the pressures, that is ε(t).

centre of the geometry (P0).

The aim is to find the values of A and B that minimizes the objective func- tion y. In other words, find the A and B that gives a pressure that is similar to the one obtained in the experiments.

The parameter A in the viscosity model represents the initial value of the viscosity at T = 298 K (25˚C) and B represents the rate of change with temperature. In other words, a large value of the parameter A will produce a high initial viscosity and hence a large initial value of the pressure. If B is large the viscosity will decrease rapidly with temperature. Values of the parameter A are limited to be between 1 · 10⁵ Pa·s and 4 · 10⁵ Pa·s and the values of B between 0.06 K⁻¹ and 0.14 K⁻¹, according to figure 3.5 below.

These values are based on a few simulations.

(27)

3.3. INVERSE MODELLING AND RHEOLOGY 19

Figure 3.5: Schematic figure of the experimental design with the chosen data points.

(28)

(29)

Chapter 4 Results

4.1 Experimental results

In figure 4.1 the pressure, at the centre of the upper surface of the charge, is presented for four mouldings with the same experimental setup, see section 3.1. In order to prevent preheating of the charge, the upper mould half begins to move downwards as soon as the charge is placed on the lower mould half. Time t = 0 is defined as the time when the upper mould half hits the SMC-charge and the pressing begins, which is in reality the point where the pressure rises. When this happens the upper mould is a distance of about 15 mm above its lowest position. This distance corresponds to the height of the SMC-charge. Three of the pressings give similar results for the pressure at P0 while the second pressing, test B, differs from the rest.

Figure 4.2 shows the pressure a distance 37.5 mm away from the centre for the same four pressings.

Figure 4.3 presents test C with the displacement of the upper mould in- cluded. It can be seen that at t = 0 when the pressing is assumed to begin the upper mould is a distance of 15 mm above the lower mould, which is the height of the SMC-charge.

4.2 Results CFD-simulations

As a second step in the response surface optimization procedure numerical simulations are performed using the commercial CFD software CFX-10.0 for a variety of conditions representing selected locations in the design space.

Since the simulations proved to be time-consuming for the finest grid, the calculations are performed on the coarsest grid. Afterwards a grid convergence study is performed to estimate the grid error.

21

(30)

Figure 4.1: Pressure at the centre of the upper surface of the SMC-charge.

Figure 4.2: Pressure at a distance 37.5 mm away from the centre of the charge.

(31)

4.2. RESULTS CFD-SIMULATIONS 23

Figure 4.3: The pressure at P0 and P1 and the displacement of the upper mould for test C.

4.2.1 Outer design

First, outer design points are chosen with values of A and B according to figure 3.5 in section 3.3.2. The corresponding y-values and response surface ˆ

y(x) is plotted in figure 4.4. This response surface experiences no minimum but the pressure curves in figure 4.5 shows that A=2.5 · 10⁵ Pa·s gives a pressure at P0 close to the measured value. The pressures from test C in the experiments are chosen to represent the experimental values of the pressures at position P0.

4.2.2 Inner design points included

From the results in the previous section new values of A are selected between 2 · 10⁵ Pa·s and 3 · 10⁵ Pa·s while the values of B remains between 0.1 K⁻¹ and 0.14 K⁻¹. Two additional design points are also included in order to further improve the response surface model (see figure 4.6). When all the data points are included in the model the response surface in figure 4.7 is obtained. An optimum value of the response surface ˆy(x), that is a minimum value, is found for A = 2.509 · 10⁵ Pa·s and B = 0.151 K⁻¹. A simulation for this optimum is performed and the corresponding response, y, is added to the response surface model.

(32)

Figure 4.4: Response surface plot for the outer data points.

Figure 4.5: Pressure as a function of time at position P0 for values of A and B according to the outer data points.

(33)

Figure 4.6: Modified experimental design with outer, inner and additional data points.

Figure 4.7: Response surface plot when all the data points are included.

(34)

Figure 4.8: Contour plot when all data points are included in the model.

Optimum found for A = 2.509 · 10⁵ Pa·s and B = 0.151 K⁻¹

Figure 4.9: Pressures for the design points with the lowest calculated value of y which means they are closest to the experimental values of the pressure.

(35)

4.2. RESULTS CFD-SIMULATIONS 27 4.2.3 Optimum for t > 0.5 s

The differences in pressure, that is (p_exp− p_sim)², between t = 0 s and t = 0.5 s will have large influence on the overall sum, while the differences for t > 0.5 s are more important for the visual fit. Therefore the summation will be carried out for t > 0.5 s instead. Also, to find a good approximation of the response surface near its optimum, only the points closest to the point with the smallest calculated value of y are included in the model. A new optimum value is then found for A = 2.521 · 10⁵ Pa·s and B = 0.141 K⁻¹. The corresponding calculated response is plotted in figure 4.10. The contour plot of the response surface is presented in figure 4.11.

Figure 4.10: Response surface plot near its optimum value. The point marked with a blue box in the surface plot represents the lowest calculated value of y. A = 2.5 · 10⁵ Pa·s and B = 0.14 K⁻¹ for this point.

4.2.4 Validation of optimum

In the final step of the optimization process the surrogate model, ˆy(x), is checked to see whether it can be used to predict the response for given values of the design variables. To validate the optimum and hence the response surface model, a simulation is performed and the calculated and predicted values are presented in figure 4.12. The red box in the figure is the calculated

(36)

Figure 4.11: Contour plot of the response surface near its optimum. The point in the middle represents the optimum.

value. It proves to be very close to the value predicted by the response surface (the green box), see also table 4.1. The values of A and B that produces the lowest calculated value of y gives the pressure that is presented, together with the measured value of the pressure, in figure 4.13 .

Table 4.1: Results near the optimum.

A (P a · s) B (K⁻¹) y (P a²) Predicted value 2.521 · 10⁵ 0.141 5.36 · 10¹⁰ Calculated value 2.521 · 10⁵ 0.141 6.14 · 10¹⁰ Lowest calculated value 2.5 · 10⁵ 0.14 5.48 · 10¹⁰

4.2.5 Model Validation

To validate if the viscosity model obtained in the optimization process is valid at other spatial locations, the pressure at a distance of 37.5 mm away from the centre is compared between the simulations and the experiments.

The result is presented in figure 4.14 below. As can be seen in the figure the pressure from the simulations differ quite a lot from the measured value of the pressure.

(37)

Figure 4.12: The calculated value of the objective function is marked by a green box in the response surface plot and the corresponding predicted value is marked with a blue box.

Figure 4.13: Pressure at P0 for the experiment and the simulation with the lowest calculated value of y.

(38)

Figure 4.14: Pressure at P1 for the experiment and the simulation with the lowest calculated value of y.

4.2.6 Results for different grids

A grid convergence study is performed to estimate the discretization error.

Three grids with 2k, 4.2k and 8.5k elements are used. The values of A and B used in the viscosity model are the values obtained for optimum. Figure 4.15 below show the result. The rise in pressure during the last 0.5 s may be an indication of a diverging solution.

4.2.7 Visualisation of the simulations

To visualise the evolution of the flow front for the first three seconds of the simulation the volume fraction of the SMC is plotted in figure 4.16. As can be seen in the pictures the interface between the two phases becomes more and more diffuse with time due to numerical diffusion.

(39)

Figure 4.15: Pressure for different grids. For the extrapolated value the grid error is added to the solution for the finest grid.

(a) t = 0 s (b) t = 1 s

(c) t = 2 s (d) t = 3 s

Figure 4.16: Pictures from the simulations showing the volume fraction of the SMC for the first three seconds of the pressing.

(40)

(41)

Chapter 5 Discussion

5.1 Assumptions and simplifications

As most mathematical models, the numerical model used in this thesis also involves some assumptions and simplifications. These simplifications limit the possibility to solve the real flow field and some of them are presented here:

Axi-symmetric flow: This is probably not true for most mouldings since the SMC-material becomes anisotropic when the initially randomly distributed glass fibres are oriented by the flow. This orientation will probably produce local differences in pressure that is not captured by the model. This simplification is made in order to simplify the model and hence reduce the simulation time and the overall pressure distribution will probably be captured.

Newtonian flow: The SMC-material is known to have a shear-thinning be- haviour, which means that the viscosity is dependent on the shear rate. This non-Newtonian behaviour is not taken into account in this model.

Material parameters: The material parameters are constant (except the vis- cosity) and are assumed not to be a function of temperature. This is true if the temperature gradients are small, but is not true in this case where the temperature goes from 150˚C at the surface to 25˚C at the centre. Also, the values of the constants are unique for every SMC-recipe. The values for the SMC used here are taken from literature not from actual measurements on the material.

One phase only: The real SMC-material consists of many phases such as resin, fibres and fill-material. Here, the material is assumed continuous consisting of only one phase, the SMC-resin. The charge consists of 6-7 layers

33

(42)

of SMC and air is actually entrapped between these layers while in the numerical model the charge is assumed to consist entirely of SMC.

Symmetric flow: In-real case moulding scenarios the SMC-charge experi- ences some degree of pre-heating before the pressing begins. This occurs since charge is placed on the lower mould half before it comes into contact with the upper mould half.

5.2 Conclusions

The inverse modelling proved to be a successful method in the process of fitting the numerical model to the experimental results, at least for the pressure at the centre of the disc. However, in order to get a viscosity model that predicts pressure in the entire disc, further improvements needs to be done.

The viscosity model seems to capture the overall behaviour of the pressure.

That is, a relatively high pressure initially which is then decreasing almost linearly with time. By changing the constants in the viscosity model the magnitude of the initial pressure and the rate of change with time can be adjusted. The simulations show an immediate increase in pressure while the pressure in the experiments reaches its peak value after 0.5s. This transient effect is probably a result of the initial air entrapment between the SMC- sheets. It takes some time before these layers are squeezed together in the pressings while in the simulations the charge is modelled as a homogenous material.

When the data points in the response surface model were too far from the optimum value, the pressure response could not be approximated with the second order response surface model. However, near its optimum this model served as a good approximation and the pressure response showed a second order behaviour.

Because compression moulding of SMC is relatively fast process, only the material closest to the mould wall will experience an increase in temperature and hence a decrease in viscosity. The simulations show that this reduction in viscosity causes a secondary flow near the mould surface. This secondary flow can only be captured if the elements closest to the wall are small enough.

This is probably the reason why the pressure is so different between the different grids used and is probably also one reason why the optimization on the coarsest grid gave poor agreement for the pressure away from the centre.

(43)

5.3. FUTURE WORK 35 Poor agreement was achieved when comparing the pressure at another spatial point namely P1 between the experiments and the simulations. This indicates that the model has to be improved before it can be used to predict the pressure distribution in the whole geometry. Comparisons between visualisations made previously of the in-mould flow [1] and the simulations however show similar flow behaviour. In other words, the numerical model serve as a good start to predict the pressure field but need further improvements in order to capture the true mould flow behaviour.

5.3 Future Work

In order to find a viscosity model that predicts the pressure distribution in the entire disc some improvements of the numerical model and the optimization needs to be done. Some suggestions of improvements are presented here.

The viscosity model can be improved by using more than two adjustable parameters. For example, different exponential functions can be used in different temperature intervals as a drastic change in viscosity is expected at a certain temperature. To improve the numerical model even further a model for the glass fibres can be included.

It is also important to have the correct values of the material parameters governing the heat transfer since the viscosity is strongly dependent on the temperature in the material. An inverse modelling approach can be used for finding these parameters too.

The optimization should be performed on a grid with small enough elements close to the wall to capture the effects from the decrease in viscosity. Another alternative is to include the error, obtained by using a coarser grid, in the objective function y:

y =^X(p_exp(t) − (p_sim(t) − ε_grid(t)))². (5.1) In order to obtain a better prediction of the pressure over the entire disc the optimization should be performed at different spatial locations at the same time. Then, the overall behaviour of the pressure is more likely captured and a better prediction at locations other than those used in the optimization can be obtained. To improve the optimization even more results from experiments and simulations with different mould closing speeds can be included. It is aslo possible to optimize for different SMC-materials.

When a good enough prediction for the pressure distribution is achieved voids can be included in the simulations to see how they move with the

(44)

flow. Then areas with a large concentration of voids can be localized. When this is known the formation of surface voids when manufacturing SMC can be counteracted. The numerical model can then be used to predict critical areas where voids are likely to appear when manufacturing SMC-parts.

(45)

Bibliography

[1] Odenberger P.T. Moulding of SMC: visualisation and inverse modelling.

Licentiate thesis 2005:34, Lule˚a University of Technology, Sweden.

[2] Kia H.G. Sheet Moulding Compounds: Science and Technology. Hanser, Munich, 1993.

[3] Tarantola A. Inverse Problem Theory. Elsevier, Amsterdam, 1987.

[4] Marjavaara B.D. CFD Driven Optimization of Hydraulic Turbine Draft Tubes using Surrogate Models. Phd thesis 2006:41, Lule˚a University of Technology, Sweden.

[5] Queipo N.V et al. Surrogate-based analysis and optimization. Progress in Aerospace Sciences, 41(1):1–28, 2005.

[6] Myers R.H. and Montgomery D.C. Response Surface Methodology. Wi- ley, New York, 2nd edition, 2002.

[7] Kundu P.K and Cohen I.M. Fluid Mechanics. Academic Press, San Diego, 3rd edition, 2004.

[8] Ferziger J.H and Peric M. Computational Methods for Fluid Dynamics.

Springer-Verlag, Berlin, Heidelberg, 3rd edition, 2002.

[10] Barnes H.A. Hutton J.F. and K. Walters. An introduction to Rheology.

Elsevier, Amsterdam, 1989.

[11] Lee C.C and Tucker C.L. Flow and heat transfer in compression mold filling. Journal of Non-Newtonian Fluid Mechanics, 24:245–264, 1987.

37

(46)

Inverse modelling of compression moulding of sheet moulding compound using CFD

M A S T E R ' S T H E S I S