# Numerical Simulation of Spin Coating Processes Involving Functionalized Carbon Nanotube Suspensions

Elías Cueto, Amon Ma, Francisco Chinesta, M. R. Mackley

► To cite this version:

Elías Cueto, Amon Ma, Francisco Chinesta, M. R. Mackley. Numerical Simulation of Spin Coating Processes Involving Functionalized Carbon Nanotube Suspensions. International Journal of Material Forming, 2008, 2, pp.89-99. hal-01008567

HAL Id: hal-01008567

https://hal.science/hal-01008567

Submitted on 24 Apr 2017

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.

Public Domain
# Numerical simulation of spin coating processes involving functionalised Carbon nanotube suspensions

E. Cueto · A. W. K. Ma · F. Chinesta · M. R. Mackley

E. Cueto
Group of Structural Mechanics and Material Modelling,
Aragón Institute of Engineering Research, I3A,
Universidad de Zaragoza, Edificio Betancourt. María de Luna, 7 50012 Zaragoza, Spain
e-mail: ecueto@unizar.es
URL: http://gemm.unizar.es

A. W. K. Ma · M. R. Mackley
Department of Chemical Engineering,
University of Cambridge, New Museums Site,
Pembroke Street, Cambridge, CB2 3RA, UK

M. R. Mackley
e-mail: mrm5@cam.ac.uk
URL: http://www.cheng.cam.ac.uk/research/groups/polymer/

F. Chinesta
Laboratoire de Mecanique des Systemes et des Procedes,
UMR 8106, CNRS-ENSAM-ESEM, Paris, France
e-mail: francisco.chinesta@paris.ensam.fr
URL: http://www.paris.ensam.fr/lmsp/

rotary diffusion coefficient accounting for randomising events. A Natural Element strategy with an updated Lagrangian framework was implemented to solve the free-surface problem involving large domain deformation and to avoid numerical problems associated with Finite Element modelling. The model reported herein couples micro-scale CNT orientation with the macroscopic suspension kinematics and it offers important insights in relation to the final properties of spin-coated CNT films as well as the processing behaviour of CNT suspensions.

**Keywords** Carbon nanotubes · Spin coating · Numerical simulation · Meshless

**Abstract** This paper reports the numerical simulation of spin coating for functionalised Carbon Nanotube (CNT) suspensions. Spin coating is a process commonly used to deposit uniform thin films onto flat substrates by means of high rotation velocity and centrifugal force. The functionalised CNTs modelled in this study were chemically treated in a way such that aggregation was prevented through electrostatic repulsion between CNTs. The functionalised CNTs in the semidilute suspensions can be modelled as rigid fibres with their orientation dictated by the flow of the solvent. The evolution of CNT orientations was simulated using a pre-averaged kinetic theory with an appropriate

## Introduction

Carbon nanotubes have generated a tremendous interest in the scientific community following a journal publication by Iijima (1991). Single-walled carbon nanotubes are one-atom thick sheets of graphite (called graphene) rolled up into seamless cylinders with diameter in the order of a nanometer. They possess impressive mechanical properties:

- Young's modulus: 1 to 5 TPa

- Tensile Strength: 13-53 GPa

- Elastic strain up to 5%.

- Density 2160 kg/m³.

Given these properties, carbon nanotubes (CNTs) have been extensively researched as possible reinforcement (and conductive elements) for high-performance nanocomposites. In addition, nanotube based transistors have been made that operate at room temperature

Work partially supported by the Spanish Ministry of Education and Science through grant number CICYT-DPI2005-08727-C02-01.
and that are capable of digital switching using a single
electron (Postma et al. 2001).

In order to design manufacturing processes more
effectively, it is important to understand and model
the behaviour of CNT suspensions and the orientation
of CNTs. Although several studies have focused on
producing polymer nanotube composites, many prac-
tical challenges remain before their potential can be
fully realized. Dispersing the nanotubes individually
and uniformly into the matrix seems to be fundamental
in producing composites with reproducible and optimal
properties.

This work is aimed at establishing a method for
the numerical simulation of CNT suspension form-
ing processes. We focus our attention on spin coating
processes, although we believe that the technique is
general enough to be applied to other forming pro-
cesses as well. Spin coating processes are used to apply
very thin, more or less uniform, films onto flat sub-
strates. In short, a drop of solution is placed in the cen-
ter of the plate, which is then rotated at high speed. The
centrifugal force makes the solution to spread on the
flat substrate until the desired thickness is achieved on
the film. The solvent is sometimes volatile in this kind
of forming processes, an effect that is not considered in
this work.

