Mathematical modeling and numerical simulations of the motion of nanoparticles in straight tube

Abstract

Nanotechnology is a very important field in science and technology, and research projects in this domain attract considerable funding. The existing research works in nanotechnology deal with chemical, physical, and biological issues or a combination of these fields, but very small number of works has been undertaken on mathematical modeling. Mathematical models can greatly reduce the time involved in experimentation, which in turn reduces the research cost. In this article, we consider the mathematical modeling of the motion of nanoparticles in a viscous flow inside straight tube. Illustrative simulations of the model are provided.

Introduction

In this article, we are working on the modeling and simulations of some important phenomena raised in biotechnology and medicine at the nanoscale. We focus on the modeling and simulations of the motion of nanoparticles in artificial straight microtubes. The case of human capillaries (see Figure 1) and more complex artificial capillaries will be the second phase of our research work (forthcoming papers).

Figure 1.

Examples of human capillaries.

Our motivations come from the great impact of nanotechnology and biotechnology on pharmacology for improving performance of many drugs and allowing the use of new drugs and new therapies. Also, nanoparticles can be manufactured and used as drug delivery systems to (1) deliver drugs specifically at places where diseased cells are, and optimally control drug release rates—preprogrammed, self-regulated, and remotely controlled delivery vehicles; (2) locate diseased cells or tumoral masses and estimate the state of the disease; and (3) carry diagnostic agents (fluorescence molecules) to enhance imaging.

Based on these motivations, we started working on the following research direction: modeling and simulations of the motion of nanoparticles in straight tubes (the case of general shapes will be the subject of forthcoming papers) and taking into account all the interactions between them and with walls and obstacles. On this direction, we obtained satisfactory numerical simulations in the simplest prototype of straight circular capillaries (see Figure 2).

Figure 2.

Straight circular capillaries.

The main task of the mathematical model used in our analysis is the non-overlapping of the nanoparticles and the obstacles inside the channel during the simulation. In a forthcoming paper, we will treat more complex situations in which we will take into account all the possible forces governing the flow inside the micro-channel (for instance, the forces exerted by the blood (or liquid) stream, gravitation, and electromagnetic interactions) and all the possible forces between nanoparticles and surfaces (for instance, Van der Waals, steric-polymer forces, electrostatic interaction (repulsive and attractive), hemodynamic forces, buoyancy, etc). Also, the shape of the channel will be taken in a more general form, that is, the walls will be a function $h(z)$ of the distance to the axis of the channel. Based on our expected numerical simulations, we will produce nanochannel with different shapes and test the obtained results and then optimize the design of the produced nanochannels.

The present work is organized as follows. The section “Notation” is stated in Appendix I and it contains the notation that will be used throughout the paper. In section “Non-overlapping algorithm,” we build and analyze the algorithm that will be used to ensure the non-overlapping of the nanoparticles. Section “Dynamical system of the motion of nanoparticles” is devoted to derive the dynamical system governing the motion of nanoparticles inside the micro-channel. In section “Illustrative simulation,” we present our satisfactory numerical simulations.

Non-overlapping algorithm

In this section, we present a numerical scheme that can be used to handle contacts between rigid disks. This numerical scheme will be our main powerful tool to simulate without overlapping the motion of nanoparticles. This algorithm has been introduced by Maury ¹ for studying inelastic collision problems and used by Maury and Venel ² for the modeling of microscopic crowd motion. For another application of this algorithm to other practical problems (to obstacle avoidance of mobile robots), we refer the reader to Hedjar and Bounkhel. ³

Consider N nanoparticles identified to inelastic disks with different radii $R_i,i=1,\dots ,N$ . The aim is to produce a numerical simulation of the motion of nanoparticles inside a given domain $\Omega \subset {\mathbb{R}}^2$ (in the present work, $\Omega$ is a straight tube) during a given interval time $[0,T]$ and subject to obstacles (static and dynamic).

We consider the configuration vector $q=(q_1,q_2,\dots ,q_N)\in {\mathbb{R}}^{2N}$ where $q_i$ is the center of the disk i with coordinates $(x_i;y_i{)}^T$ . To ensure the non-overlapping of nanoparticles and obstacles, we consider the following set of admissible configurations

$$Q=\{q\in {\mathbb{R}}^{2N}:D_{\mathrm{ij}}(q)\ge 0,\forall i,j\}$$

