Evolution of mass distribution in walls of rigid polyurethane foams

Accepted Manuscript Evolution of mass distribution in walls of rigid polyurethane foams Pavel Ferkl, Iveta Kršková, Juraj Kosek PII: DOI: Reference:

S0009-2509(17)30635-8 https://doi.org/10.1016/j.ces.2017.10.024 CES 13856

To appear in:

Chemical Engineering Science

Received Date: Accepted Date:

27 February 2017 19 October 2017

Please cite this article as: P. Ferkl, I. Kršková, J. Kosek, Evolution of mass distribution in walls of rigid polyurethane foams, Chemical Engineering Science (2017), doi: https://doi.org/10.1016/j.ces.2017.10.024

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Evolution of mass distribution in walls of rigid polyurethane foams Pavel Ferkl, Iveta Krˇskov´a , Juraj Kosek∗ Department of Chemical Engineering, University of Chemistry and Technology, Prague, Czechia

Abstract In this work, we describe a theoretical model for the simulation of reactive foaming of polyurethane (PU) consisting of three parts – reaction kinetics, foam expansion and wall evolution, which are coupled together. The advantage of this approach is that it provides comprehensive details about the development of PU foam – from the evolution of temperature and foam density to morphology features like bubble size, wall thickness and strut shape at the same time. The performance of the model is evaluated by analysing the model predictions for realistic process conditions and comparing to experimental data. It is demonstrated that the increased viscosity of the reaction mixture will lead to foams with thicker walls. The presented model can serve as a useful tool for the optimization of PU foaming process. Keywords: Foams, Foaming process, Mathematical modelling, Foam morphology, Wall thickness, Struts

1

1. Introduction

2

Owing to the heat insulation and mechanical properties, polyurethane (PU)

3

foams are widely used in construction, automotive and other industries. The ∗ Corresponding

author Email address: [email protected] (Juraj Kosek)

Preprint submitted to Elsevier

October 20, 2017

4

final application properties of PU foam are heavily dependent on its morphology

5

– density, size of cells, thickness of walls, etc. Thus, it is very important to have

6

precise control over the foaming process to be able to create foam with desirable

7

properties.

8

Polyurethane foams are created by the so-called reaction foaming, during

9

which the polymerization and the foam expansion take place simultaneously.

10

Overall, the process is very complex. The progress of reaction and evaporation

11

of blowing agents dramatically changes the physical properties like the viscosity

12

of reaction mixture and solubility of blowing agents. Meanwhile the morphology

13

of the foam evolves from a system of small isolated spherical bubbles to a low

14

density foam, where the bubbles are polyhedral and most of the condensed phase

15

is located inside walls and struts. Thus, the research papers often focus only on

16

one isolated aspect of the whole process.

17

The reaction kinetics is usually simplified in terms of two overall reactions,

18

i.e., the gelling reaction and the blowing reaction (Baser and Khakhar, 1994a,b).

19

The viscosity of reaction mixture is often described as a function of temperature

20

and conversion by empirical Castro-Macosko model (Castro and Macosko, 1982).

21

These models can be incorporated to CFD mould-filling simulations (Seo and

22

Youn, 2005; Bikard et al., 2007; Geier et al., 2009; Samkhaniani et al., 2013),

23

which can predict the evolution of temperature and foam density. However,

24

they usually lack the information about the bubble size distribution.

25

The evolution of bubble size can be calculated by one of two main ap-

26

proaches. First, using the assumption that the process is entirely controlled

27

by chemical reaction. The concentration of dissolved blowing agent is assumed

28

to be always equal to equilibrium concentration and the rest of the blowing

29

agent is immediately evaporated into bubbles. And second, the more general

30

bubble–shell model (Kim and Youn, 2000), which includes also the effect of

3

31

diffusion limitation. Recently, Karimi and Marchisio (2015) proposed a model,

32

which can predict the evolution of foam density, temperature and bubble size

33

distribution using the population balance equation and successfully coupled it

34

with bubble–shell model in Ferkl et al. (2016). However, even this model lacks

35

additional important aspects of foam morphology like the thickness of walls and

36

content of polymer in struts. This information is critical for the estimation of

37

mechanical (Chen et al., 2015) or heat insulating (Ferkl et al., 2017) properties

38

of the final foam.

39

The prediction of wall thickness in foams is often limited to static liquid

40

foams (Bhakta and Khilar, 1991) and based on Reynolds equation, which pro-

41

vides only average film thickness. Harikrishnan and Khakhar (2009) extended

42

Reynolds equation by a simple geometrical relation to include also the important

43

effect of bubble expansion. Schwartz and Roy (2003) developed a model, which

44

simulates the simultaneous wall drainage and wall stretching and predicts the

45

evolution of wall thickness profile. The model takes advantage of the symmetry

46

of the computational domain created by three bubbles of equal size. Implement-

47

ing the effect of bubble size distribution requires a spatially three-dimensional

48

simulation of multiple growing bubbles. The models based on diffuse interface

49

approach cannot precisely simulate the wall thickness (K¨orner et al., 2005) and

50

sharp interface models require complex adaptive re-meshing schemes (Mitrias

51

et al., 2017). To the best of our knowledge no sharp interface model for foam

52

morphology evolution was published in open literature yet.

53

The main goal of this work is to extend the work of Ferkl et al. (2016) and

54

Schwartz and Roy (2003) and to develop a comprehensive model that will be

55