From the simulation point of view, spin coating
processes involving CNTs suspensions present some
very important challenges. First of all, tracking the
evolution of the free surface of the liquid, that under-
goes very large deformations, remains a problem for
traditional, Eulerian or Arbitrary Lagrangian Eulerian
Finite Element techniques (Donea 1983), although
some early references are available on the literature
(Emslie et al. 1958). To this end, Lagrangian frame-
works, in which the nodes of the simulation move with
the material velocity, seem to be an appropriate choice
for the simulation. The family of meshless methods
(Belytschko et al. 1994; Liu et al. 1995; Sukumar et al.
1998) seem to be an appealing choice, since they do
not rely on a mesh (and consequently mesh distortion
does not affect the accuracy of the results). In this work
we employ a particular meshless method known as
the Natural Element method and the method was prov-
en to be capable of solving numerical problems that
involve large distortion (Martinez et al. 2004; González
et al. 2007).

The other remaining challenge is the constitutive
modelling of CNT suspensions. In a first attempt, func-
tionalized CNTs are considered as rigid fibres with very
large aspect ratio (between 10³⁻¹⁰⁴, possibly up to 10⁶,
but typically in the order of hundreds after disper-
sion), whose orientation is governed by the flow of the

solvent, which is Newtonian in this study. A quadratic
closure relation is then applied to the preaveraged
kinetic theory model. At the macroscopic scale, the
flow kinematics of the suspension is governed by the
momentum and mass conservation equations, which
are further coupled with the micro-scale orientation of
the CNTs. The CNT orientation and the suspension
kinematics are updated at each time step as discussed
in Section Discrete model.

Some experimental measurements have been per-
formed on the rheology of the CNTs suspensions, in or-
der to characterize its behaviour. The CNTs modelled
in this paper were chemically treated such that they
did not aggregate and no mesostructures were observed
optically.

The paper includes a brief description of the Natural Element Method as well as the main assumptions involved the numerical modelling of spin coating. Some numerical examples are given to illustrate the potential of the proposed process simulation technique.

Spin coating processes

As mentioned before, spin coating processes consist of
making a drop of liquid spread on a surface by means of
centrifugal acceleration. This acceleration is achieved
using a disk that rotates at a high angular velocity, ω,
see Fig. 1.

The difficulties in the numerical simulation of spin
coating processes arise from (1) the tremendous defor-
mation of the domain (in spin coating, the boundary
can deform from a liquid droplet to thin films with a
thickness of 30 – 50 µm) and (2) the multiscale nature
of the problem.

In this work, the rotating plate is assumed to be
perfectly horizontal, so that there is no radial com-
ponent of gravity. Following Emslie et al. (1958), we
assume that Coriolis forces are negligible, compared
to centrifugal forces. To model the geometry of the
drop, axial symmetry around the rotation axis is also
assumed. A system of reference that rotates with the
disc was used in describing the equations of motion.

<figure><img src="image_1.png" /><figcaption>**Fig. 1** Schematic representation of spin coating processes</figcaption></figure>
## Numerical modelling of the spin coating process

CNT suspensions are modelled as short fibre suspensions described by the following equations (see Martínez et al. 2003, and references therein)

- The balance of momentum equations, where we only consider the centrifugal forces as in the treatment of Emslie et al. (1958)

$$ \mathrm{Div}\sigma = -\rho\omega \times (\omega \times \mathbf{r}) = \mathbf{f}_r \quad (1) $$

where $\sigma$ is the stress tensor, $\rho$ represents the density, $\times$ denotes the tensor product, and $\omega$ the rotation speed of the spin coater as schematically drawn in Fig. 1.

- The incompressibility condition

$$ \mathrm{Div}\mathbf{v} = 0 \quad (2) $$

where $\mathbf{v}$ represents the velocity field.

- The constitutive equation, with a quadratic closure relation for the fourth order orientation tensor and other simplifying assumptions, results

$$ \sigma = -p\mathbf{I} + 2\eta\{\mathbf{D} + N_p\mathrm{Tr}(a\mathbf{D})a\} \quad (3) $$

where $p$ denotes the pressure, $\mathbf{I}$ is the unit tensor, $\eta$ is the equivalent suspension viscosity, $\mathbf{D}$ is the strain rate tensor, $N_p$ a scalar parameter that depends on both the tube concentration and its aspect ratio, and $a$ is the second order orientation tensor defined by

$$ \mathbf{a} = \oint \rho \otimes \rho \Psi(\rho) d\rho \quad (4) $$

where $\rho$ is the unit vector defining the CNT axis direction, $\otimes$ denotes the dyadic product, and $\Psi(\rho)$ is the orientation distribution function, that should satisfy the normality condition

$$ \oint \Psi(\rho) d\rho = 1 \quad (5) $$

