Direct numerical simulation of noninvasive channel healing in electrical field

Abstract

Noninvasive channel healing is a new idea to repair the broken pipe wall, using external electric fields to drive iron particles to the destination. The repair can be done in the normal operation of the pipe flow without any shutdown of the pipeline so that this method can be a potentially efficient and safe technology of pipe healing. However, the real application needs full knowledge of healing details. Numerical simulation is an effective method. Thus, in this research, we first established a numerical model for noninvasive channel healing technology to represent fluid–particle interaction. The iron particles can be attached to a cracking area by external electrostatic forces or can also be detached by mechanical forces from the fluid. When enough particles are permanently attached on the cracking area, the pipe wall can be healed. The numerical criterion of the permanent attachment is discussed. A fully three-dimensional finite difference framework of direct numerical simulation is established and applied to different cases to simulate the full process of channel healing. The impact of Reynolds number and particle concentration on the healing process is discussed. This numerical investigation provides valuable reference and tools for further simulation of real pipe healing in engineering.

Introduction

Piping systems are very important in petroleum transportation engineering. After long-period operations, the pipelines usually face the problem of corrosion, which causes cracks on the pipe wall and increases pipelining risks. ¹ Standard treatment of this problem in engineering currently is repair or replacement of the damaged segments, which involves high capital cost. ² This requires shutdown of the whole transportation system, including the pipeline and affiliated pumps, heat exchangers and power systems. Interruption of petroleum transportation has great impact on downstream oil refinery factories and users, resulting in large economic loss usually. For pipelines transporting waxy oil, whose pour point is higher than regular oil, the oil gels because of lower temperature during the period of shutdown. Gelled oil blocks the pipeline so that restart fails to bring second hazard to the pipeline. Expensive treatment, economic loss and second hazard motivate to derive a strategy of pipeline repair without shutdown. Zohdi proposed a concept of “noninvasive” pipe healing using charged iron particles guided by electromagnetic field from the exterior of the pipeline, which can heal the pipe wall resulting in no shutdown so that the cost may be largely reduced. However, this idea is only in a primary stage with so many issues not being addressed, including particle attachment and detachment, particle–fluid–wall interaction in an electronic field, and so on. Some results on biological systems can be a good reference.^3–5 Also, some studies have been carried out by Zohdi and colleagues,^6,7 focusing on particle adhesion to walls, but the simple notion of a detachment stress threshold may be inadequate. Recently, Mukherjee et al. ⁸ provided a numerical framework for the idea of the noninvasive pipe healing technique. They treat the particle–wall interaction very well in laminar flow. The flow status in real pipelines for petroleum engineering must be turbulence. Therefore, we propose a numerical model for channel healing under the condition of three-dimensional turbulent flow and do the direct numerical simulation (DNS) in this paper.

Numerical model

Fluid flow model

The pipeline transporting petroleum is usually as long as hundreds to thousands of kilometers. A simplification in engineering is done using cross-sectional average velocity instead of local velocity so that calculation can be done for this length of pipeline. This assumption is not suitable for the simulation of pipe healing because the spatial distribution of flow field, which should be considered, is neglected. Thus, a fully three-dimensional (3D) flow model with corresponding governing equations should be established. However, the 3D model requires a dense mesh to resolve small eddies, so that the total grid number for a several-kilometer-long pipeline can be too large to exceed the storage of computers (even if super computers). Hence, flow domain must be simplified first. Limited shorter domain with periodic boundary condition is a good choice, as shown in Figure 1. The boxes with symbol “+” and “–” represent the positive and negative electrodes installed outside of channel walls. They are located on the center of x–z plane. The computational domain has side lengths of 7.5h, 2h, 2.5h in the x, y, z directions, respectively (h is the half height of the channel). The crack is assumed at the central area with the dimension of 1.5h × 0.04h × 0.5h (x × y × z) on the bottom wall. Periodic flow boundary condition is imposed on the x and z directions.

Figure 1.

Illustration of channel flow with external electric field.

In this computational domain, the governing equations describing the three-dimensional, unsteady incompressible viscous fluid flow are continuity equation and Navier–Stokes (N-S) equation in tensor forms:

• Continuity equation

$$\frac{\partial \mathbf{u}}{\partial \mathbf{x}}=0$$

(1)

