# Numerical modelling of thermoplastic CO2 assisted extrusion foaming: Exploring the process-microstructure relationship

Mustafa Altinisik, Margot Chauvet, Martial Sauceau, Raffaele D'Elia, Romain Sescousse, Fabien Baillon, Guilhem Michon

► To cite this version:

Mustafa Altinisik, Margot Chauvet, Martial Sauceau, Raffaele D'Elia, Romain Sescousse, et al.. Numerical modelling of thermoplastic CO2 assisted extrusion foaming: Exploring the process-microstructure relationship. ESAFORM 2025 - 28th International ESAFORM Conference on Material Forming, May 2025, Paestum, Italy. pp.2391-2398, 10.21741/9781644903599-258 . hal-05272232

HAL Id: hal-05272232

https://imt-mines-albi.hal.science/hal-05272232v1

Submitted on 22 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical modelling of thermoplastic CO₂ assisted extrusion foaming: Exploring the process-microstructure relationship

Mustafa ALTINISIK¹,ª, Margot CHAUVET³,ᵇ, Martial SAUCEAU²,ᶜ, Raffaele D'ELIA¹,ᵈ, Romain SESCOUSSE²,ᵉ, Fabien BAILLON²,ᵉ
Guilhem MICHON¹,⁹

¹Institut Clément Ader (ICA), UMR-CNRS 5312, "INSA, UPS, Mines Albi, ISAE" 31077, Toulouse, France

²Centre RAPSODEE IMT Mines Albi, CNRS, Université de Toulouse, 81013, Albi, France

³IPREM, Technopôle Helioparc, 2 avenue P. Angot, 64053 Pau Cedex 9, France

ªmustafa.altinisik@isae-supaero.fr, ᵇmargot.chauvet@univ-pau.fr,

ᶜmartial.sauceau@mines-albi.fr, ᵈraffaele.delia@mines-albi.fr,

ᵉromain.sescousse@mines-albi.fr, Ṡfabien.baillon@mines-albi.fr,

⁹guilhem.michon@isae-supaero.fr

**Keywords:** Supercritical Fluid, Thermoplastic Polymer, Simulation, Acoustic, Foaming

**Abstract.** This study presents a comprehensive numerical modeling approach for thermoplastic foams processing, produced via supercritical CO₂-assisted extrusion foaming. The model focuses on the nucleation and growth of bubbles within the polymer matrix, employing the influence volume approach to simulate these dynamics. Key phenomena such as gas diffusion, nucleation rate, bubble growth, and the evolution of influence and non-influence volumes are captured to predict foam microstructure and cellular density. Results highlight the model's capability to predict equilibrium bubble sizes ranging from micrometers to hundreds of micrometers, offering insights into foam behavior and microstructure control. The final size distribution will be used to create on Paraview® the microstructure that will be compared to experimental observations in the future.

## Introduction

Sandwich structures are widely used in lightweight applications due to their enhanced strength-to-weight ratio. Nomex® honeycombs are commonly employed as core materials because of their excellent mechanical properties, ease of processing, and compliance with fire, smoke, and toxicity standards required by transportation regulations. However, their limitations in sound insulation, damping capabilities, and recyclability remain significant concerns. A promising alternative involves the use of multifunctional thermoplastic foams developed through supercritical CO₂-assisted extrusion foaming technology. This environmentally friendly, non-flammable solvent/solute serves as a physical blowing agent during the extrusion process, enabling precise control over the foaming process and, consequently, the final foam microstructure and functionality. The porosity, size, and density of the bubbles, as well as the type of porosity—whether closed or open cells—determine the foam's multifunctionality, as illustrated in Table 1. To assess the foam's functionality, it is essential to model not only the nucleation, growth, and coalescence processes but also the microstructure itself. The nucleation and growth processes have been previously modeled by Shafi et al. [1], whose work predicts bubble size distribution in LDPE-N₂ (Low-Density Polyethylene) systems using a batch approach. They assume that the molten polymer contains a finite amount of gas and, by introducing the "influence volume" approach, successfully account for nucleation and growth phenomena simultaneously under isobaric conditions. Similarly, Taki et al. [2] employed this approach while incorporating a continuously decreasing pressure profile, which allows for the consideration of property evolutions, such as

<figure><img src="image_2.png" /></figure>
viscosity, diffusion coefficient, nucleation, and growth rates. Shimoda et al. [3] further refined Shafi’s model by including the continuous evolution of properties and accounting for polymer flow during the extrusion process.