If $\Psi(\rho) = \delta(\rho - \hat{\rho})$, with $\delta(\cdot)$ the Dirac's function, all the orientation probability is concentrated in the direction defined by $\hat{\rho}$, and the corresponding orientation tensor is $\hat{\mathbf{a}} = \hat{\rho} \otimes \hat{\rho}$. In this case the quadratic closure relation becomes exact. However, when the tubes are not perfectly aligned in a certain direction the quadratic closure is no more exact and it constitutes an approximation that is widely used in the context of short fibre suspensions, see Martínez et al. (2003) and references therein.

From a physical point of view, we can consider that the eigenvalues of the second order orientation tensor ($\mathbf{a}$) represent the probability of finding the CNT in the direction of the corresponding eigenvectors.

- With a quadratic closure relation the orientation equation is expressed as

$$ \frac{da}{dt} = \Omega a - a\Omega + k(\mathbf{D}a + a\mathbf{D} - 2\mathrm{Tr}(a\mathbf{D})a) - 6D_r \left( a - \frac{\mathbf{I}}{3} \right) \quad (6) $$

$\mathbf{D}$ and $\Omega$ are the symmetric and skew-symmetric components of $\mathrm{Grad}\mathbf{v}$, $k$ is a constant that depends on the nanotube aspect ratio $r$ (fibre length to fibre diameter ratio): $k = (r^2 - 1)/(r^2 + 1)$ with $k \approx 1$ given the high aspect ratio of the CNTs. $D_r$ is a rotary diffusion coefficient accounting for randomising events such as Brownian motion and fibre-fibre interaction in the suspension.

The flow model is defined in the volume occupied by the fluid at time $t$, $\Omega_f(t)$. On its boundary, $\Gamma_f(t) \equiv \partial\Omega_f(t)$ either the velocity or the traction is imposed:

$$ \mathbf{v}(\mathbf{x} \in \Gamma_1) = \mathbf{v}_g \quad (7) $$

or

$$ \sigma \mathbf{n}(\mathbf{x} \in \Gamma_2) = \mathbf{F}_g \quad (8) $$

with $\Gamma_1 \cup \Gamma_2 = \Gamma_f(t)$, $\Gamma_1 \cap \Gamma_2 = \emptyset$, and where $\mathbf{n}(\mathbf{x})$ is the unit outwards vector, defined on the boundary at the point $\mathbf{x}$.

In terms of the initial condition, the CNTs are assumed to be isotropically orientated (i.e., $\mathbf{a} = \mathbf{I}/3$). This model is solved on the domain defined by the nodes as they evolve. More details can be found at Martínez et al. (2003) and references therein.

In this modelling, the following form of the orientation vector is used for describing CNTs suspended within the solvent:

$$ \rho = \begin{pmatrix} r \\ \theta \\ z \end{pmatrix} \quad (9) $$

with $\theta \in [0, 2\pi]$.

The velocity field is expressed in cylindrical coordinates, after the imposition of axial symmetry:

$$ \mathbf{v} = \begin{pmatrix} v_r \\ \omega \cdot r \\ v_z \end{pmatrix} \quad (10) $$

where $\omega$, as mentioned before, represents the angular velocity of the spin coater. If, following Emslie et al. (1958), we assume a system of reference rotating with the disk, we arrive at

$$ \mathbf{v} = \begin{pmatrix} v_r \\ 0 \\ v_z \end{pmatrix} \quad (11) $$
In this case, tensors **D** and **Ω** come from an axisym-
metric flow, so their expression in cylindrical coordi-
nates, as mentioned before, gives:

$$
\mathbf{D} = \begin{pmatrix}
\frac{\partial v_r}{\partial r} & \frac{1}{2} \left( \frac{1}{r} \frac{\partial v_r}{\partial \theta} + \frac{\partial v_\theta}{\partial r} - \frac{v_\theta}{r} \right) & \frac{1}{2} \left( \frac{\partial v_r}{\partial z} + \frac{\partial v_z}{\partial r} \right) \\
\frac{1}{2} \left( \frac{1}{r} \frac{\partial v_r}{\partial \theta} + \frac{\partial v_\theta}{\partial r} - \frac{v_\theta}{r} \right) & \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} + \frac{v_r}{r} & \frac{1}{2} \left( \frac{1}{r} \frac{\partial v_z}{\partial \theta} + \frac{\partial v_\theta}{\partial z} \right) \\
\frac{1}{2} \left( \frac{\partial v_r}{\partial z} + \frac{\partial v_z}{\partial r} \right) & \frac{1}{2} \left( \frac{1}{r} \frac{\partial v_z}{\partial \theta} + \frac{\partial v_\theta}{\partial z} \right) & \frac{\partial v_z}{\partial z}
\end{pmatrix}, \quad (12)
$$

and, provided that the flow is axially symmetric ($v_{\theta} = 0$,
$v_{r} = v_{r}(r, z)$, $v_{z} = v_{z}(r, z)$), and by combining the above
equation with Eq. 11 **D** reduces to