• N-S equation

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\nabla )\mathbf{u}=-\frac{1}{{\rho }_f}\nabla p+\frac{{\mu }_f}{{\rho }_f}\Delta \mathbf{u}+\mathbf{f}$$

(2)

where u is the fluid velocity vector with three components u, v, and w in x, y, and z directions. $\nabla p$ is the pressure gradient, ρ_f is the density of fluid, μ_f is the dynamic viscosity of fluid, and t is time. x is the Cartesian axis direction vector representing the three directions (x, y, z). f is the force per unit mass imposed on fluid. The particles used for pipe healing are tiny powder so that gravity can be neglected and particles can be assumed spherical. Moreover, it is apparent that the particles can be accumulated to the bottom wall more easily with gravity. The addition of gravity may cover the effect of electrostatic force on the accumulation, which is the main target of this study. Thus, neglect of gravity is necessary in this study. After this discussion, f only contains fluid–particle interaction force. According to Newton’s third law, we can derive expression of f as follows:

$$\mathbf{f}=\frac{\sum _{\frac{i=1}{i\in V_{\mathbf{u}}}}^N{\mathbf{F}}_{\mathrm{fi}}^{\mathrm{drag}}}{{\rho }_f{V}_{\mathbf{u}}}$$

(3)

where ${\mathbf{F}}_{\mathrm{fi}}^{\mathrm{drag}}$ is the drag force of the fluid imposed by the ith particle in control volume $V_{\mathbf{u}}$ . Thus, all the N particles in $V_{\mathbf{u}}$ have an overall drag force $\sum _{\frac{i=1}{i\in V_{\mathbf{u}}}}^N{\mathbf{F}}_{\mathrm{fi}}^{\mathrm{drag}}$ to the corresponding fluid. ${\rho }_f{V}_{\mathbf{u}}$ is the mass of the fluid in the control volume. According to Newton’s third law of motion, the ith particle is also driven by a drag force ${\mathbf{F}}_{\mathrm{pi}}^{\mathrm{drag}}$ imposed by the fluid in $V_{\mathbf{u}}$ . These relations are illustrated in Figure 2. The cubic box represents a control volume of fluid, that is, an element of mesh $V_{\mathbf{u}}=\Delta x_i\Delta y_j\Delta z_k$ . Circles represent particles in the volume. ${\mathbf{F}}_{\mathrm{fi}}^{\mathrm{drag}}$ and ${\mathbf{F}}_{\mathrm{pi}}^{\mathrm{drag}}$ are the forces imposed on the fluid and particle, respectively, with the following relation:

$${\mathbf{F}}_{\mathrm{fi}}^{\mathrm{drag}}=-{\mathbf{F}}_{\mathrm{pi}}^{\mathrm{drag}}$$

(4)

Figure 2.

Illustration of fluid–particle interaction in one control volume.

Particle motion model

The above discussions show it is important to deal with the fluid–particle interaction force ${\mathbf{F}}_{\mathrm{pi}}^{\mathrm{drag}}$ on each particle. If the impact of particle shape on flow field is considered, small-scale irregular mesh around every particle should be used so that the mesh density is very high to bring huge computational load. Actually, shape influence of particles is usually low when the size of the particle is small. The size for pipe healing is indeed small for tiny iron powder particles. Thus, the following assumptions are made: (1) a particle only has influence on the nearby flow field and does not affect the fluid flow far from itself, (2) rotation of the particle is neglected, and (3) the particles in fluid flow mesh are considered as moving points with mass but without volume. The following relations can be obtained:

• Particle motion equation

$$m_p\frac{d{\mathbf{v}}_{\mathrm{pi}}}{\mathrm{dt}}={\mathbf{F}}_{\mathrm{pi}}^{\mathrm{drag}}+{\mathbf{F}}_{\mathrm{pi}}^{\mathrm{elec}}$$

(5)

• Drag force

$${\mathbf{F}}_{\mathrm{pi}}^{\mathrm{drag}}=\frac{1}{2}{\rho }_f{C}_D\pi r_p^2({\mathbf{u}}_{\mathrm{pi}}-{\mathbf{v}}_{\mathrm{pi}})\Vert {\mathbf{u}}_{\mathrm{pi}}-{\mathbf{v}}_{\mathrm{pi}}\Vert$$