able to predict the evolution of temperature, foam density, bubble size and wall

56

thickness during PU foaming. This will provide deeper insight into the influence

57

of reaction kinetics and viscosity on the final foam morphology and enable the

4

58

optimization of the foaming process.

59

2. Governing equations

60

There are three main processes happening during the PU foaming: (a)

61

the polymerization producing polyurethane, (b) the foam expansion driven by

62

growth of bubbles and (c) the formation of foam morphology represented on the

63

bubble-scale by the distribution of polymer between walls and struts. For the

64

sake of clarity, we dedicate a section to describe the main equations for each

65

process. However, we must note that the processes are highly interconnected

66

and only together they provide the full picture.

67

2.1. Kinetic model

68

The main purpose of the kinetic model is to calculate the production rate

69

of the carbon dioxide and to provide the evolution of temperature. However,

70

changes of other physical properties like the viscosity of the reaction mixture

71

or the solubility of the blowing agents are also connected to the results of the

72

kinetic model.

73

The detailed polymerization scheme is quite complex. Thus, it is a com-

74

mon practice to simplify it to just two global reactions, which are often called

75

“gelling” and “blowing”. The gelling reaction between isocyanate and polyol,

76

which creates urethane bonds, can be written as:

77

78

R1 – NCO + R2 – OH −−→ R1 – NH – CO – O – R2 isocyanate

79

80

polyol

urethane bond

The blowing reaction between isocyanate and water, which creates urea bonds, can be written as:

81

82

2 R – NCO + H2 O −−→ RNH – CO – NHR + CO2 isocyanate

urea bond

water

5

83

Baser and Khakhar (1994b,a) showed that the gelling reaction is of the sec-

84

ond order and the blowing reaction is of the first order. Taking into account the

85

concentrations of the reactive end groups, it can be shown that the conversion

86

of polyol groups Xp will change according to:   dXp Ep = Ap exp − (1 − Xp )(cc,0 − 2cw,0 Xw − cp,0 Xp ), dt Rg T

87

(1)

and the conversion of water molecules Xw will evolve according to:   dXw Ew = Aw exp − (1 − Xw ), dt Rg T

(2)

88

where t is the time, Ap and Aw are the pre-exponential factors, Ep and Ew are

89

the activation energies, Rg is the gas constant, T is the temperature and cc,0 ,

90

cp,0 and cw,0 are the initial concentrations of isocyanate groups, polyol groups

91

and water, respectively.

92

The enthalpy balance can be also simplified. At the start of the foaming the

93

foam temperature differs only slightly from the room temperature. And at the

94

end of the foaming the foam is a very good heat insulator. Both of these facts

95

cause that the foam loses only minimal amount of heat to the outside. Thus,

96

the adiabatic enthalpy balance can be written as: N

dT (−∆Hp ) cp,0 dXp (−∆Hw ) cw,0 dXw X (−∆Hv,i ) dwi = + + , dt ρrm cp,f dt ρrm cp,f dt cp,f dt i=1

(3)

97

where ∆Hp and ∆Hw are the reaction enthalpies of the gelling and blowing

98

reactions, ρrm is the density of the liquid mixture undergoing polymerization,

99

cp,f is the thermal capacity of the foam, N is the number of blowing agents,

100

∆Hv,i is the heat of evaporation for the ith blowing agent and wi is the mass

101

fraction for the ith blowing agent in the gas phase with respect to the foam.

6

102

2.2. Foam expansion model

103

The main objective of the foam expansion model is to predict the evolution

104

of bubble size and foam density during foaming. This process can be simulated

105

by the bubble–shell model originally developed by Amon and Denson (1984).

106

The model was successfully applied to expansion of polyurethane (Kim and

107

Youn, 2000), polystyrene (Tuladhar and Mackley, 2004) and metal (Sahu et al.,

108

2014) foams.

109

The model uses several simplifying assumptions: (i) the bubbles are spher-

110

ically symmetric, (ii) each bubble is surrounded by an effective shell of the

111

reaction mixture, whose volume is determined by the number of bubbles in the

112

initial mixture (see Fig. 1), (iii) there is no transport of blowing agents or reac-

113

tants between the shells, and (iv) there is a negligible coalescence of gas bubbles.

bubble S(t)

Rb(t)

shell

Figure 1: Illustration of an isolated bubble as considered in the bubble–shell model. 114

115

In case of PU foaming, the initial reaction mixture is vigorously mixed for

116

several seconds. This results in small air bubbles being whipped into the liquid.

117

Thus, the system never reaches large enough supersaturation, which would lead

118

to nucleation of additional bubbles. The foam expansion model simulates the

119

process only after this initial mixing step.

120

The bubble growth can be mathematically described through coupled mo-

121

mentum and mass transfer equations. The dynamics of the growth of a spherical 7

122

123

bubble will result in a purely radial velocity field, and thus the momentum balance can be simplified according to Feng and Bertelo (2004): "  2 # N X d2 Rb 3 dRb 2γ 4µrm dRb pi + pair − prm = ρrm Rb + + + , 2 dt 2 dt Rb Rb dt i=1

(4)

124

where pi is the partial pressure of the ith blowing agent in the bubble, prm is the

125

pressure in the reaction mixture, Rb is the actual bubble radius, γ is the surface

126

tension and µrm is the viscosity of the reaction mixture. The terms on the right

127

hand side represent inertial, surface tension and viscous forces, respectively.