$$
\mathbf{D} = \begin{pmatrix}
\frac{\partial v_r}{\partial r} & 0 & \frac{1}{2} \left( \frac{\partial v_r}{\partial z} + \frac{\partial v_z}{\partial r} \right) \\
0 & \frac{v_r}{r} & 0 \\
\frac{1}{2} \left( \frac{\partial v_r}{\partial z} + \frac{\partial v_z}{\partial r} \right) & 0 & \frac{\partial v_z}{\partial z}
\end{pmatrix} . \quad (13)
$$

Similarly, the vorticity tensor takes the following
form:

$$
\Omega = \begin{pmatrix}
0 & 0 & \frac{1}{2} \left( \frac{\partial v_r}{\partial z} - \frac{\partial v_z}{\partial r} \right) \\
0 & 0 & 0 \\
-\frac{1}{2} \left( \frac{\partial v_r}{\partial z} - \frac{\partial v_z}{\partial r} \right) & 0 & 0
\end{pmatrix} \quad (14)
$$

An isotropic orientation is defined in this model by
the uniform distribution

$$
\psi(\rho) = \frac{1}{4\pi} \tag{15}
$$

giving rise to an orientation tensor, for isotropic
distributions,

$$
a = \frac{1}{3} \mathbf{I} \tag{16}
$$

These conditions are imposed at the initial time step
only.

**Rheology of functionalized CNTs suspensions**

A Fokker-Planck (FP) based orientation model was
used to describe the steady shear rheological responses
of the chemically treated CNT suspensions. Figure 2a
shows the model fitting of surface-treated CNT sus-
pensions with three different weight concentrations
(0.05%, 0.2% and 0.33%). Taking the 0.3% CNT sus-
pension as an example, the evolution of the apparent

viscosity ηₐ in the orientation model is controlled by
the values of Dₚ and Nₚ and the best constant- Dₚ
fit was obtained for Dₚ = 0.005s⁻¹ and Nₚ = 7. In
Fig. 2b, different values of Dₚ were used to illus-
trate the sensitivity of the model on the fitting pa-
rameter Dₚ. In general, there is a good agreement

<figure><img src="image_2.png" /><figcaption>**Fig. 2** a Orientation model fitting of 0.05%, 0.2% and 0.33% surface-treated CNTs suspended in epoxy resin. N<sub>p</sub> and D<sub>r</sub> are the fitting parameters. b Sensitivity of the model to D<sub>r</sub> given the experimental data of 0.33% treated CNT suspension and N<sub>p</sub> = 7</figcaption></figure>
between $\eta_a$ predicted by the model and the experimental data, but for high $\dot{\gamma}_R$, the predicted $\eta_a$ is slightly higher than the experimental value. For instance, at shear rate $\dot{\gamma}_R = 60s^{-1}$, $\eta_a$ of the 0.3% suspension was found experimentally to be $10.5 Pa \cdot s$, but the orientation model predicted a viscosity of $11.3 Pa \cdot s$ (with a constant $D_r = 0.005s^{-1}$) and this gives an error of about 8%. Although there is an error of a few percents in predicting high-shear-rate data, the use of a constant $D_r$ should be sufficient to provide a reasonably good estimation of $\eta_a$ for general engineering problems, and this also saves the effort in determining the exact relationship between $D_r$ and $\dot{\gamma}_R$. However, in cases where a more accurate description of high-shear viscosity is needed, the dependence of $D_r$ on should be carefully evaluated (see for example, Folgar and Tucker 1984; Larson 1999). Moreover, our recent modelling work suggested that the FP based orientation model cannot satisfactorily describe the shear-thinning characteristics for untreated CNT suspensions where optically resolvable CNT aggregates are present (see for example Rahatekar et al. 2006). For untreated CNT suspensions, a new aggregation/orientation model should be used and in that model, diffusion coefficient depends not only on the orientation of CNTs, but also on the entanglement state of CNTs. Detailed model formulation and fittings for both treated and untreated CNT suspensions can be found in the article entitled "The Rheological Modelling of CNT Suspensions" in this special issue.

## Basics of the natural element method

As mentioned in the introductory section, we have chosen to employ a meshless method, the Natural Element Method, to perform the simulation, mainly because it allows for a Lagrangian description of the flow kinematics, convenient for solving free-surface flow problems with large surface deformation. The basics of the method are described next.

<figure><img src="image_3.png" /><figcaption>**Fig. 3** Delaunay triangulation and Voronoi diagram of a cloud of points</figcaption></figure>

## Natural neighbour interpolation