(6)

where $m_p$ , ${\mathbf{v}}_{\mathrm{pi}}$ , and $r_p$ are the mass, velocity, and radius, respectively, of the ith particle. All particles are assumed to have same mass and radius. ${\mathbf{F}}_{\mathrm{pi}}^{\mathrm{elec}}$ is the force imposed on the ith particle by an external electric field. ${\mathbf{v}}_{\mathrm{pi}}$ is the velocity of the ith particle. ${\mathbf{u}}_{\mathrm{pi}}$ is the fluid velocity at the position of the ith particle. “‖·‖” means L2 norm. $C_D$ is the drag coefficient. To calculate $C_D$ , local Reynolds number based on particle radius should be defined as follows

$$R{e}_{\mathrm{pi}}=\frac{2r_p{\rho }_f\Vert {\mathbf{u}}_{\mathrm{pi}}-{\mathbf{v}}_{\mathrm{pi}}\Vert }{{\mu }_f}$$

(7)

• Stokes’ drag law can be used for the calculation of $C_D$

$$C_D=\frac{24}{R{e}_{\mathrm{pi}}}$$

(8)

Stokes’ law only satisfies $R{e}_{\mathrm{pi}}\in [0,1]$ , which is very low speed laminar flow. In real pipe flow, the flow status must be turbulent even in the low-speed viscous sublayer. Thus, a correlation should be used for wider range of Reynolds number. There are numerous reasonable corrections to the drag coefficient C_D. Haider and Levenspiel ⁹ find an accurate estimation model compared to the available experimental data for the pipe-healing system, in which the coefficient C_D are defined as

$$C_D=\frac{24}{R{e}_{\mathrm{pi}}}(1+0.1806{\mathrm{Re}}_{\mathrm{pi}}^{0.6459})+\frac{0.4251R{e}_{\mathrm{pi}}}{R{e}_{\mathrm{pi}}+6880.95}$$

(9)

The electric force imposed on the ith particle can be calculated from the definition of electrostatic force

$${\mathbf{F}}_{\mathrm{pi}}^{\mathrm{elec}}=\frac{k_e\mathrm{Qq}}{r_i^2}\frac{{\mathbf{r}}_i}{r_i}$$

(10)

where $k_e$ is the electrostatic coefficient, and Q and q are the charges carried by the wall and each particle, respectively. ${\mathbf{r}}_i$ is the vector from the center (x_c, y_c, z_c) of the bottom wall to the particle, ${\mathbf{r}}_i=(x_c-x_{\mathrm{pi}})\mathbf{i}+(y_c-y_{\mathrm{pi}})\mathbf{j}+(z_c-z_{\mathrm{pi}})\mathbf{k}$ . $r_i$ is the length of the vector ${\mathbf{r}}_i:r_i=\Vert {\mathbf{r}}_i\Vert$ .

It is not necessary to calculate ${\mathbf{v}}_{\mathrm{pi}}$ using equation (5) directly because the mass of every particle is same, that is, $m_p={\rho }_p(4/3)\pi r_p^3$ . Therefore, mass of the particle in equation (5) should be divided to obtain the particle motion equation per mass

$$\frac{d{\mathbf{v}}_{\mathrm{pi}}}{\mathrm{dt}}={\mathbf{f}}_{\mathrm{pi}}^{\mathrm{drag}}+{\mathbf{f}}_{\mathrm{pi}}^{\mathrm{elec}}$$

(11)

where ${\mathbf{f}}_{\mathrm{pi}}^{\mathrm{drag}}$ and ${\mathbf{f}}_{\mathrm{pi}}^{\mathrm{elec}}$ are drag force per unit mass and electrostatic force per unit mass, respectively. They have following expressions

$${\mathbf{f}}_{\mathrm{pi}}^{\mathrm{drag}}=\frac{{\mathbf{F}}_{\mathrm{pi}}^{\mathrm{drag}}}{m_p}=\frac{3}{8r_p{d}_p}{C}_D({\mathbf{u}}_{\mathrm{pi}}-{\mathbf{v}}_{\mathrm{pi}})\Vert {\mathbf{u}}_{\mathrm{pi}}-{\mathbf{v}}_{\mathrm{pi}}\Vert$$

(12)