128

Partial pressure of air in the bubble pair can be calculated as follows: pair

129

T = pair,0 T0



3 (5)

where Rb,0 is the initial radius and pair,0 can be calculated as follows: pair,0 = prm +

2γ Rb,0

(6)

The material balance for the ith blowing agent in the bubble can be written

130

131

Rb,0 Rb

as: d dt



pi Rb3 Rg T

 =

3Di Rb2

∂ci , ∂r r=Rb

(7)

132

where Di is the diffusion coefficient of the blowing agent in the reaction mixture,

133

ci is the molar concentration of the blowing agent in the reaction mixture and r

134

is the spatial coordinate. Here we assume that the species are perfectly mixed

135

inside the bubble, and thus the resistance to mass transfer is entirely on the

136

side of the reaction mixture.

137

138

Finally, the mass balance for the ith blowing agent in the reaction mixture can be written as: ∂ci R2 dRb ∂ci Di ∂ + 2b = 2 ∂t r dt ∂r r ∂r

8

 r

2 ∂ci

∂r

 + ri ,

(8)

139

where ri is the reaction source term, which can be expressed as:    cw,0 dXw if i = CO2 dt ri =   0 if i 6= CO

(9)

2

140

We assume that the concentration at the bubble–shell interface is given by

141

the Henry’s law and that the blowing agent is not transported across the outer

142

shell boundary. Thus, the boundary conditions for eq. (8) are set according to:

143

ci |r=Rb = Hi pi ,

(10)

∂ci = 0, ∂r r=S

(11)

and

144

where Hi is the Henry constant and S is the outer radius of the shell. The size

145

of the shell S is a function of time, but it can be directly calculated from the

146

bubble radius and the initial bubble and shell sizes as: S=

147

q 3 3 , Rb3 + S03 − Rb,0

and the foam density ρf can be determined at any time through:  ρf =

148

149

(12)

R3 1 − 3b S

 ρrm .

(13)

Arefmanesh and Advani (1991) showed that it is advantageous to introduce the following substitution: y = r3 − Rb3 ,

(14)

150

which transforms the time varying domain size [Rb , S] to a constant one [0, S03 −

151

3 Rb,0 ]. Here S0 is the initial outer radius of the shell, which is related to the

152

number density of bubbles in the initial mixture, nb,0 through the following

153

expression:  S0 =

3 4πnb,0 9

1/3 .

(15)

154

2.3. Wall evolution model

155

At the start of the foaming, bubbles are far apart and isolated from each

156

other. However, as they grow they start to influence each other. When two

157

bubbles get close to each other, they deform and create a circular film between

158

themselves. This film will extend as the bubbles continue to grow. This is illustrated in Fig. 2.

(a)

(b)

Figure 2: Schematic illustration of (a) two spherical isolated bubbles and (b) formation of a wall between two close bubbles. 159

160

As the close packing of spherical bubbles is advancing during foaming, more

161

than two bubbles get into contact throughout the time. The interaction of two

162

bubbles creates a wall, the interaction of three bubbles creates a strut and the

163

interaction of four bubbles creates a node. The Wall evolution model presented

164

in this paper considers walls and struts and the dynamics of mass redistribution

165

between them. However, it neglects the nodes. The computational domain used

166

in this work is shown in Fig. 3.

167

The distribution of reaction mixture between the wall and the strut is influ-

168

enced by two factors. First, the capillary forces, which are given by the surface

169

tension and the curvature of the liquid–gas interface, decrease the pressure in

170

the strut. This causes the flow of reaction mixture from the wall to the strut

171

(wall drainage effect). And second, the growth of the bubbles stretches the wall

172

and causes the flow of reaction mixture from the strut to the wall. 10

π/6

Rb

symmetry lines

wall

computational domain

gas cell

Rd strut

Figure 3: Schematic illustration of the computational domain for the Wall evolution model.

173

The mathematical model is inspired by the work of Schwartz and Roy (2003),

174

who used Cartesian instead of cylindrical coordinates, incorporated Marangoni

175

forces, used more simplified description of capillary forces, extended the com-

176

putational domain differently, and thus employed different set of boundary con-

177

ditions.

178

The evolution of bubble radius Rb in time is calculated using the Foam

179

expansion model presented in Section 2.2. The domain size Rd (shown in Fig. 3)

180

is then defined as: π
Rd = Rb + h|r=0 − h|r=Rd tan , 6

(16)

181

where h is the film half thickness and r = 0 is the centre of the film. Thus,

182

h|r=0 is half of the initial bubble separation distance.

183

The evolution of wall thickness profile can be calculated using the lubrication

11

184

theory as: (Joye et al., 1992) ∂h 1 1 ∂ = ∂t 3µrm r ∂r



3



rh

∂prm ∂r

 ,

(17)

185

where prm is the pressure in the reaction mixture and µrm is the viscosity of the

186

reaction mixture, which is also a function of time.

187

The pressure is calculated as: pg − prm = γ

1 ∂ (r sin(α)) + Π, r ∂r

(18)

188

where pg is the pressure in gas, γ is the surface tension, Π is the disjoining

189

pressure and sin(α) is defined as: ∂h sin(α) = ∂r

 1+

∂h ∂r

2 !−1/2 .

(19)

190

The disjoining pressure, which incorporates the effect of van der Waals at-

191

tractive forces, electrostatic repulsion of ionic surfactants and steric repulsion

192