The vast majority of meshless methods are based on the use of scattered data approximation techniques to construct the approximating spaces of the Galerkin method. These techniques must have, of course, low sensitivity to mesh distortion, as opposed to Finite Element methods. Among these techniques, the Natural Element Method employs any instance of Natural Neighbour interpolation (Sibson 1981; Hiyoshi and Sugihara 1999) to construct trial and test functions. Prior to the introduction of these interpolation techniques, it is necessary to define some basic concepts.

The model will be constructed upon a cloud of points with no connectivity established a priori among them. We will call this cloud of points $N = \{n_1, n_2, ..., n_M\} \subset \mathbb{R}^d$, and there is a unique decomposition of the space into regions such that each point within these regions is closer to the node with which the region is associated than to any other in the cloud. This kind of space decomposition is called a Voronoi diagram of the cloud of points and each Voronoi cell is formally defined as (see Fig. 3):

$$T_I = \{ \mathbf{x} \in \mathbb{R}^d : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) \quad \forall J \neq I \}, \quad (17)$$

where $d(\cdot, \cdot)$ is the Euclidean distance function.

The dual structure of the Voronoi diagram is the Delaunay triangulation, obtained by connecting nodes that share a common ($d-1$)-dimensional facet. While the Voronoi structure is unique, the Delaunay triangulation is not, there exist some degenerate cases in which there are two or more possible Delaunay triangulations (consider, for example, the case of triangulating a square in 2D, as depicted in Fig. 3 (right)). Another way to define the Delaunay triangulation of a set of nodes is by invoking the *empty circumcircle* property, which means that no node of the cloud lies within the circle covering a Delaunay triangle. Two nodes sharing a facet of their Voronoi cell are called *natural neighbours* and hence the name of the technique.
In order to define the natural neighbour co-ordinates it is necessary to introduce some additional concepts. The second-order Voronoi diagram of the cloud is defined as

$$ T_{IJ} = \{ \mathbf{x} \in \mathbb{R}^d : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) < d(\mathbf{x}, \mathbf{x}_K) \\ \forall J \neq I \neq K \}. \quad (18) $$

The most common natural neighbour interpolation method, is the Sibson interpolant (Sibson 1980, 1981). Consider the introduction of the point $\mathbf{x}$ in the cloud of nodes. Due to this introduction, the Voronoi diagram will be altered, affecting the Voronoi cells of the natural neighbours of $\mathbf{x}$. Sibson (1980) defined the natural neighbour coordinates of a point $\mathbf{x}$ with respect to one of its neighbours $I$ as the ratio of the cell $T_I$ that is transferred to $T_x$ when adding $\mathbf{x}$ to the initial cloud of points to the total volume of $T_x$. In other words, if $\kappa(\mathbf{x})$ and $\kappa_I(\mathbf{x})$ are the Lebesgue measures of $T_x$ and $T_{xI}$ respectively, the natural neighbour coordinates of $\mathbf{x}$ with respect to the node $I$ is defined as

$$ \phi_I(\mathbf{x}) = \frac{\kappa_I(\mathbf{x})}{\kappa(\mathbf{x})}. \quad (19) $$

In Fig. 4 the shape function associated to node 1 may be expressed as

$$ \phi_I(\mathbf{x}) = \frac{A_{abfe}}{A_{abcd}}. \quad (20) $$

It is straightforward to prove that natural element (NE) shape functions (see Fig. 5) form a partition of unity

<figure><img src="image_4.png" /><figcaption>Fig. 4 Definition of the natural neighbour coordinates of a point $\mathbf{x}$</figcaption></figure>

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

<figure><img src="image_6.png" /><figcaption>Fig. 5 Typical function $\phi(\mathbf{x})$</figcaption></figure>

(Babuška and Melenk 1997), as well as some other properties like positivity (i.e., $0 \le \phi_I(\mathbf{x}) \le 1 \forall I, \forall \mathbf{x}$) and strict interpolation:

$$ \phi_I(\mathbf{x}_J) = \delta_{IJ}. \quad (21) $$

Sibson interpolants have some remarkable properties that help to construct the trial and test functional spaces of the Galerkin method (see Sukumar et al. 1998; Hiyoshi and Sugihara 1999).

Besides properties like continuity and smoothness (everywhere except at the nodes for Sibson interpolants), Sibson interpolants possess linear completeness (i.e., exact reproduction of a linear field) and it is noteworthy that they are able to exactly interpolate prescribed essential (Dirichlet) boundary conditions, a condition usually not fulfilled by other meshless methods.

## The $\alpha$-shapes-based NE method