$${\mathbf{f}}_{\mathrm{pi}}^{\mathrm{elec}}=\frac{{\mathbf{F}}_{\mathrm{pi}}^{\mathrm{elec}}}{m_p}=\frac{k_e\mathrm{Qq}}{r_i^3}\frac{3}{4\pi r_p^3{\rho }_p}{\mathbf{r}}_i$$

(13)

where ${\rho }_p$ is the density of iron particle, and $d_p$ is the density ratio of iron to fluid ( $d_p={\rho }_p/{\rho }_f$ ).

Particle–wall interaction model

Particle-wall interaction includes two parts: attachment on the wall driven by electrostatic force and collision with wall. For the first part, an energy criterion is proposed for the accumulation of particles: a particle can be accumulated on the target wall if the electrostatic energy is larger than the momentum energy of the particle. In this article, only the y direction is the wall-normal direction so that only the component of velocity in this direction is considered in the momentum energy

$$E_{\mathrm{vi}}=\frac{1}{2}{m}_p{v}_{\mathrm{pi}}^2$$

(14)

where $E_{\mathrm{vi}}$ is the momentum energy of the ith particle, and $v_{\mathrm{pi}}$ is the velocity component of the ith particle in y direction. The expression of electrostatic energy can be obtained from the following integration

$$E_{\mathrm{pi}}={\int }_{r_i}^{+\infty }\Vert {\mathbf{F}}_{\mathrm{pi}}^{\mathrm{elec}}\Vert d{r}_i={\int }_{r_i}^{+\infty }\frac{k_e\mathrm{Qq}}{r_i^3}\Vert {\mathbf{r}}_i\Vert d{r}_i={\int }_{r_i}^{+\infty }\frac{k_e\mathrm{Qq}}{r_i^2}d{r}_i=k_e\mathrm{Qq}{\int }_{r_i}^{+\infty }\frac{1}{r_i^2}d{r}_i=\frac{k_e\mathrm{Qq}}{r_i}$$

(15)

If $E_{\mathrm{vi}}\le E_{\mathrm{pi}}$ , the particle should be attached on the wall. If $E_{\mathrm{vi}}>E_{\mathrm{pi}}$ , the particle should be detached from the wall.

For the second part, elastic collision is the simplest assumption. However, it does not reflect the real collision which contains energy lost. Elastic collision without energy loss may also cause the particles take periodic motion between the target wall and the bulk flow, so that no particles can be accumulated on the wall. Thus, nonelastic collision is considered

$$v_{\mathrm{pi}}=-C_R{v}_{\mathrm{pi}}$$

(16)

where $C_R$ is the collision coefficient with the value 0.597 for iron.

Dimensionless governing equations

All the governing equations are nondimensionalized using the following dimensionless variables

$${\mathbf{u}}^{+}=\frac{\mathbf{u}}{u_{\tau }},{\mathbf{v}}_p^{+}=\frac{{\mathbf{v}}_p}{u_{\tau }},{\mathbf{x}}^{*}=\frac{\mathbf{x}}{h},r_p^{*}=\frac{r_p}{h},t^{*}=\frac{t}{h/u_{\tau }},p^{+}=\frac{p}{{\rho }_f{u}_{\tau }^2}=\overline{p^{+}}+p^{+*},R{e}_{\tau }=\frac{{\rho }_f{u}_{\tau }h}{{\mu }_f}$$

(17)

$$C=\frac{k_e\mathrm{Qq}}{{\rho }_p{u}_{\tau }^2{h}^4}$$

(18)

where the superscript (⁺) and (^*) denote the dimensionless variables; ${\mathbf{u}}^{+}$ and ${\mathbf{v}}_p^{+}$ represent the dimensionless velocity vector of fluid and particles respectively; $t^{*}$ and $p^{+*}$ represent the dimensionless pressure and the dimensionless time, respectively; $r_p^{*}$ represents the dimensionless radius; $x^{*}$ represents the dimensionless Cartesian coordinates; and C is a new dimensionless variable representing the ratio between electrostatic force and inertial force.

• Dimensionless flow equations

$$\frac{\partial {\mathbf{u}}^{+}}{\partial {\mathbf{x}}^{*}}=0$$

(19)