(1)

where $D_{\mathrm{ij}}(q)=\parallel q_i-q_j\parallel -(R_i+R_j)$ is the distance between the inelastic disks i and j. Note that the mathematical interpretation of the non-overlapping of the inelastic disks i and j at time t is the non-negativity of the distance $D_{\mathrm{ij}}(t)$ . Starting from a given admissible configuration $q(t)$ at time t, we wish to guarantee that $q(t+h)\in Q$ after a small value of unit of time $h>0$ , that is, the new configuration vector $q(t+h)$ is admissible. For this purpose, we define the admissible velocity set

$$C_h(q)=\{v\in {\mathbb{R}}^{2N}:D_{\mathrm{ij}}(q)+h\nabla D_{\mathrm{ij}}(q)·v\ge 0,\mathrm{for}\mathrm{all}i<j\}$$

(2)

Here, $\nabla D_{\mathrm{ij}}(q)$ denotes the gradient vector of $D_{\mathrm{ij}}(q)$ . Assume that $q_n$ is a given admissible configuration at time $t_n\in [0,T]$ , that is, $q_n:=q(t_n)\in Q$ . The next configuration $q_{n+1}$ after h unit of time is $q_{n+1}=q(t_n+h)$ . Using the first-order approximation of the distance function $D_{\mathrm{ij}}(·)$ , we have

$$D_{\mathrm{ij}}(q(t_n+h))=D_{\mathrm{ij}}(q(t_n))+h\nabla D_{\mathrm{ij}}(q_n)·\stackrel{·}{q}(t_n)$$

(3)

Taking the real velocity at that time $t_n$ lying in $C_h(q_n)$ yields

$$D_{\mathrm{ij}}(q_n)+h\nabla D_{\mathrm{ij}}·\stackrel{·}{q}(t_n)\ge 0$$

and hence $q(t_n+h)\in Q$ .

Let $U(q)=(v_1,v_2,\dots ,v_{N_1})$ (with $v_i=(v_{x_i},v_{y_i})$ ) be the velocity of all nanoparticles obtained after solving a system of differential equations (see section “Dynamical system of the motion of nanoparticles” for the the construction of such system). $U(q)$ is called the desired velocity of all nanoparticles. To ensure the non-overlapping of all the nanoparticles and obstacles, we consider the following iterative scheme:

Assume that $q_n\in Q$ , that is, the configuration vector $q_n$ at time $t_n$ is admissible. Then, the next configuration $q_{n+1}$ will be defined by the following relation

$$q_{n+1}:=q_n+h{w}^{n+1}$$

(4)

where $w^{n+1}$ is the projection of the desired velocity at time $t_n$ over the admissible velocities set at the same time, that is, $w^{n+1}\in C_h(q_n)$ and satisfies the convex constrained minimization problem $({\mathcal{P}}_n)$

$$\parallel w^{n+1}-U(q_n){\parallel }^2=\mathrm{Mi}{n}_{v\in C_h(q_n)}\parallel v-U(q_n){\parallel }^2$$

The numerical scheme (4) can be seen geometrically as $w^{n+1}$ is the nearest direction to the desired velocity $U(q_n)$ in the admissible velocity set $C_h(q_n)$ ensuring the non-overlapping of nanoparticles and ensuring that the real velocities will be very close to the desired velocities. It is very important to note that the convex minimization problem ${\mathcal{P}}_n$ can be solved numerically by any suitable solvers for the determination of the value of $w^{n+1}$ . In this article, we combine the previous algorithm with a kinematic model of nanoparticles (that we describe in the next section) to produce a numerical simulation of a given huge number N of nanoparticles with the main task avoiding, during the time interval, the overlapping of nanoparticles and the avoidance of some given obstacles.

Remark 1

The main characterizations of the above algorithm are as follows: as explained before, the first one is the non-overlapping between the disks, which is ensured by the setup of the algorithm. However, the second characterization, which is the most requested in numerical analysis, is the continuity of the trajectory of the nanoparticles generated by the algorithm. This property of smoothness has been proved by Maury and Venel ² (see also section 3 in Hedjar and Bounkhel ³ ).

Dynamical system of the motion of nanoparticles

In this section, we are going to build the system of differential equations of the motion of nanoparticles inside a straight tube.