The identification of the free surface in an updated Lagrangian flow simulation deserves some comments. Using meshless methods, in which models are constructed by a set of nodes only, boundary tracking can be performed by employing different strategies. In particular, we have employed *shape constructors* to perform this task. Shape constructors are geometrical entities that transform finite point sets into a multiply connected shape in general. In particular, we employ $\alpha$-shapes (Edelsbrunner and Mücke 1994). $\alpha$-shapes define a one-parameter family of shapes $S_\alpha$ (being $\alpha$ the parameter), ranging from the “coarsest” to the “finest” level of detail. $\alpha$ can be seen, precisely, as a measure of this level of detail.
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Evolution of the family of α-shapes of a cloud of points representing drop impacting on a flat surface. Shapes $S_{0.03}$ (a), $S_{0.05}$ (b), $S_{0.5}$ (c), $S_{2.0}$ (d), $S_{100.0} \sim S_\infty$ (e) are depicted</figcaption></figure>

Details about the formal definition of the family of α-shapes can be found in Edelsbrunner and Mücke (1994). In brief, the use of α-shapes to define the boundary of the domain relies in the choice of the level of detail needed to represent the domain, which is always an analyst's decision. It is obvious that the minimum nodal spacing parameter, say $h$, should be chosen so as to reproduce at least that level of detail $\alpha$.

α-shapes provide a means to eliminate from the triangulation those triangles or tetrahedra whose size is bigger than the specified level of detail. This criterion is

very simple: just perform a filtration to eliminate those triangles (tetrahedra) whose circum-radius is bigger than the level of detail, $\alpha$.

In Fig. 6 an example of the previously presented theory is presented. It represents some instances of the finite set of shapes for a cloud in a intermediate step of the simulation of drop impacting onto a flat surface. The experimental solution of this problem, together with some snapshots of its numerical simulation, are shown in Figs. 7 and 8.

Note that the key question in using α-shapes is not to find the precise value of $\alpha$ for a given configuration of the nodal cloud. Instead, we must set the problems in terms of *what level of detail we are interested in taking into account* for a particular geometry.

The use of α-shapes, however, has another relevant influence in the NE Method. As demonstrated in Cueto et al. (2000), the construction of natural neighbour interpolation on an α-shape of the domain alters the distance measure. Natural neighbour interpolation is performed on the basis of Voronoi diagrams, which employ euclidean distance measure in their most general form. This leads to some lack of interpolation along non-convex boundaries. This interpolation is recovered if we construct the natural neighbour interpolants over an α-shape of the domain.

Thus, the use of α-shapes in the construction of updated Lagrangian simulations of fluid flow provides an appealing way to track the boundary of domain while ensuring appropriate interpolation of essential boundary conditions, that can be imposed directly in the discrete system of equations, as in the Finite Element Method.

## Discrete model

A mixed $C^0 - C^{-1}$ NE interpolation has been applied to discretize the strong form of the problem introduced before. More details on the application of mixed NE interpolations in incompressible and nearly incompressible media can be found in Cueto et al. (2001).

In the formulation here presented, a $C^0$ interpolation scheme – smooth everywhere except at the nodes – has been chosen in the velocity field approximation,

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Sequence of the drop deformation under very low ambient pressure. Photos courtesy of Lei Xu and Sidney Nagel</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8 Evolution of the free surface in the drop impact</figcaption></figure>

whereas a discontinuous $C^{-1}$ interpolation has been used in the pressure approximation:

$$v^h(\mathbf{x}) = \sum_{I=1}^{n} \phi_I(\mathbf{x}) v_I \qquad (22)$$

$$p^h(\mathbf{x}) = \sum_{I=1}^{n} \psi_I(\mathbf{x}) p_I = \sum_{I=1}^{n} \frac{1}{n} p_I \qquad (23)$$

where $v_I$ and $p_I$ represent the nodal velocities and pressures, respectively, and $n$ is the number of natural neighbours of the point $\mathbf{x}$ under consideration.

This kind of approximation does not verify the LBB condition. However, it has been shown that its behaviour is very similar to that of the bilinear velocity-constant pressure finite element (see González et al. 2004a). No spurious modes or locking have been observed in all the simulations computed in this paper.

Thus, with the fluid domain $\Omega_f(t)$ extracted at time $t$ from the cloud of nodes by using the $\alpha$-shape technique, as described in the previous section, and the velocity and pressure NE interpolation defined by Eqs. 22 and 23, we can proceed to a standard discretisation of the mixed variational formulation of the flow equations

$$\int_{\Omega_f(t)} \sigma : \mathbf{D}^* d\Omega = \int_{\Omega_f(t)} f_r v^* d\Omega \qquad (24)$$

$$\int_{\Omega_f(t)} \mathrm{Div} v p^* d\Omega = 0 \qquad (25)$$

with

$$\sigma = -p\mathbf{I} + 2\eta \{\mathbf{D} + N_p \mathrm{Tr}(a\mathbf{D})a\} \qquad (26)$$

where a null traction is assumed on the flow front and a prescribed velocity is enforced on the other part of the fluid domain boundary. Essential boundary conditions involve setting the velocity of the nodes in contact with the rotating plate to zero.