$$\frac{\partial {\mathbf{u}}^{+}}{\partial t}+({\mathbf{u}}^{+}·\nabla ){\mathbf{u}}^{+}=\mathbf{e}-\nabla p^{+*}+\frac{1}{R{e}_{\tau }}\Delta {\mathbf{u}}^{+}+{\mathbf{f}}^{+}$$

(20)

where $\mathbf{e}={[1,0,0]}^T$ , only x direction has value.

• Dimensionless particle motion equations

$$\frac{d{\mathbf{v}}_{\mathrm{pi}}^{+}}{d{t}^{*}}={\mathbf{f}}_{\mathrm{pi}}^{+\mathrm{drag}}+{{\mathbf{f}}^{+}}_{\mathrm{pi}}^{\mathrm{elec}}$$

(21)

$${\mathbf{f}}_{\mathrm{pi}}^{+\mathrm{drag}}=\frac{3}{8r_p^{*}{d}_p}{C}_D({\mathbf{u}}_{\mathrm{pi}}^{+}-{\mathbf{v}}_{\mathrm{pi}}^{+})\Vert {\mathbf{u}}_{\mathrm{pi}}^{+}-{\mathbf{v}}_{\mathrm{pi}}^{+}\Vert$$

(22)

$${{\mathbf{f}}^{+}}_{\mathrm{pi}}^{\mathrm{elec}}=\frac{C}{V_p^{*}{r}_i^{*3}}{\mathbf{r}}_i^{*}$$

(23)

• Dimensionless electrostatic energy on the particle positions

$$E_{\mathrm{pi}}^{+}=\frac{C}{r_i^{*}}$$

(24)

• The kinetic energy of particles after collision with walls

$$E_{\mathrm{vi}}^{+}=\frac{1}{2}{V}_p^{*}{v}_{\mathrm{pi}}^{+2}$$

(25)

where

$$V_p^{*}=\frac{4}{3}\pi r_p^{*3}$$

Numerical results of channel healing

The numerical model is solved by the fractional step method with grid number 64 × 64 × 64 and dimensionless time step $\Delta t^{*}={10}^{-4}$ . The main numerical method is DNS. Reynolds numbers Re_τ = 10 and Re_τ = 150 are selected. The dimensionless radius of each iron particle is 0.01. Total number of iron particles added into the flow is set to 1000 and 10000. The dimensionless size of the crack is 1.5 × 0.04 × 0.5. The initial condition is set to be fully developed laminar flow for the case Re_τ = 10 and turbulent flow for the case Re_τ = 150. The boundary condition is the periodic condition in the x and z directions, and the no-slip boundary condition in the y direction. When the total volume of the particles accumulated in the crack equals to the volume of the crack, the channel is considered to be healed so that the computation is completed. It should be noted that the dimensionless depth of the crack is very tiny (0.04) compared with the dimensionless height of the channel (2) because of which we assume the flow field is not affected by the crack. This assumption makes the computation of the flow field much more easily but needs to be improved in future.

First, we should check the flow status. The mean velocity profiles for the case Re_τ = 10 and Re_τ = 150 are shown in Figure 3. It is obvious that the mean velocity profile for the case Re_τ = 10 is far below the turbulent three layers, ¹⁰ while the mean velocity profile for the case Re_τ = 150 agrees very well with the typical three-layer-mean-velocity profile of turbulent flow. Moreover, the bulk mean velocities (U_m⁺) and the bulk mean Reynolds number (Re_m) are listed in Table 1. Re_m for the cases Re_τ = 10 and Re_τ = 150 are 67 and 4443, respectively. Thus, the two designed cases represent laminar flow and turbulent flow in the real pipeline.

Figure 3.

Mean velocity profile of fluid flow under different friction Reynolds number.

Table 1.

Bulk mean velocity and Reynolds number.

	Um+	Rem
Reτ = 10	3.3	67
Reτ = 150	14.8	4443

Second, it is important to study the motion of particles and their influences on the flow field in the healing process, which can clarify the influence of the noninvasive healing method on the transporting fluid in pipeline. From the left columns in Figures 4–7, we can see that the flow field does not change much under different numbers of particles and different Reynolds numbers. This indicates that the noninvasive healing method is not sensitive to the inputting particle numbers and does not apparently affect the flow field, because of which the fluid transportation is not affected by this healing method, which is the main advantage of this method that has been verified by numerical results.