<table><caption>Table 1 – Factors determining foams functionality [4]</caption><thead><tr><th>Functional Applications</th><th>Cell density</th><th>Cell size</th><th>Expansion ratio</th><th>Porosity morphology</th></tr></thead><tbody><tr><td>Mechanical strength</td><td>--</td><td>++</td><td>--</td><td>Closed porosity</td></tr><tr><td>Weight reduction</td><td>+</td><td>+</td><td>++</td><td>/</td></tr><tr><td>Thermal insulation</td><td>--</td><td>+</td><td>+</td><td>Closed porosity</td></tr><tr><td>Sound insulation</td><td>-</td><td>+</td><td>+</td><td>Open porosity</td></tr></tbody></table>

The objective for this modeling is to implement the model by Shafi et al. [5] to predict cellular density. This will help to identify the kinetics of nucleation and growth, providing a better understanding of the phenomenological changes during foaming. Initially, the flow within the die was considered, as was done by Shimoda et al. but due to the complexity of the model developed by Shafi et al. [6] only the pressure changes in the polymer to represent extrusion foaming was considered. This approach is therefore based on the method developed by Taki et al. for pressure change modeling.

## Model description and formulation

In this approach, an influence volume, shown in Fig.1, $V_s(t)$ is considered around each bubble. The influence volume approach involves considering that there is a certain volume around a bubble, within which the nucleation of another bubble cannot occur. The influence volume $V_s$ depends on time and can be calculated using Eq. 1:

$$V_s(\tau) = \frac{4}{3}(S^3(\tau) - R^3(\tau)) \qquad (1)$$

where $\tau$ is the time since the nucleation of the bubble in [s], R the radius of the bubble in [m] and S, the radial position where the gas concentration equals the nucleation concentration threshold, which is the limit where bubble growth cannot occur due to gas depletion. As the bubble grows, the concentration gradient propagates radially within the polymer, leading to an increase in the influence volume. Knowledge of the operational pressure profile allows for the calculation of the bubble's radius, the pressure inside the bubble, and the concentration profile in the polymer surrounding it, based on the governing growths equations.

<figure><img src="image_3.png" /><figcaption>Figure 1 – Influence volume approach [1]</figcaption></figure>
Based on the influence volumes of each bubble, it is possible to define a non-influence, $V_L(t)$, within which new bubbles can form. The non-influence volume is calculated by using the following Eq. 2:

$$V_L(t) = V_L(0) - \int_0^t V_L(t')J(t')V_s(t-t')dt' \quad (2)$$

where $J(t)$ is the nucleation rate in [bubble.m⁻³.s⁻¹] and $V_L(0)$ is the initial non-influence volume in [m³]. Nucleation is the formation of a new phase from an initial phase, such as the formation of bubbles in a liquid. Classical Nucleation Theory (CNT) predicts bubble nucleation [7]. According to Gibbs et al., bubble nucleation occurs either within a continuous liquid phase (i.e., homogeneous nucleation) or at a liquid/liquid or liquid/solid interface (i.e., heterogeneous nucleation) [7]. At this stage, only homogeneous nucleation is considered, and nucleation rate, is determined using the following Eq. 3:

$$J(t) = C_a(t)N_a \sqrt{\frac{2\gamma}{\pi M}} \exp\left(-\frac{16\pi\gamma^3}{3k_B T (P_{\text{bubble}}(t) - P_{\text{polymer}}(t))^2}\right) \quad (3)$$

where $C_a$ is the average concentration of dissolved gas in [mol.m⁻³], $N_a$ is Avogadro's number in [mol⁻¹], $\gamma$ is the surface tension in [N.m⁻¹], $M$ is the molar mass of the blowing agent in [g.mol⁻¹], $k_b$ is the Boltzmann constant in [J.K⁻¹], $T$ is the temperature of the polymer in [K], $P_{\text{bubble}}$ is the pressure inside the bubble in [Pa], and $P_{\text{polymer}}$ is the pressure in the continuous phase (the polymer) in [Pa].

Using the influence volume and the nucleation rate, it is possible to determine the number of nucleated bubbles $N_b$ in [bubble.m⁻³], which is calculated with Eq. 4:

$$\frac{dN_B}{dt} = J(t)V_L(t) \quad (4)$$

The average concentration of dissolved gas is calculated by the following Eq. 5:

$$C_a V_L = c_0 V_{L,0} - \int_0^t \frac{4\pi}{3} J(t') V_L(t') \frac{P_{\text{bubble}}(t-t')R^3(t-t')}{RT} dt' \\ - \int_0^t J(t') V_L(t') \int_{R(t-t')}^{S(t-t')} 4\pi r^2 c(r, t-t') dr dt' \quad (5)$$

where $R$ is the bubble radius, $r$ is the radial coordinate as shown in Fig. 1, $\Re$ is the universal gas constant, $c_0$ the initial concentration of dissolved gas in the polymer phase and $c$ the concentration of dissolved gas in the polymer phase. The growth of bubbles depends on the diffusion coefficient of the blowing agent, which is concentration-dependent, the gradual cooling of the foam during expansion, the influence of the blowing agent on the polymer's viscosity, and the release of the blowing agent from the foam surface as Ramesh et al. established [8]. The accumulation of gas molecules adsorbed on the bubble surface is negligible, the mass conservation law implies that the change in mass within the gas bubble is balanced by the mass of gas diffusing into the bubble through its surface. Furthermore, considering that the gas inside the bubble follows the ideal gas law, the pressure within the bubble can therefore be related to the concentration gradient at the bubble surface as shown in Eq 6. It is assumed that gas diffusion coefficient in the polymer/gas solution remains constant during growth, the diffusion equation can be written as Eq. 7. Finally, the evolution of the bubble radius is therefore determined by a balance between the differential pressure, which promotes growth, and the interfacial tension, which opposes it. This equation is commonly used to model bubble growth in polymer systems under the effect of gas diffusion and internal pressure as described in Eq. 8.
$$ \frac{d}{dt} \left( \frac{4\pi P_{\text{Bubble}} R^3}{3 \mathfrak{R} T} \right) = 4\pi R^2 D \left. \frac{\partial c}{\partial r} \right|_{r=R} \qquad (6) $$

$$ \frac{\partial c}{\partial t} + \frac{\dot{R} R^2}{r^2} \frac{\partial c}{\partial r} = \frac{D}{r^2} \frac{\partial}{\partial t} \left( r^2 \frac{\partial c}{\partial r} \right) \qquad (7) $$

$$ \frac{dR}{dt} = \frac{(P_{\text{bubble}} - P_{\text{polymer}})R}{4\eta} - \frac{\gamma}{2\eta} \qquad (8) $$

Here, D is the diffusion coefficient in [m².s⁻¹], η is the viscosity of the polymer in [Pa.s] and K_h is the Henry's law constant (also called the solubility coefficient) in [mol.N⁻¹.m⁻¹]. The initial and boundary conditions for the concentration are as follows in Eq. 9, Eq. 10, and Eq. 11:

$$ c(r, 0) = c_i(r) \qquad (9) $$

$$ c(R, t) = c_R(t) = K_H P_{\text{bubble}}(t) \qquad (10) $$

$$ c(\infty, t) = c_0 = K_H P_{\text{bubble}}^0 \qquad (11) $$

Here, $P_{\text{bubble}}^0$ is the initial gaz pressure in the bubble. The numerical simulation explication is described in the following part.

## Numerical simulation

Fig. 2 describes the flowchart of the functional algorithm. The model is initialized by assuming that the non-influence volume $V_L$ is equal to the total volume $V_T$, as there are no pre-existing bubbles. The concentration of dissolved gas in the polymer is known and equals $c_0$, which implies that the initial bubble pressure, $P_{\text{bubble}}^0$ is also known. After initialization, the first step is to calculate the initial nucleation rate $J_0$. From this, the threshold concentration $c_s$ can be determined, which will be reached when the nucleation and growth of various bubble populations consume most of the dissolved gas (such that $J_s \le 0.01 J_0$). Following initialization, the program enters its main loop. This loop discretizes time into intervals of $\Delta t$. For each time step, it is checked whether the non-influence volume is sufficient to continue ($V_L \ge 0.01 V_T$). If this condition is met and if we have a saturated environment ($J > J_s$), a new bubble population appears, all born with a radius R and a concentration $c_R$. Nucleation is immediately followed by bubble growth, leading to changes in the radius R(t) and concentration $c_R(t)$. At the end of each time step, it is necessary to calculate the influence volume $V_s$ around each bubble (the larger the bubbles are at their formation, the greater the influence volume). Knowing the number of bubbles in each population, the new non-influence volume $V_L$ can be calculated with the following Eq. 12. as well as the average gas concentration $C_a$ in the continuous phase. The new nucleation rate J is then calculated for the next time step.

If the non-influence volume $V_L$ becomes negligible ($V_L < 0.01 V_T$), the main loop is exited. This is followed by a correction step for the total number of bubbles, to eliminate any populations that result in a negative non-influence volume (excess populations) or a negligible influence volume (in other words, bubble populations that may have formed but are ultimately too small to be viable). One of the assumptions of this modeling approach is that bubbles do not interact with each other and therefore do not coalesce. As a result, at the final time $t_f$, when the non-influence volume has
become too small and nucleation has stopped, the number of bubbles also becomes constant and known. However, the bubbles continue to grow by consuming the dissolved gas within the influence volume $V_s(t)$ until they reach equilibrium. A correction of the final equilibrium radius is performed using a global material balance across the influence volumes of each bubble population. In the end, with the equilibrium radius $R$ of the bubbles in each population and the number of bubbles in each population known, it becomes possible to calculate key characteristics of the formed foam, such as the bubble size distribution and the apparent global porosity of the polymer.

<figure><img src="image_4.png" /><figcaption>Figure 2 – Flowchart of the functional algorithm</figcaption></figure>

## Results

To validate the model, semi-crystalline polyethylene parameters was used to compare with results from Mao et al [9] by using the following Table 2 parameters taken in Shafi's et al literature. To achieve the same results, it was necessary to multiply the nucleation rate by a factor of $10^7$. It should be noted that the model by Shafi et al. refers to a nucleation prefactor whose value is never explicitly stated. The hypothesis put forward regarding this nucleation factor is that it represents the heterogeneous nucleation not accounted for in the model. Also, all figures were plot with simulation done with $\Delta t$ : 0.5 s.

<table><caption>Table 2 – Simulation data for a LDPE-N₂ system [1]</caption><thead><tr><td>Parameters</td><td>Value</td></tr></thead><tbody><tr><td>Diffusion coefficient (const)</td><td>4.26 × 10<sup>-9</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Solubility coefficient (const)</td><td>3.61 × 10<sup>-5</sup> mol N<sup>-1</sup> m<sup>-1</sup></td></tr><tr><td>Viscosity (const)</td><td>4.9 × 10<sup>4</sup> Pa.s</td></tr><tr><td>Surface tension (const)</td><td>11.5 × 10<sup>-3</sup> N m<sup>-1</sup></td></tr><tr><td>Initial bubble pressure</td><td>1.11 × 10<sup>7</sup> Pa</td></tr><tr><td>Polymer pressure (const)</td><td>1.01 × 10<sup>5</sup> Pa</td></tr><tr><td>Temperature (const)</td><td>423.8 K</td></tr></tbody></table>

Fig. 3 illustrates the nucleation rate and the non-influence volume as a function of time, as well as the number of bubbles as a function of size. Fig. 3A shows that the nucleation rate remains constant for over 0.3 seconds before dropping sharply as the gas concentration in the dissolved medium decreases. The non-influence volume initially, shown in Fig. 3B, remains constant but then decreases significantly starting at 0.25 seconds. Fig. 3C shows that the pressure inside nucleated bubble decrease rapidly from $90 \times 10^5$ to $8 \times 10^5$ Pa. This is due to gas depletion as shown in Fig. 3D where gas concentration in bubbles drop.
<figure><img src="image_5.png" /><figcaption>Figure 3 – (A) Nucleation rate during simulation; (B) Non-influence volume during simulation; (C) Bubble pressure during simulation; (D) Gas concentration during simulation</figcaption></figure>

As for Fig. 4, populations that nucleate near last $\Delta t$ have less and smaller bubbles in it due to the decrease of non-influence volume.

<figure><img src="image_6.png" /><figcaption>Figure 4 – Bubble radius during simulation</figcaption></figure>

Fig. 5 shows the equilibrium bubble size as a function of bubble age. The bubbles on the far right are the oldest and therefore the largest, while those on the left are the most recent and the smallest. It is observed that the bubble size ranges from 0.5 µm to 238 µm. By knowing the number of bubbles and their final equilibrium radius, it is possible to calculate an average equilibrium radius, which is 119.25 µm.
<figure><img src="image_7.png" /><figcaption>Figure 5 – Final bubbles sizes at equilibrium</figcaption></figure>

Table 3 and Fig. 6 illustrates how Δt impacts the simulation. A decrease in Δt significantly increases the number of bubble populations. However, this increase in the number of populations does not correspond to a rise in the total number of bubbles. The reduction in Δt affects the characteristic sizes of equilibrium bubbles, leading to an increase in the average radius size. The impact of Δt diminishes at lower values, but as Δt increases, the results become more difficult to predict.

<table><caption>Table 3 – Δt influence</caption><thead><tr><th>Δt (s)</th><th>Nucleation start at (s)</th><th>Nucleation end at (s)</th><th>Populations</th><th>Bubbles</th><th>R<sub>min</sub> (μm)</th><th>R<sub>max</sub> (μm)</th><th>R<sub>average</sub> (μm)</th></tr></thead><tbody><tr><td>0,005</td><td>0,110</td><td>0,510</td><td>81</td><td>1,2134 x 10<sup>13</sup></td><td>0,5</td><td>215</td><td>107,75</td></tr><tr><td>0,01</td><td>0,11</td><td>0,54</td><td>44</td><td>1,2342 x 10<sup>13</sup></td><td>0,5</td><td>238</td><td>119,25</td></tr><tr><td>0,05</td><td>0,150</td><td>0,450</td><td>10</td><td>1,1869 x 10<sup>13</sup></td><td>5</td><td>285</td><td>145</td></tr><tr><td>0,1</td><td>0,200</td><td>0,4</td><td>3</td><td>1,1561 x 10<sup>13</sup></td><td>18</td><td>135</td><td>76,5</td></tr></tbody></table>

The convergence criteria for the different steps have been satisfied, and the independence of the results with Δt over a certain range appears to be reassuring. Overall, the results remain consistent for Δt<0.015. However, for Δt≥0.015 s, the values are overestimated, corresponding to apparent cumulative values due to discretization errors. In conclusion, using a very small Δt (<0.005 s) is unnecessary, as it has little effect on bubble size or count while increasing computation time.

## Conclusion

This study presents a comprehensive numerical modeling approach to predict the behavior of thermoplastic foams during supercritical CO₂-assisted extrusion foaming. The key focus is on understanding the relationship between the foaming process and the resulting microstructure, with an emphasis on nucleation and bubble growth dynamics. By implementing the influence volume approach, the model successfully captures the essential phenomena of homogeneous bubble nucleation and growth. This approach considers both the influence and non-influence volumes, enabling a detailed prediction of bubble population dynamics and equilibrium characteristics. The nucleation rate remains constant initially before declining as gas concentration in the polymer decreases. This aligns with the consumption of dissolved gas during bubble formation. Bubble sizes at equilibrium range from micrometers to hundreds of micrometers, with an average size influenced by simulation parameters like the time step Δt. A reduction in Δt increases the number of populations but has minimal effect on the total number of bubbles or average bubble size.
<figure><img src="image_8.png" /><figcaption>Figure 6 – Number of bubbles for different $\Delta t$</figcaption></figure>

Conversely, an increase in $\Delta t$ results in fewer populations but larger bubble sizes. The model was validated against literature data, requiring a correction factor for the nucleation rate due to an unspecified prefactor in the referenced model. The model enables the calculation of critical foam properties, including bubble size distribution and porosity. These outputs are essential for tailoring the microstructure to specific applications. Future work will involve correlating experimental observations with the model, incorporating heterogeneous nucleation, coalescence phenomena to improve the model's applicability to more complex systems and building of the realistic microstructure by depicting a 1 mm³ representative volume element.

## References

[1] M.A. Shafi, J.G. Lee, and R.W. Flumerfelt, Prediction of cellular structure in free expansion polymer foam processing. Polym Eng Sci. 36 (1996) 1950-1959. https://doi.org/10.1002/pen.10591

[2] K. Taki, H. Hayashizaki, K. Fukada, ISIJ International, 54 (2014), 2493-2502
https://doi.org/10.2355/isijinternational.54.2493

[3] M. Shimoda, I. Tsujimura, M. Tanigaki et M. Ohshima, Polymeric Foaming Simulation for Extrusion Processes, J. Cellular Plastics, 37.6 (2001) 517-536, https://doi.org/10.1106/W4C0-CAG6-H3FM-LTJP

[4] G. Wu, P. Xie, H. Yang, A review of thermoplastic polymer foams for functional applications. J Mater Sci, 56 (2021) 11579-11604. https://doi-org/10.1007/s10853-021-06034-6

[5] M. A. Shafi and R. W. Flumerfelt, Chemical Engineering Science, 52.4 (1997) 627-633
Initial bubble growth in polymer foam processes. https://doi.org/10.1016/S0009-2509(96)00434-4

[6] M. A. Shafi, K. Joshi and R. W. Flumerfelt, Bubble size distributions in freely expanded polymer foams, Chemical Engineering Science, 52.4 (1997) 635-644,
https://doi.org/10.1016/S0009-2509(96)00433-2

[7] J. W. Gibbs. The Scientific Papers of J. Willard Gibbs, Volume 1. Dover Publications Inc., (1961) 43

[8] N.S. Ramesh. « Chapter 7, Foam Growth in Polymers ». Foam Extrusion - Principles and Practice, Edited by S.-T. Lee and Chul B. Park. (2014) 213-237

[9] D. Mao, J. R. Edwards and A. Harvey. Prediction of foam growth and its nucleation in free and limited expansion. Chem. Engineering Science, 61.6 (2006), 1836-1845
https://doi.org/10.1016/j.ces.2005.07.026