The orientation equation is solved at each time increment. With the flow kinematics known at time $t$, $v^t(\mathbf{x})$, the position of nodes can be updated at the same time that the fibres orientation evolution, given by Eq. 6 is computed by using the method of characteristics. The simplest explicit updating consists in writing

$$\mathbf{x}_I^{t+\Delta t} = \mathbf{x}_I^t + \mathbf{v}_I^t \Delta t, \quad \forall I \qquad (27)$$

and

$$\mathbf{a}_I^{t+\Delta t} = \mathbf{a}_I^t + \left\{ \begin{aligned} & \mathbf{\Omega}_I^t \mathbf{a}_I^t - \mathbf{a}_I^t \mathbf{\Omega}_I^t + k \mathbf{D}_I^t \mathbf{a}_I^t + k \mathbf{a}_I^t \mathbf{D}_I^t \\ & - 2k \mathrm{Tr}(\mathbf{a}_I^t \mathbf{D}_I^t) \mathbf{a}_I^t - 6D_r \left( \mathbf{a}_I^t - \frac{\mathbf{I}}{3} \right) \end{aligned} \right\} \Delta t, \quad \forall I \qquad (28)$$

where $\mathbf{D}_I^t$ and $\mathbf{\Omega}_I^t$ are the symmetric and skew-symmetric components of the gradient of velocity tensor, at time $t$ in the node $\mathbf{x}_I$, respectively.

Finally, in the kinematics resolution stage, the fibre orientation described by the second order orientation tensor $\mathbf{a}$ is assumed to be known at the nodes at present time step, $\mathbf{a}_I^t$. The value of $\mathbf{a}$ at the integration points
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Snapshots of the velocity field at 1st (a), 10th (b), 20th (c) and 30th (d) time steps. $\Delta t = 0.01s$</figcaption></figure>

used to evaluate Eqs. 24 and 25 is computed by using the NE interpolation

$$ \mathbf{a}^t(\mathbf{x}) = \sum_{I=1}^{I=n} \phi_I(\mathbf{x}) \mathbf{a}_I^t \qquad (29) $$

The only difficulty in applying Eq. 28 to update the orientation is related to the non-derivability of the NE shape functions at their definition nodes. Thus, we can evaluate the velocity gradient tensor from the expression

$$ \mathbf{v}_{\cdot,k}^t(\mathbf{x}) = \sum_{I=1}^{I=n} \phi_{I,k}(\mathbf{x}) \mathbf{v}_I^t \qquad (30) $$

where the subscripts k denotes the spatial derivative with respect to the k-coordinate. The velocity gradient can be computed everywhere except at the nodes, since $\phi_{I,k}(\mathbf{x}_I)$ is not defined. A possible solution for this relies in the use of Stabilised Conforming Nodal integration techniques (Chen et al. 2001; Gonzalez et al. 2004b) or to use some kind of projection from integration points

and then averaging (Zienkiewicz and Zhu 1992; Cueto et al. 2003).

## Numerical results

The previously presented model was applied to the simulation of spinning a drop of CNT suspension, with the assumption that the initial shape is Gaussian. The model was composed of 977 nodes under axisymmetric assumptions for the flow (but not for the orientation field, which is three-dimensional). Non-slip boundary conditions were assumed at the horizontal plane, whereas symmetry boundary conditions are assumed at the axis of symmetry. Some snapshots for the simulated evolution of the velocity field for different time steps are shown in Fig. 9.

The results on the orientation field for the CNTs are shown in Figs. 10 and 11. It can be noticed how CNTs tend to adopt an horizontal orientation in the r-z plane, due to the gradient of velocities in this plane, very close to pure shear due to the non-slip boundary
<figure><img src="image_11.png" /></figure>

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

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

<figure><img src="image_14.png" /><figcaption>**Fig. 10** Orientation field for the CNT suspension at some intermediate steps of the simulation (r-z orientation field). 1st (a), 10th (b), 20th (c) and 30th (d) time steps</figcaption></figure>

conditions. On the θ-z plane, however, and by virtue of the diverging velocity field, CNTs tend to align to the θ direction, although the ellipses are not as distorted as those on the r-z plane and the degree of alignment is less significant.

## Conclusions

A model has been presented for the numerical simulation of spin coating of CNTs suspensions. The model is composed of a meshless (NE) strategy for the flow that allows for a proper description of the free-surface flow and a micro-macro description of the nanotube scale that allows for a suitable description of CNT orientation field.

The use of a meshless strategy in an updated Lagrangian framework, combined with an α-shape approach, to describe the large motion of the free surface seems to be an appealing choice, due to its simplicity. In this way since nodes move with material velocity, the orientation field can be stored at the nodes and be updated easily by a simple algorithm based on the method of characteristics.

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

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

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