of surface molecules, is widely used in the literature as the explanation of the

193

drainage slowdown in liquid soap films. The disjoining pressure isotherm influ-

194

ences the stability of the film (Bhakta and Ruckenstein, 1997). It can have sev-

195

eral maxima in polymer systems due to supramolecular forces (Bergeron, 1999).

196

In PU systems, cell opening (controlled by the overall polymerization/foaming

197

recipe) was linked to separation of urea domains (Yasunaga et al., 1996; Zhang

198

et al., 1999). First, small pinholes are created in the dimple region, where the

199

wall is thinnest, which then leads to partially or fully open cells depending on

200

viscosity. The wall rupture onset depends on critical film thickness (Manev and

201

Nguyen, 2005), which in our system depends on surfactant and type of reactants.

202

The film rupture is a fast process. Thus, the cell opening can be incorporated

203

into the model by checking at each step if the wall is thicker than the critical

204

film thickness at all places.

12

205

To the best of our knowledge, in polyurethane systems, it was used only by

206

Schwartz and Roy (2003). They suggested the following model for the disjoining

207

pressure " Π = B1

208

209

210

213

N −1

 −

h∗ h

M −1 # 

 h∗ −C , h

(20)

where N = 4, M = 3, C = 0.5, B1 = 3 Pa and h∗ = 10 nm. We need four boundary conditions for the system of eqs. (17) and (19). At the center of the film we use no flow boundary conditions ∂h = 0, ∂r r=0 ∂p = 0. ∂r r=0

211

212

h∗ h

(21)

(22)

At the other domain boundary we specify the inclination of the film and the flux Q between the wall and the strut

214

π ∂h = tan , ∂r r=Rd 6

(23)

∂p = Q. ∂r r=Rd

(24)

215

The flux Q in eq. (24) is not a priory known, instead we use Newton’s method

216

at each time step to calculate Q with the condition that the overall mass in strut

217

and wall is constant.

218

When the porosity of foam increases to porosity of close packed spheres

219

it signifies that bubbles come into contact. At this time the Wall evolution

220

simulation starts. We assume that at this time the bubbles are separated by

221

small distance given by model parameter 2h0 and that the bubbles are still

222

spherical at this point.

223

Overall, the mathematical model consists of (N +1) partial differential equa-

224

tions (PDE) – for ci and h (eqs. (8) and (17)) and (N + 4) ordinary differential

225

equations (ODE) – for Xp , Xw , T , R and pi (eqs. (1) to (4) and (7)). Here

13

226

N is the number of blowing agents considered in simulation. The PDEs are

227

discretized in space using the Finite volume method. The resulting system of

228

ODEs is integrated in time using the ODEPACK library (Hindmarsh, 1980).

229

3. Material properties

230

In this paper, we consider water blown PU foams, i.e., water reacts with

231

isocyanates and produces CO2 according to reaction introduced in Section 2.1.

232

Foaming conditions and material properties of this system were taken from lit-

233

erature (Baser and Khakhar, 1994b; Ferkl et al., 2016) and they are summarized in Table 1. Table 1: Summary of operating conditions and material properties.

Property

value

pamb (Pa) ρrm (kg m−3 ) cp,pol (J kg−1 K−1 ) cp,CO2 (J kg−1 K−1 ) MCO2 (kg mol−1 ) DCO2 (m2 s−1 ) HCO2 (mol m−3 Pa−1 ) γ (N m−1 )

1.0 × 105 1100 1800 837 0.044 4.4 × 10−10 1.1 × 10−4 0.025

234

235

The viscosity of the reacting mixture is not constant during the PU foaming

236

process. It decreases with increasing temperature, increases with the conver-

237

sion of polyols and approaches infinity when the gel point is reached. This

238

behaviour is traditionally described by the following semi-empirical Newtonian

239

model (Castro and Macosko, 1982):  µrm = A exp

E Rg T



Xp,g Xp,g − Xp

B+CXp (25)

240

where A = 1.6 × 10−7 Pa s, E = 44.9 × 103 J mol−1 , B = 1.29, C = 1.86 are

241

constants determined from experiments (Lee and Kim, 1991), and Xp,g is the

242

conversion of polyols at the gel point, which is taken in this work equal to 0.6. 14

243

Kinetic parameters of gelling and blowing reactions were published by Baser and Khakhar (1994b) and are summarized in Table 2. Table 2: Summary of kinetic parameters.

Ap (m3 mol−1 s−1 )

Ep (J mol−1 )

−∆Hp (J mol−1 )

Aw (s−1 )

Ew (J mol−1 )

−∆Hw (J mol−1 )

1.735

4.04 × 104

7.07 × 104

1.39 × 103

3.27 × 104

8.60 × 104

244

245

4. Results and discussion

246

In this study we focused on the evolution of walls and struts for water blown

247

PU foams. The recipe is summarized in Table 3 and Fig. 4 shows the evolution of

248

important properties. It can be seen that the predicted evolution of temperature

249

and foam density is in a good agreement with experiments (Baser and Khakhar,

250

1994b). Unfortunately, experimental measurements of evolution of bubble size

251

and wall thickness are still unavailable. They can be determined only from the

252

final foam morphology with the help of image analysis tools and SEM and micro-

253

CT measurements (Nistor et al., 2016). The moment when the foam density

254

decreases to 300 kg m−3 is used as the starting time for the drainage simulation.

255

In this particular case, it was at 47 s.