Figure 4.

Flow field (left) and particle positions (right) at different moments of channel healing in laminar flow (1000 particles). a) 0 $\Delta t^{*}$ : 40 particles attached, (b) 200000 $\Delta t^{*}$ : 378 particles attached, (c) 400000 $\Delta t^{*}$ : 544 particles attached, (d) 1000000 $\Delta t^{*}$ : 881 particles attached, (e) 1470000 $\Delta t^{*}$ : 1000 particles attached (unhealed).

Figure 5.

Flow field (left) and particle positions (right) at different moments of channel healing in laminar flow (10000 particles). (a) 0 $\Delta t^{*}$ : 400 particles attached, (b) 150000 $\Delta t^{*}$ : 3563 particles attached, (c) 350000 $\Delta t^{*}$ : 5636 particles attached, (d) 450000 $\Delta t^{*}$ : 6705 particles attached, (e) 495000 $\Delta t^{*}$ : 7166 particles attached (healed).

Figure 6.

Flow field (left) and particle positions (right) at different. (a) 0 $\Delta t^{*}$ : 40 particles attached, (b) 20000 $\Delta t^{*}$ : 550 particles attached, (c) 60000 $\Delta t^{*}$ : 879 particles attached, (d) 80000 $\Delta t^{*}$ : 934 particles attached, (e) 100000 $\Delta t^{*}$ : 967 particles attached, (f) 180000 $\Delta t^{*}$ : 1000 particles attached (unhealed).

Figure 7.

Flow field (left) and particle positions (right) at different moments of channel healing in turbulent flow (1000 particles) moments of channel healing in turbulent flow (10000 particles). (a) 0 $\Delta t^{*}$ : 400 particles attached, (b) 4000 $\Delta t^{*}$ : 2346 particles attached, (c) 12000 $\Delta t^{*}$ : 4198 particles attached, (d) 16000 $\Delta t^{*}$ : 4863 particles attached, (e) 20000 $\Delta t^{*}$ : 5572 particles attached, (f) 36000 $\Delta t^{*}$ : 7166 particles attached (healed).

Comparing Figure 4 with Figure 6, we find that the pipe wall cannot be healed for laminar flow and turbulent flow inputting 1000 particles because the number of particles is not enough to fill the crack. This can be expected according to the volume of particles and the volume of the crack. However, the accumulation of the 1000 particles takes a much longer time than the accumulation of 7166 particles for inputting 10000 particles, for both the laminar flow and turbulent flow (Figures 5 and 7). This phenomenon shows that higher concentration of particles is helpful for channel healing. The reason is that more frequent collisions in higher concentration make the particles more easily overcome the resistance of fluid to reach the wall. Thus, high concentration of particles should be used in the real noninvasive healing process.

It can be seen in Figures 5 and 7 that the final healing time for the same inputting 10000 particles is much shorter in turbulent flow than in laminar flow. The number of accumulated particles in the whole process is shown in Figure 8 (N_P is the number of accumulated particles). It is further verified that the particles are accumulated much faster in turbulent flow field than in laminar flow field at any moment. Their total dimensionless time is about 3.4 and 49.4, respectively. Turbulent flow can accelerate the healing process about 15 times. This is probably because the strong convection in turbulence makes the particles have larger velocity to overcome the flow resistance in the wall-normal direction more easily. Therefore, a high Reynolds number should be used for faster healing.

Figure 8.

Healing time under different friction Reynolds numbers when inputting 10000 particles.

Conclusion

A numerical model is established for the new concept of channel healing via adding small iron particles driven by an external electric field. Fluid–particle interaction in electric field and particle–wall interaction are discussed in this model. A criterion is proposed for numerically judging channel healing according to the balance of electrostatic energy and momentum energy. The direct numerical simulations are made under different numbers of particles (1000 and 10000) and different flow status (Re_τ = 10 for laminar flow and Re_τ = 150 for turbulent flow). The numerical results show that the fluid transportation process is not affected in the healing process, verifying the noninvasive advantage of this kind of healing method. Key parameters of the noninvasive healing are studied. Higher concentration of particles and higher Reynolds number are helpful for faster channel healing because the two factors make the particles have more frequent collisions and higher convection due to which they have much larger velocity to overcome the wall-normal resistance from the fluid.