Consider the motion of N nanoparticles in a viscous flow inside a tube of width a and length l and bounded by two horizontal walls as shown in Figure 4. Let $(x,y,z)$ be a right-handed system of rectangular Cartesian coordinates, of which the x axis points to the axis of the tube and has the direction of the flow, and the z axis points to the fluid normal to the wall as shown in Figure 3. For simplicity, we assume that the motion of the nanoparticle is constrained to the xz-plane, that is, the nanoparticle center is always in the $y=0$ plane, and there is no rotation of the nanoparticle about the z-axis. The nanoparticle, therefore, has two degrees of freedom: its motion is fully described by specifying $(x_i,z_i)$ as functions of time t, where $x_i$ and $z_i$ are coordinates of the center of nanoparticle i.

Figure 3.

Nanoparticles inside a straight tube.

The trajectory of each nanoparticle is governed by the forces exerted by the flow (blood stream) and the gravitation. Forces acting on nanoparticles include hemodynamic forces, buoyancy, Van der Waals interactions between nanoparticles each-to-other and between nanoparticles and walls.

To derive the equations describing the motion on the x-axis and on the z-axis, we apply (as suggested by Zhao et al. ⁴ ) the balance principle of forces acting on the nanoparticles. ⁵ The balance of forces acting on the nanoparticles requires, for each nanoparticle i, a system of two equations on the x-axis and the z-axis. So, on the z-axis, we have, for any nanoparticle i, the following equation

$$[h_5^{w_1}(z_i)+h_5^{w_2}(z_i)]·v_{z_i}-\sum _{j\ne i}{F}_{\mathrm{ij}}-F_{\mathrm{iw}}=\frac{4\pi }{3}\mathrm{gR}\Delta \rho$$

(5)

where $F_{\mathrm{ij}}$ is the force exerted over nanoparticle i due to van der Waals’ interaction energy between the nanoparticle i and nanoparticle j, ⁵ that is

$$F_{\mathrm{ij}}=\frac{{\mathrm{AR}}_i^3}{D_{\mathrm{ij}}^4}$$

where $D_{\mathrm{ij}}$ is the distance between the nanoparticles i and j, and $F_{\mathrm{iw}}$ is the force exerted over nanoparticle i due to van der Waals’ interaction energy between the nanoparticle i and the walls, that is

$$F_{\mathrm{iw}}=|F_{i{w}_1}-F_{i{w}_2}|$$

where $F_{i{w}_k}$ is the force exerted over nanoparticle i due to van der Waals’ interaction energy between the nanoparticle i and the wall k $(k=1,2)$ and given by ⁵

$$F_{i{w}_k}=\frac{{\mathrm{AR}}_i^3}{{\mathrm{dw}}_{k,i}^4}$$

where $d{w}_{k,i}$ is the distance between the nanoparticle i and the wall k $(k=1,2)$ (see Figure 4), that is, $d{w}_{2,i}=a-z_i-R_i$ and $d{w}_{1,i}=a+z_i-R_i$ . Here, $v_{x_i}=d{x}_i/\mathrm{dt}$ and $v_{z_i}=d{z}_i/\mathrm{dt}$ . Clearly, $F_{\mathrm{iw}}=0$ , whenever the nanoparticle i is on the axis of the tube, since $F_{i{w}_1}=F_{i{w}_2}$ . The function $h_5^{w_k}$ is called the hemodynamic resistant force^5,6 induced by the wall $w_k$ , $(k=1,2)$ , on the z-axis which is given by ⁴

$$h_5^{w_k}(z_i)=6\pi \mu R_i[-\frac{1}{1-\frac{R_i}{z_i}}+\frac{1}{8}\log (1-\frac{R_i}{z_i})]$$

(6)

Figure 4.

Motion of 500 nanoparticles inside a straight microtube.

On x-axis, we have, for any nanoparticle i, the following equation

$$[h_1^{w_1}(z_i)+h_1^{w_2}(z_i)]·v_{x_i}=-[h_6^{w_1}(z_i)+h_6^{w_2}(z_i)]·S$$

(7)

where S is the shear rate of the flow and the functions $h_1^{w_k}$ and $h_6^{w_k}$ are also called the hemodynamic resistant force induced by the wall $w_k$ , $(k=1,2)$ , on x-axis which is given by ⁴

$$h_1^{w_k}(z_i)=\frac{48}{15}\pi \mu R_i[\ln (1-\frac{R_i}{R_i+d{w}_{k,i}})+\frac{R_i}{R_i+d{w}_{k,i}}]$$