256

The evolution of reaction mixture viscosity is shown in Fig. 4e. In the begin-

257

ning, the viscosity decreases, because the increase in temperature outweighs the Table 3: Recipe for the studied foam, initial values of parameters.

Name

Value −3

Concentration of isocyanate ci,0 (mol m Concentration of polyol cp,0 (mol m−3 ) Concentration of water cw,0 (mol m−3 ) Temperature T0 (K) Bubble radius Rb,0 (µm) Bubble number density nb,0 (m−3 ) Initial separation of bubbles 2h0 (µm)

15

)

4363 3252 555 305 10.0 1.4 × 1011 5.0

258

increase in polymer concentration. However, as we get close to the gel point,

259

the viscosity quickly starts to increase. The evolution of wall thickness profile

260

is shown in Fig. 4f. It can be seen that the reaction mixture is drained rapidly

261

from the center of the film, which leads to the formation of the wall. After some

262

time, the thickness at the centre of the wall stabilizes, but the wall is further

263

stretched by the bubble growth. The wall develops the so-called dimple shape,

264

where the wall is thinnest near the strut. This dimple shape was experimen-

265

tally observed in PU foams by Zhang et al. (1999). Figure 5 shows the same

266

results, but with the view of the whole strut reconstructed from the computa-

267

tional domain. It can be seen that the strut shape is qualitatively comparable to

268

that of the real foam (see Fig. 6). Moreover, the predicted thickness of walls is

269

around 1 µm and this corresponds to experimental findings (Ahern et al., 2005;

270

Pardo-Alonso et al., 2013).

271

Figure 7 shows the predicted influence of the initial water content on the

272

development of foam morphology. The lines represent 1, 2 and 3 % mass fraction

273

of water with respect to polyol (Baser and Khakhar, 1994b). It can be seen

274

that higher amount of water leads to faster reaction generating CO2 , and thus

275

to faster foam density reduction, and bubble radius growth. Also, the foam

276

sooner reaches the critical point when the bubbles start to deform from the

277

spherical shape and the walls start to form. However, it can be seen that the

278

final wall thickness is relatively similar in all cases. The walls in foams created

279

with higher amount of blowing agent are slightly thicker at the centre, but they

280

are also slightly thinner near the transition to the strut.

281

The quality of mixing of the initial mixture determines the amount of initial

282

air bubbles, which are whipped into the mixture. The effect of initial number

283

density of bubbles is shown in Fig. 8. Smaller number of initial bubbles will re-

284

sult in the delayed onset of foaming and larger bubbles, because larger amount

16

1

460 440 polyol water

Temperature T (K)

Conversion X

0.8 0.6 0.4 0.2 0 0

model experiment

420 400 380 360 340 320

50

100 Time t (s)

150

300 0

200

50

(a)

Bubble radius R (μm)

Foam density ρfoam (kg m-3)

model experiment

800 600 400 200 50

100 Time t (s)

150

150

200

250 200 150 100 50 0 0

200

50

(c)

100 Time t (s)

(d) 12 Wall thickness h (μm)

Reaction mixture viscosity μrm (Pa s)

200

300

1000

22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 0

150

(b)

1200

0 0

100 Time t (s)

50

100 Time t (s)

150

8 6 4

(e)

Rd

2 0 0

200

t = 47 s t = 57 s t = 85 s t = 123 s t = 161 s t = 199 s

10

25

50 75 100 125 Radial coordiante r (μm)

150

175

(f)

Figure 4: Calculated evolution of conversion, temperature, foam density, bubble radius, viscosity, and wall thickness for default values of parameters. The experimental data were taken from Baser and Khakhar (1994b). The black dotted line shows the final size of computational domain Rd .

17

t=47 s t=62 s t=199 s

150

Dimensions in ( m)

100 50 0 50 100 150 50

0 50 100 Dimensions in ( m)

150

Figure 5: The evolution of strut shape.

Figure 6: SEM image of a strut and a part of the wall. The thickness of the wall is approximately 0.8 µm.

18

6

300

5

Wall thickness h (μm)

Bubble radius R (μm)

350

250 200 150

cw,0 = 298 mol m-3 cw,0 = 555 mol m-3 cw,0 = 780 mol m-3

100 50 0 0

50

100 Time t (s)

150

4 3 2 1 0 0

200

cw,0 = 298 mol m-3 cw,0 = 555 mol m-3 cw,0 = 780 mol m-3

25

50 75 100 125 Radial coordiante r (μm)

(a)

150

175

(b)

Figure 7: Calculated evolution of bubble radius and final wall thickness profile in dependence on initial water concentration.

500 400

6

nb,0 = 0.25 × 1011 m-3 nb,0 = 1.00 × 1011 m-3 nb,0 = 4.00 × 1011 m-3

Wall thickness h (μm)

Bubble radius R (μm)

600

300 200 100 0 0

50

100 Time t (s)

150

5 4 3 2 1 0 0

200

(a)

nb,0 = 0.25 × 1011 m-3 nb,0 = 1.00 × 1011 m-3 nb,0 = 4.00 × 1011 m-3

50

100 150 200 Radial coordiante r (μm)

250

(b)

Figure 8: Calculated evolution of bubble radius and final wall thickness profile in dependence on initial number density of bubbles nb,0 .

285

of blowing agent will be available for each bubble. Larger bubbles will also de-

286

velop with thicker walls between themselves. The discussion of initial dispersion

