JEC COMPOSITES MAGAZINE - Issue #114 - July 2017 - 56

TECHNOLOGY software
n, and the Nakamura kinetics function, K. These two new material properties are added in the Physics Builder. Note that the
user will eventually have to input the thermal dependency of K
in the model (see example in the Application section). In terms
of parameter input for the Nakamura model, only two scalars,
αmin and αmax, are required.
The exotherm, Q, (Equation 2) can be computed automatically
in the Physics Builder by adding a domain feature node that includes the domain variable exotherm. In addition to this node,
two new material properties (ρm and φm) are needed.
Coefficient form PDE
The right-hand side of Equation (1) is implemented by adding a
new domain variable, RHS. Then, the ordinary differential equation (1) is solved by adding a coefficient form equation in the
Physics Builder.
The newly developed physics node is then added in the local
physics library using the Physics Manager. It is now readily available to the COMSOL user to be added in multiphysical models.

Application: simulating the cooling
of a thermoplastic composite plate
In this section, the quiescent cooling of a thermoplastic composite square is modelled in COMSOL. Heat transfer is handled using a heat transfer node whereas crystallization physics is added
using the newly developed node. Full coupling shown in Figure
1 is taken into account.
Modelling
A 2x2cm square is sketched in a 2D model in COMSOL.
Standard heat transfer in a solid node is added. The thermal
material properties are adapted from the literature and used to
model an orthotropic carbon/PEEK composite [4, 9, 10]. The
material's thermal properties depend on the temperature and
the degree of crystallization. The properties are usually characterized versus temperature for amorphous (α=0) and fully crystalline (α=1) states. A mixture law is then used to determine the
thermo- and crystallo-dependent properties.
An initial temperature of 300°C is imposed in the whole domain.
Symmetry is imposed on the left and lower boundaries. A convective heat flux at an external temperature of 230°C with an
exchange coefficient h=20Wm-2K-1 is considered on the upper
and right boundaries. The heat source arising from the crystallization exotherm is also added automatically by calling the exotherm variable.
The crystallization node is added in the model and allows the
resolution of Nakamura crystallization kinetics at each point.
The missing material properties automatically appear in the
material node. Values from the literature [9] were implemented.
Note that the temperature dependency is accounted for.
An initial value of the degree of crystallization is set to 0 in the
whole domain and used to model an initially fully molten material. Note that even if imposing α=0, crystallization will be
initiated artificially thanks to the non-zero value of αmin. An

jec composites magazine / N°114 July 2017

Fig. 4: Predicting the cooling and solidification of a composite plate. The
surface plot represents the degree of crystallization, the lines are isothermals and the arrows the heat flux. The plot is at time t=90s

exotherm computation domain feature is also added for a fully
coupled resolution, as depicted in Figure 1.
A default triangular mesh using a "fine" size is used. The time-dependent study is performed between 0 and 150 seconds. Solutions
are stored every 10 seconds. The solving methods are the default
methods proposed by COMSOL, except that a maximum time
step size of 1 second is imposed in the time-dependent solver.
Indeed, crystallization is an event that will occur while cooling
down the material. Large time steps may result in missing this
event. Correct time steps, featuring automatic time stepping with
a maximum time step size, is thus required [4].
Results
The temperature and degree of crystallization fields are obtained
versus time. This information is shown in Figure 9 at time t = 90
s. They prove very useful to better understand a forming process
(such as injection or compression moulding for instance). The
time at which the material is crystallized to the core is of particular interest. This corresponds to the time when the material is
fully solid. This is the end of the forming phase and the part can
eventually be ejected [11].
Discussion
The novel crystallization physics node proves useful to efficiently implement a coupled heat transfer and crystallization case.
Nonetheless, some numerical considerations should be discussed.
First, the above resolution is efficient and accurate because the
crystallization transition zone is wide compared to the mesh and
part size. In the case of thick or large geometries and/or faster
cooling cycles, crystallization fronts may appear. An accurate
modelling of such a sharp transition zone would require an automatic remeshing method, as described in [4]. This can be implemented in COMSOL but the team is waiting for future versions
to provide a better remeshing algorithm.
Then, the present test case considers a quiescent material. In the

Table of Contents for the Digital Edition of JEC COMPOSITES MAGAZINE - Issue #114 - July 2017