(8)

and

$$\begin{array}{l}{h}_6^{w_k}(z_i)=6\pi \mu R_i\\ \left[1+\frac{9}{16}\left(\frac{R_i}{R_i+d{w}_{k,i}}\right)+0.13868{\left(\frac{R_i}{R_i+d{w}_{k,i}}\right)}^{1.4829}\right]\end{array}$$

(9)

Equations (5) and (7) are, respectively, the force balance in the x- and z-directions. In these equations, $R_i$ is the radius of nanoparticle i, ${\Delta }_{\rho }$ is the difference between mass densities of the nanoparticles and the fluid, and g is the gravitational acceleration. Combining equations (5) and (7), we get the following differential system driving the nanoparticles

$$\{\begin{array}{c}\frac{d{x}_i}{\mathrm{dt}}=\frac{-[h_6^{w_1}(z_i)+h_6^{w_2}(z_i)]·S}{h_1^{w_1}(z_i)+h_1^{w_2}(z_i)}\hfill \\ \frac{d{z}_i}{\mathrm{dt}}=\frac{\frac{4\pi }{3}\mathrm{gR}\Delta \rho +\sum _{j\ne i}{F}_{\mathrm{ij}}+F_{\mathrm{iw}}}{h_5^{w_1}(z_i)+h_5^{w_2}(z_i)}\hfill \end{array}$$

(10)

Remark 2

1. For a spherical particle, the hydrodynamic resistance function $h_i$ has been given in the literature.^7–9 They are given either in the form of tabulated results of numerical computation or in the form of asymptotic expressions for a sphere very close to or very far away from the bounding wall.

2. The approximations of the hemodynamic resistant function $h_i^{w_k}$ are derived using a method proposed by Zhao et al.^4,5

3. The two equations for the force balance contain two unknowns, namely, the two velocity components $v_x$ and $v_z$ .

4. The function $h_i$ does not depend on the translational velocity $v_{x_i}$ along x axis.

Illustrative simulation

In our simulations, we assume that $N=500$ , $l=7\mu \mathrm{m}$ , and $h=0.0075$ ; the radius of a nanoparticle i is taken to be a random value in the nanointerval $[5\mathrm{nm};10\mathrm{nm}]$ , that is, the minimum possible value of radius is $R_i=5\mathrm{nm}$ and the maximum possible value of $R_i=10\mathrm{nm}$ ; $T=10$ ; $A=5e-21$ ; $g=9.81e9$ ; $\mu =3.4e-3$ ; $\Delta \rho =1e3$ ; and $a=0.5\mathrm{\mu }\mathrm{m}$ . We start with an initial configuration $q_0$ at time $t=0$ . The initial configuration is a random distribution of all the nanoparticles. In order to plot the next configuration $q_1$ at time $t=1$ , we solve the differential system presented in section “Dynamical system of the motion of nanoparticles” to get the desired velocity of all the nanoparticles by getting the values of $v_i=(v_{x_i},v_{y_i})$ and so the value of $U(q_0)$ . In order to avoid the overlapping of the disks (representing the nanoparticles) in the next configuration $q_1$ , we project the obtained value of $U(q_0)$ on the admissible velocity set $C_h(q_0)$ . Once the projection is done, we get the vector $w^1$ , ensuring the non-overlapping, we plot the next configuration via the simple formula $q_1=q_0+h{w}^1$ and so on until the final time $t=10$ . Our numerical results are shown in Figure 4.

Discussion

The first image in Figure 4 shows the initial situation (taken randomly) of the nanoparticles inside the microtube with diameter $a=0.5\mathrm{\mu }\mathrm{m}$ . During the time values $t\in \{1,\dots ,10\}$ , the motion of the nanoparticles appears exactly as we expected, that is, the more the particles are far from the walls, the more they are faster, in other words, the speed is directly proportional to the distance to walls.

Conclusion

In this work, we considered the simple situation (straight tube and only van der Waals’ interaction) in order to setup an approach to produce a realistic numerical simulation for a huge number of moving nanoparticles. Our next task which seems to be doable, due to the flexibility of our iterative scheme, is to insert all the possible forces acting on the nanoparticles, some existing obstacles (fixed or moving), and also the general form of the function defining the shape of the tube.