287

of bubbles and its effect on the final foam morphology is in our work simplified

288

because effects of Ostwald ripening and bubble coalescence are not taken into

289

account.

290

The typical range of values for surface tension in PU foaming process is

291

between 15 and 35 mN m−1 and can be modified by surfactants in practical

292

formulations. Higher surface tension creates larger capillary pressure difference

293

between the wall and the strut. Thus, higher surface tension leads to thinner

294

walls (see Fig. 9). However, it must be noted that higher surface tension also

19

Wall thickness h (μm)

6 γ = 15 mNm-1 γ = 25 mNm-1 γ = 35 mNm-1

5 4 3 2 1 0 0

25

50 75 100 125 Radial coordiante r (μm)

150

175

Figure 9: Calculated final wall thickness profile in dependence on surface tension γ.

6 h0 = 1 μm h0 = 5 μm h0 = 10 μm

5 4

Wall thickness h (μm)

Wall thickness h (μm)

6

3 2 1 0 0

10

20 30 40 Radial coordiante r (μm)

50

5

3 2 1 0 0

60

(a)

h0 = 1 μm h0 = 5 μm h0 = 10 μm

4

25

50 75 100 125 Radial coordiante r (μm)

150

175

(b)

Figure 10: Calculated wall thickness profile for several values of initial inter-bubble distance h0 (a) 30 s after bubble contact and (b) at the end of the simulation.

295

decreases the initial number density of bubbles in the mixing step.

296

Model parameter h0 determines the distance between bubbles at the start of

297

the wall evolution simulation (when foam porosity decreases to porosity of close

298

packed spheres). The influence of this parameter on the predicted wall thickness

299

and strut content is shown in Fig. 10. Initially, the average wall thickness is

300

larger for higher value of parameter h0 . However, the final wall thickness is

301

practically the same for all three chosen cases. The final wall thickness is thus

302

not sensitive to the chosen value of h0 .

20

12

A = 1.6 × 10-7 Pa s A = 8.0 × 10-7 Pa s A = 40.0 × 10-7 Pa s

100

Wall thickness h (μm)

Reaction mixture viscosity μrm (Pa s)

1000

10

1 0

50

100 Time t (s)

150

8 6 4 2 0 0

200

(a)

A = 1.6 × 10-7 Pa s A = 8.0 × 10-7 Pa s A = 40.0 × 10-7 Pa s

10

25

50 75 100 125 Radial coordiante r (μm)

150

175

(b)

Figure 11: Calculated evolution of reaction mixture viscosity and final wall thickness profile in dependence on viscosity parameter A (cf. eq. (25)).

303

Thin walls and high strut content in PU foams are caused by the low viscosity

304

of the reaction mixture (Mills, 2007). The calculated evolution of viscosity and

305

final wall thickness profile for three cases of viscosity parameter A (cf. eq. (25))

306

are shown in Fig. 11. As can been seen the increased viscosity truly leads to

307

thicker walls because the drainage is hindered by the increased resistance in the

308

film.

309

It should be noted that the viscosity of reaction mixture can be different at

310

various places in the foam, because temperature hot spots can be created in PU

311

foams even in relatively simple mould geometries (Karimi et al., 2016). These

312

will decrease the reaction mixture viscosity and could lead to locally thinner

313

walls and broader bubble size distribution due to increased coalescence, which

314

would result in the spatial variation of thermal and mechanical properties.

315

5. Conclusions

316

In this work, mathematical model for the simulation of PU foam morphology

317

development is presented. It combines reaction kinetics modelled by the end

318

group model with the simulation of bubble growth using the bubble–shell model

319

and with the simulation of wall evolution using thin film hydrodynamics. The

320

coupled model offers an opportunity to study the foaming process in great detail

21

321

as it provides the evolution of temperature, concentration of chemical species,

322

foam density, bubble size and wall thickness profile at the same time.

323

Morphology and thus quality of the foam is influenced by many parameters.

324

The amount of blowing agent and number density of initial bubbles are critical

325

for the final foam density and cell size, whereas the viscosity of reaction mixture

326

and the growth of bubbles are responsible for the evolution of walls. Wall

327

thickness profiles are essentially important for application properties of closed-

328

cell PU foams as they determine mechanical and heat insulation properties and

329

also the rate of blowing agent permeation out of the foam. Both wall drainage

330

and wall stretching by bubble growth are contributing to the final wall thickness

331

profiles.

332

The model captures the most important phenomena responsible for the foam

333

production and the simulation results show acceptable level of agreement with

334

the real process. In the future, the proposed methodology can be further im-

335

proved by including the effect of Marangoni convection and more realistic de-

336

scription of rheological behaviour close to the gel point. Further work will

337

focus on implementation of the CFD macro-scale model for the simulation of

338

industrial-scale mould filling applications.

339

Acknowledgements

340

The work was supported by Czech Science Foundation (Project No. GA14-

341

18938S). The financial support from specific university research (MSMT No.

342

20-SVV/2017) is gratefully acknowledged. We also thank Andra Nistor for

343

providing SEM images.

344

References

345

Ahern, A., Verbist, G., Weaire, D., Phelan, R., Fleurent, H., aug 2005. The

346

conductivity of foams: a generalisation of the electrical to the thermal case. 22

347

Colloids and Surfaces A: Physicochemical and Engineering Aspects 263 (1-3),

348

275–279.

349

URL

350

S0927775705000920