<figure><img src="image_18.png" /><figcaption>**Fig. 11** Orientation field for the CNT suspension at some intermediate steps of the simulation (r-θ plane). 1st (a), 10th (b), 20th (c) and 30th (d) time steps</figcaption></figure>

Other possibilities, such as the kinetic description derived from the FP equation or its equivalent stochastic (Itô) counterpart could also be applied in this same framework. This remains an important part of ongoing research and the findings will be reported in a future paper.

## References

* Babuška I, Melenk JM (1997) The partition of unity method. Int J Numer Methods Eng 40:727-758

* Belytschko T, Lu YY, Gu L (1994) Element-free galerkin methods. Int J Numer Methods Eng 37:229-256

* Chen JS, Wu CT, Yoon S, You Y (2001) A stabilized conforming nodal integration for galerkin mesh-free methods. Int J Numer Methods Eng 50:435-466

* Cueto E, Doblaré M, Gracia L (2000) Imposing essential boundary conditions in the natural element method by means of density-scaled α-shapes. Int J Numer Methods Eng 49-4: 519-546

* Cueto E, Martínez MA, Doblaré M (2001) El método de los Elementos Naturales en Elasticidad compresible y cuasi-incompresible. Boletín Técnico del Instituto de Materiales y Modelos Estructurales. Universidad Central de Venezuela 39 (In Spanish)
Cueto E, Sukumar N, Calvo B, Martínez MA, Cegoñino J, Doblaré M (2003) Overview and recent advances in natural neighbour Galerkin methods. Arch Comput Methods Eng 10(4):307-384

Donea J (1983) Arbitrary Lagrangian-Eulerian finite element methods. In: Belytschko T, Hughes TJR (eds) Computer methods for transient analyses. Elsevier, Amsterdam, pp 474-516

Edelsbrunner H, Mücke E (1994) Three dimensional alpha shapes. ACM Trans Graph 13:43-72

Emslie AG, Bonner FT, Peck LG (1958) Flow of a viscous liquid on a rotating shell. J Appl Phys 29(5):858-863

Folgar F, Tucker CL III (1984) Orientation behaviour of fibers in concentrated suspensions. J Reinf Plast Compos 3:98-119

González D, Cueto E, Doblaré M (2004a) Volumetric locking in natural neighbour galerkin methods. Int J Numer Methods Eng 61(4):611-632

Gonzalez D, Cueto E, Martinez MA, Doblare M (2004b) Numerical integration in natural neighbour galerkin methods. Int J Numer Methods Eng 60(12):2077-2104

González D, Cueto E, Chinesta F, Doblaré M (2007) A natural element updated Lagrangian strategy for free-surface fluid dynamics. J Comput Phys 223(1):127-150

Hiyoshi H, Sugihara K (1999) Two generalizations of an interpolant based on Voronoi diagrams. Int J Shape Model 5(2):219-231

Iijima S (1991) Single-shell carbon nanotubes of 1-nm diameter. Nature 363:603-605

Larson RG (1999) The structure and rheology of complex fluids. Oxford University Press, Oxford

Liu WK, Jun S, Li S, Adee J, Belytschko T (1995) Reproducing kernel particle methods. Int J Numer Methods Eng 38:1655-1679

Martínez MA, Cueto E, Doblaré M, Chinesta F (2003) Fixed mesh and meshfree techniques in the numerical simulation of injection processes involving short fiber suspensions. J Non-Newton Fluid Mech 115:51-78

Martinez MA, Cueto E, Alfaro I, Doblare M, Chinesta F (2004) Updated Lagrangian free surface flow simulations with natural neighbour galerkin methods. Int J Numer Methods Eng 60(13):2105-2129

Postma HWC, Teepen T, Yao Z, Grifoni M, Dekker C (2001) Carbon nanotube single-electron transistors at room temperature. Science 293:5527

Rahatekar SS, Koziol KKK, Butler SA, Elliott JA, Mackley M, Windle AH (2006) Optical microstructure and viscosity enhancement for an epoxy resin matrix containing multi-wall carbon nanotubesorientation behaviour of fibers in concentrated suspensions. J Rheol 50:599-610

Sibson R (1980) A vector identity for the dirichlet tesselation. Math Proc Camb Philos Soc 87:151-155

Sibson R (1981) A brief description of natural neighbour interpolation. In: Barnett V (ed) Interpreting multivariate data. Wiley, New York, pp 21-36

Sukumar N, Moran B, Belytschko T (1998) The natural element method in solid mechanics. Int J Numer Methods Eng 43(5):839-887

Zienkiewicz OC, Zhu JZ (1992) The superconvergent patch recovery and a posteriori error estimates. Part I: the recovery technique. Int J Numer Methods Eng 33:1331-1364