http://linkinghub.elsevier.com/retrieve/pii/

351

Amon, M., Denson, C. D., sep 1984. A study of the dynamics of foam growth:

352

Analysis of the growth of closely spaced spherical bubbles. Polymer Engineer-

353

ing and Science 24 (13), 1026–1034.

354

URL http://doi.wiley.com/10.1002/pen.760241306

355

Arefmanesh, A., Advani, S. G., 1991. Diffusion-induced growth of a gas bubble

356

in a viscoelastic fluid. Rheologica Acta 30 (3), 274–283.

357

URL http://link.springer.com/10.1007/BF00366641

358

Baser, S. A., Khakhar, D. V., apr 1994a. Modeling of the dynamics of R-11

359

blown polyurethane foam formation. Polymer Engineering and Science 34 (8),

360

632–641.

361

URL http://doi.wiley.com/10.1002/pen.760340804

362

Baser, S. A., Khakhar, D. V., apr 1994b. Modeling of the Dynamics of Water and

363

R-11 blown polyurethane foam formation. Polymer Engineering and Science

364

34 (8), 642–649.

365

URL http://doi.wiley.com/10.1002/pen.760340805

366

Bergeron, V., may 1999. Forces and structure in thin liquid soap films. Journal

367

of Physics: Condensed Matter 11 (19), R215–R238.

368

URL http://stacks.iop.org/0953-8984/11/i=19/a=201?key=crossref.

369

b643cc4f833da948943d3f60099bee3d

370

371

Bhakta, A., Ruckenstein, E., 1997. Decay of standing foams: drainage, coalescence and collapse. Advances in Colloid and Interface Science 70, 1–124.

23

372

URL

373

S0001868697000316

374

http://linkinghub.elsevier.com/retrieve/pii/

Bhakta, A. R., Khilar, K. C., aug 1991. A network simulation for drainage of

375

static foam columns. Langmuir 7 (8), 1827–1832.

376

URL http://pubs.acs.org/doi/abs/10.1021/la00056a041

377

Bikard, J., Bruchon, J., Coupez, T., Silva, L., nov 2007. Numerical simulation

378

of 3D polyurethane expansion during manufacturing process. Colloids and

379

Surfaces A: Physicochemical and Engineering Aspects 309 (1-3), 49–63.

380

URL

381

S0927775707003044

382

http://linkinghub.elsevier.com/retrieve/pii/

Castro, J. M., Macosko, C. W., mar 1982. Studies of mold filling and curing in

383

the reaction injection molding process. AIChE Journal 28 (2), 250–260.

384

URL http://doi.wiley.com/10.1002/aic.690280213

385

Chen, Y., Das, R., Battley, M., jan 2015. Effects of cell size and cell wall

386

thickness variations on the stiffness of closed-cell foams. International Journal

387

of Solids and Structures 52, 150–164.

388

URL

389

S0020768314003692

http://linkinghub.elsevier.com/retrieve/pii/

390

Feng, J. J., Bertelo, C. A., 2004. Prediction of bubble growth and size distri-

391

bution in polymer foaming based on a new heterogeneous nucleation model.

392

Journal of Rheology 48 (2), 439.

393

URL http://link.aip.org/link/JORHD2/v48/i2/p439/s1{&}Agg=doi

394

Ferkl, P., Karimi, M., Marchisio, D. L., Kosek, J., jul 2016. Multi-scale mod-

395

elling of expanding polyurethane foams: Coupling macro- and bubble-scales.

396

Chemical Engineering Science 148, 55–64.

24

397

URL

398

S0009250916301543

http://linkinghub.elsevier.com/retrieve/pii/

399

Ferkl, P., Toulec, M., Laurini, E., Pricl, S., Fermeglia, M., Auffarth, S., Eling,

400

B., Settels, V., Kosek, J., nov 2017. Multi-scale modelling of heat transfer in

401

polyurethane foams. Chemical Engineering Science 172, 323–334.

402

URL

403

S0009250917304244

http://linkinghub.elsevier.com/retrieve/pii/

404

Geier, S., Winkler, C., Piesche, M., sep 2009. Numerical Simulation of Mold

405

Filling Processes with Polyurethane Foams. Chemical Engineering & Tech-

406

nology 32 (9), 1438–1447.

407

URL http://doi.wiley.com/10.1002/ceat.200900202

408

Harikrishnan, G., Khakhar, D. V., 2009. Modeling the dynamics of reactive

409

foaming and film thinning in polyurethane foams. AIChE Journal 56 (2),

410

522–530.

411

URL http://doi.wiley.com/10.1002/aic.12002

412

413

414

Hindmarsh, A. C., 1980. ODEPACK library, http://www.netlib.org/odepack/. URL http://www.netlib.org/odepack/ Joye, J. L., Miller, C. A., Hirasaki, G. J., 1992. Dimple formation and behavior

415

during axisymmetrical foam film drainage. Langmuir 8 (12), 3083–3092.

416

URL

417

0-0026976357{&}partnerID=tZOtx3y1

http://www.scopus.com/inward/record.url?eid=2-s2.

418

Karimi, M., Droghetti, H., Marchisio, D. L., feb 2016. Multiscale Modeling

419

of Expanding Polyurethane Foams via Computational Fluid Dynamics and

420

Population Balance Equation. Macromolecular Symposia 360 (1), 108–122.

421

URL http://doi.wiley.com/10.1002/masy.201500108

25

422

Karimi, M., Marchisio, D. L., jul 2015. A Baseline Model for the Simulation of

423

Polyurethane Foams via the Population Balance Equation. Macromolecular

424

Theory and Simulations 24 (4), 291–300.

425

URL http://doi.wiley.com/10.1002/mats.201500014

426

Kim, C., Youn, J. R., feb 2000. Environmentally Friendly Processing of

427

Polyurethane Foam for Thermal Insulation. Polymer-Plastics Technology and

428

Engineering 39 (1), 163–185.

429

URL http://www.tandfonline.com/doi/abs/10.1081/PPT-100100022

430

K¨ orner, C., Thies, M., Hofmann, T., Th¨ urey, N., R¨ ude, U., oct 2005. Lattice

431

Boltzmann Model for Free Surface Flow for Modeling Foaming. Journal of

432

Statistical Physics 121 (1-2), 179–196.

433

URL http://link.springer.com/10.1007/s10955-005-8879-8

434

Lee, S. S., Kim, S. C., aug 1991. Analysis of unsaturated polyester-polyurethane

435

interpenetrating polymer networks. I: Reaction kinetics and viscosity model.

436

Polymer Engineering and Science 31 (16), 1182–1188.

437

URL http://doi.wiley.com/10.1002/pen.760311605

438

Manev, E. D., Nguyen, A. V., jun 2005. Critical thickness of microscopic thin

439

liquid films. Advances in colloid and interface science 114-115, 133–46.

440

URL http://www.ncbi.nlm.nih.gov/pubmed/15936287

441

442

Mills, N. J., 2007. Polymer foams handbook: engineering and biomechanics applications and design guide. Butterworth-Heinemann.

443

Mitrias, C., Jaensson, N. O., Hulsen, M. A., Anderson, P. D., jun 2017. Direct

444

numerical simulation of a bubble suspension in small amplitude oscillatory

445

shear flow. Rheologica Acta 56 (6), 555–565.

446

URL http://link.springer.com/10.1007/s00397-017-1009-0

26

447

Nistor, A., Toulec, M., Zubov, A., Kosek, J., feb 2016. Tomographic Recon-

448

struction and Morphological Analysis of Rigid Polyurethane Foams. Macro-

449

molecular Symposia 360 (1), 87–95.

450

URL http://doi.wiley.com/10.1002/masy.201500113

451

Pardo-Alonso, S., Sol´ orzano, E., Brabant, L., Vanderniepen, P., Dierick,

452

M., Van Hoorebeke, L., Rodr´ıguez-P´erez, M., may 2013. 3D Analysis of

453

the progressive modification of the cellular architecture in polyurethane

454

nanocomposite foams via X-ray microtomography. European Polymer Jour-

455

nal 49 (5), 999–1006.

456

URL

457

S0014305713000232

http://linkinghub.elsevier.com/retrieve/pii/

458

Sahu, S., Gokhale, A., Mehra, A., jan 2014. Modeling nucleation and growth of

459

bubbles during foaming of molten aluminum with high initial gas supersatu-

460

ration. Journal of Materials Processing Technology 214 (1), 1–12.

461

URL

462

S0924013613002355

http://linkinghub.elsevier.com/retrieve/pii/

463

Samkhaniani, N., Gharehbaghi, A., Ahmadi, Z., jun 2013. Numerical simulation

464

of reaction injection molding with polyurethane foam. Journal of Cellular

465

Plastics 49 (5), 405–421.

466

URL http://cel.sagepub.com/cgi/doi/10.1177/0021955X13485594

467

Schwartz, L., Roy, R., aug 2003. A mathematical model for an expanding foam.

468

Journal of Colloid and Interface Science 264 (1), 237–249.

469

URL

470

S0021979703004259

http://linkinghub.elsevier.com/retrieve/pii/

471

Seo, D., Youn, J. R., aug 2005. Numerical analysis on reaction injection

472

molding of polyurethane foam by using a finite volume method. Polymer 27

473

46 (17), 6482–6493.

474

URL

475

S0032386105006397

http://linkinghub.elsevier.com/retrieve/pii/

476

Tuladhar, T., Mackley, M., dec 2004. Experimental observations and modelling

477

relating to foaming and bubble growth from pentane loaded polystyrene

478

melts. Chemical Engineering Science 59 (24), 5997–6014.

479

URL

480

S0009250904004920

http://linkinghub.elsevier.com/retrieve/pii/

481

Yasunaga, K., Neff, R. A., Zhang, X. D., Macosko, C. W., 1996. Study of cell

482

opening in flexible polyurethane foam. Journal of Cellular Plastics 32 (5),

483

427–447.

484

URL

485

0-0030245318{&}partnerID=tZOtx3y1

http://www.scopus.com/inward/record.url?eid=2-s2.

486

Zhang, X. D., Davis, H. T., Macosko, C. W., 1999. New cell opening mechanism

487

in flexible polyurethane foam. Journal of Cellular Plastics 35 (5), 458–476.

488

URL

489

0-0032643099{&}partnerID=tZOtx3y1

http://www.scopus.com/inward/record.url?eid=2-s2.

28

Highlights

• • • •

Modelling tool for the prediction of polyurethane foam morphology. Coupling of reaction kinetics, bubble growth and wall evolution models. Demonstrated effect of surface tension and viscosity on final morphology. Wall thickness profiles and strut shapes reproduced in simulated foams.

1