The Development of a Wax Layer on the Interior Wall of a Circular Pipe Transporting Heated Oil -- The Effects of Temperature Dependent Wax Conductivity

In this paper we develop and significantly extend the thermal phase change model, introduced in [12], describing the process of paraffinic wax layer formation on the interior wall of a circular pipe transporting heated oil, when subject to external cooling. In particular we allow for the natural dependence of the solidifying paraffinic wax conductivity on local temperature. We are able to develop a complete theory, and provide efficient numerical computations, for this extended model. Comparison with recent experimental observations is made, and this, together with recent reviews of the physical mechanisms associated with wax layer formation, provide significant support for the thermal model considered here.


Introduction
The transport of oil in long subsea pipelines occurs in many oil-producing fields. In a recent paper [12], one of the key challenges arising in subsea development of oil fields, namely the formation of paraffinic wax deposits on the inside of the pipe wall, was considered. This wax deposit happens when the temperature of the pipe wall falls below the wax solidification temperature, generally in the range 35 − 40 • C. Deposition of wax on the pipe wall has significant operational consequences since the inner diameter of the pipe will decrease thereby reducing the transport capacity of the pipe for a given driving pressure. A more detailed discussion of this phenomena is given in [12]. In particular [12] was concerned with the introduction of a thermal phase change model (first proposed by Schulkes in [14]) to capture the fundamental mechanism in the wax layer deposition process accurately. The associated mathematical model was both fundamental and analysed in detail in [12]. The outcomes were encouraging, in that a number of key observations in the wax layer deposition process, which were at odds with the historic material diffusion mechanism, were now fully accounted for in the thermal phase change model. This model has had further recent support in the independent works in [1] and [10] 1 .
In [12], a thermal approach towards modelling the deposition of paraffinic wax on the inside of the pipe wall was proposed. The basic assumption introduced in [12] is that the principal mechanism leading to wax deposition is a thermal phase change process. This led to a simple model based on the fundamental balances of heat conduction from the heated oil, through the evolving solid wax layer and pipe wall, and finally into the cooling fluid surrounding the pipe. The mathematical formulation of this model gave rise to a free boundary problem (referred to as [IBVP] in [12]) of generalised Stefan type. This fundamental problem was analysed in considerable detail in [12]. A number of key salient features arising from the model were identified (see [12, section 8, p.119-122]) with the intention of qualifying and quantifying the basis of the thermal model with detailed experiments performed by Hoffman & Amundsen [7] and Halstensen et al [6] amongst others. The experiments of Hoffman & Amundsen [7] show that with a constant flow rate, the wax layer reaches an equilibrium height after sufficient time with this time decreasing as the oil temperature increases. This basic feature is not predicted by the molecular diffusion model that is widely applied in the oil and gas industry [15] [16] and a "shear-stripping" mechanism has to be introduced to match experimental observation with model predictions. A modelling approach based on the thermal phase change mechanism as outlined in [12] shows that some fundamental experimental observations can be explained without the need to include rather speculative physical mechanisms. Very recently, the extensive review by Mehrotra et al [11] has given a thorough consideration of experimental evidence, which provides significant and substantial support for the thermal phase change mechanism introduced in [12], as the principal and key mechanism in the process of wax deposition on the interior wall of pipes transporting heated-oil. Most recently, a further review has been provided by Van der Geest et al [17], which again supplies detailed and critical experimental support for the thermal phase change mechanism. In addition to these reviews, a recent thesis has addressed this issue experimentally, Mahir [10], and has also provided significant support for the thermal mechanism proposed in [12] 2 . With the very recent emergence of this significant and compelling support for the basic thermal model introduced in [12], the purpose of the present paper is to investigate how the mathematical model in [12] can be further developed to accommodate more detailed features associated with wax phase change. Specifically, in this paper, we develop the model in [12] to account for the dependence of the thermal conductivity of the paraffinic wax solid on local temperature. This effect can be significant in many paraffinic wax materials when in solid state and is therefore an effect which should be investigated in the modelling process. With this inclusion in the mathematical model developed in [12], the significant change is that the associated free boundary problem [IBVP] becomes nonlinear. In particular, the partial differential equation, rather than being the associated linear heat conduction equation, becomes a nonlinear (in fact uniformly quasi-linear) strictly parabolic equation, with, in addition, the mixed Robin-type linear boundary condition on the solid interior pipe wall, and the latent heat boundary condition on the free solid wax layer surface, both adopting a generalised nonlinear form. The principal aim of the current paper is to fully investigate the effects of the inclusion of the temperature dependent solid wax thermal conductivity in the model developed in [12]. In Section 2, we review and extend the thermal phase change model and it's mathematical formulation. This is followed in Section 3 by an extensive study of associated qualitative results for [IBVP]. These structural results are then complemented by consideration of the nature of the solution to [IBVP] as t → 0 + and t → ∞. In Section 5 we devote attention to developing a complete and tractable theory for [IBVP] when the parameter ε is small, a case which often pertains in physical applications. In Section 6 we consider numerical solutions to [IBVP] for comparison with the theory developed earlier, whilst Section 7 gives a qualitative comparison with experiments presented in [7]. Finally we end with a discussion in Section 8.

The Model
In this section, following [12], we formulate the thermal phase change model which was discussed in section 1. The heated oil is in uniform flow through a long, straight section of pipe, with circular cross section, and internal radius R. We first observe that when the wax layer thickness on the interior pipe wall h is very much smaller than the pipe radius R (h R), which is often the case in applications, then we can reduce the problem to a planar geometry. Following [12] we restrict attention to a long section of the pipe which is remote from entry and exit effects. Under these circumstances, we may consider all properties in the model to be independent of axial distance along the pipe. Thus all dependent variables in the model are functions of x and t alone, where x measures normal distance from the inner pipe wall towards the pipe axis, and t is time. The situation is illustrated in Figure 1. The pipe is surrounded by an aligned circular coolant jacket. The fluid in the coolant jacket is maintained at constant temperature T c , with the coolant jacket having width d c . The thickness of the pipe wall is d p , with the outer pipe wall located at x = −d p and the inner pipe wall located at x = 0. The temperature within the pipe wall is denoted by T p . The solid wax layer is initiated, when t = 0, at the inner pipe wall x = 0, with its upper surface at x = h. The temperature within the solid wax layer is denoted by T , with T h being the constant temperature of the wax formation, at x = h. The temperature of the oil flowing in the pipe, when h < x ≤ R, is taken as constant, and represented by T o , an approximation which observations confirm to be reasonable in applications [1]. Throughout, we consider the situation when, We now consider the temperature field in the pipe wall. An application of Fourier's law gives where k p , c p and ρ p are the conductivity, specific heat capacity and density of the pipe wall respectively. With t s (as given in (2.12) and (2.18)) being the time scale associated with wax layer formation, the thickness of the pipe wall is thin, and such that  Consequently, the temperature in the pipe wall is in a quasi-steady state, and equation (2.2) may be approximated by subject to the following boundary conditions, Here k c is the coolant conductivity, Nu c is the Nusselt number for the coolant flow and k w (T ) is the temperature dependent solid wax conductivity. Condition (2.4) represents continuity of heat flux across the exterior wall of the pipe, which is in contact with the coolant jacket, whilst condition (2.5) represents continuity of heat flux from the lower boundary of the wax layer into the interior pipe wall, and condition (2.6) represents continuity of temperature at the lower boundary of the wax layer and the inner pipe wall. In allowing the solid wax conductivity to be temperature dependent, we write for −∞ < T < ∞, where k h w is the conductivity of the solid wax at the temperature of wax solidification T = T h , and D : R → R is the dimensionless solid wax conductivity, representing the variation in solid wax conductivity with, In general D : R → R will be taken as smooth, bounded and bounded above zero. Thus, throughout, we will consider D : R → R to satisfy the conditions: after which boundary condition (2.5) then requires, which is a condition on T at x = 0, the lower boundary of the solid wax layer at the inner pipe wall. We are now in a position to consider the temperature field T in the solid wax layer. Fourier's law requires that, subject to the boundary conditions, where ρ w , c w and H w are the density, specific heat capacity and latent heat of the solid wax respectively, whilst k o and Nu are the conductivity of oil and the Nusselt number for the oil flow respectively. The conditions (2.9) and (2.10) express that the outer surface of the solid wax layer must be at the wax solidification temperature, and that the difference in heat flux across this interface balances the latent heat required for solid wax formation. For convenience we nondimensionalise the free boundary problem for T (x, t) and h(t) given by (2.7) -(2.10). We introduce the dimensionless variables, with the scales T s , x s and t s chosen as, On substituting (2.11) and (2.12) into (2.7) -(2.10) we obtain the non-dimensional form of the free boundary problem as, where the two dimensionless parameters ε and k are given by, . (2.17) Here primes have been dropped for ease of notation. The parameter ε measures the ratio of the time scale for heat conduction in the wax layer to the time scale for wax layer growth, whilst the parameter k measures the ratio of heat extracted from the wax layer, by cooling, to the heat the wax layer gains from the oil. This is a free boundary problem for u(x, t) with 0 ≤ x ≤ h(t), t > 0 and h(t) ≥ 0. It is worth noting that typical values, estimated in Schulkes [14] and Kaye and Laby [8], for the scales in (2.12), are T s ∼ 10 • C, x s ∼ 0.16m, t s ∼ 2 × 10 6 s ∼ 20 days, (2.18) with the dimensionless parameters,

The Free Boundary Problem [IBVP]
The free boundary problem associated with the mathematical model introduced in section 2 ((2.13) -(2.16)) may be written fully as, The problem, (3.1) -(3.6), will be referred to as [IBVP]. For any T > 0, the following subsets of R 2 are also introduced, namely, A reformulation of [IBVP] (with (R1) and (R2)) in terms of coupled integral equations is given by Friedman [5] and used by Schatz [13] and Cannon and Hill [3] to study the regularity of solutions to [IBVP]. It is established by Cannon and Hill [3] that with u : D ∞ → R and h : [0, ∞) → R being a solution to [IBVP], (with (R1) and (R2)) then, in fact, which requires (D3), in particular. A consequence of (R1) and (3.7) is also, Before proceeding to analyse [IBVP] further, we first consider steady state solutions associated with [IBVP].

Steady State Solutions to [IBVP]
A steady state solution to [IBVP] is a solution to [IBVP] which is independent of t. Thus, u s : , and satisfies the boundary value problem, It is now convenient to introduce F : R → R given by, where F (0) = 0 and, (3.14) Observe, via (D2) and (D3), that F (X) is strictly increasing with X, and F ∈ C 4 (R). Therefore the inverse F −1 : R → R exists, with F −1 ∈ C 4 (R), via (3.14). The boundary value problem (3.9) -(3.12) can now be written as, An integration of (3.15) gives, where A, B ∈ R are constants. Applying (3.16) and (3.18) we obtain, Finally, applying (3.17) and rearranging, we require, In a steady state we must have h s > 0. Therefore, it follows from (3.21) and (3.22) that we have established: It is worth noting that (3.22) may be written as, is the mean value of D(λ) over the interval λ ∈ k −1 , 1 . Also, we may rewrite (3.21) explicitly as An examination of h s = h s (k) in (3.23) establishes the following properties, ) and is strictly monotone increasing, We give a qualitative sketch of h s (k) against k in Figure 2. Next, an examination of (3.25) establishes the following properties of u s : [0, h s ] → R, It is worth noting from (iv) that inflection points will occur in the graph of u s (x) for

Qualitative Theory for [IBVP]
In this subsection we determine the principal structure and qualitative properties of the solution to the free boundary problem [IBVP]. To begin with, we have, Proof. Consider u : D ∞ → R on the compact region D T (for any T > 0). Then u is continuous on D T since u is continuous on D ∞ (via (R2)). Also, we have that u x , u t and u xx all exist and are continuous on D ∞ and so are continuous on D T (via (R2)). Also, we have from (3.1) that and so, Now suppose that u is not non-negative on D T . Then using (3.26)-(3.28), we may apply the Weak Parabolic Minimum Principle [18] to conclude that there exists a point (x * , t * ) ∈ ∂D T such that, via condition (3.2) and (3.29). Then, via (D2), we conclude that, However, since u achieves its infimum at (0, t * ) then u x (0, t * ) ≥ 0, which contradicts (3.30). Thus, we conclude that u must be non-negative on D T . This holds for any T > 0 and so u( However, via (3.26)-(3.28), it follows from the Weak Parabolic Maximum Principle [18] that there exists a point (x * , t * ) ∈ ∂D T such that, via (3.31). Now suppose that (x * , t * ) ∈ ∂D T L , so that (x * , t * ) = (0, t * ), and from (3.2) and (3.32) we have that, via (D2). However, since u achieves its supremum at (0, t * ) then u x (0, t * ) ≤ 0, which contradicts (3.33). Hence we conclude that (x * , t * ) ∈ ∂D T R , and so, via Therefore, u ≤ 1 on D T . Since this holds for any T > 0, we conclude that u(x, t) ≤ 1 for all (x, t) ∈ D ∞ . Consequently, we have, as required.
We next refine these inequalities in, and it follows from the Strong Parabolic Minimum Principle [18] that, and it follows from the Strong Parabolic Maximum Principle [18] that, where now, Also, via condition (3.2) with regularity condition (R2), we have from (3.34), and so, via (D2), as required.
As a consequence, we have the refinement, Proof. We introduce the functionv : with v : D ∞ → R as introduced in (3.53) and with µ ∈ R to be chosen. From (3.54) it follows that, Hence, after substituting (3.73)-(3.75) into (3.58), we obtain, Setting T > 0 and recalling (3.59) and (3.60), we have, We now choose, Thus, from (3.76) and (3.77) together with (3.72), we have, Then, via (3.71), we have thatv is continuous on D T and thatv x ,v t andv xx exist and are continuous on D T . Hence, with (3.78), we may apply the Strong Parabolic Maximum Principle [18] tov on It then follows, via (3.3) and regularity conditions (R1) and (R2) that, However, this contradicts (3.2). Therefore Since this holds for any T > 0, we have, as required.
The regularity condition (R2) with (D3) then immediately allows for, Corollary 3.7. Let u : D ∞ → R be a solution to [IBVP]. Then, In addition, we have, Proof. Let u : D ∞ → R and h : [0, ∞) → R be a solution to [IBVP]. Then from Proposition 3.6 we have, which is only possible when k > 1.
We now obtain bounds on h in the following, Next take any t > 0. It follows from the mean value theorem with (3.37), (3.38) and (3.40) that, with 0 <θ < 1. Then, via Proposition 3.6, (3.90) and (3.91), from which we obtain, via (3.3), However, via (3.2) and Corollary 3.7, Therefore, via (3.92) and (3.13), withū(x) monotone increasing for x ∈ [0,h], and, Further, the bounds obtained on u x , u t and h t , and consequently bounds on u xx together with the Ascoli-Arzela Theorem, allow for a deduction thatū x andū xx exist and are continuous, and moreover,ū : [0,h] → R andh must satisfy problem (3.9)-(3.12), and so are steady state solutions to [IBVP]. Henceū = u s andh = h s , as discussed in subsection 3.1. It is convenient to summarize the results in this subsection in the following, . We recall that the limit in (vi) is from below, whilst the limit in (vii) is from above. Also, we note that in physical terms the requirement that k > 1, for a solution to [IBVP] to exist, requires that the cooling process has to be sufficiently strong in order for the development of a wax layer to initiate. From Theorem 3.11 (i) -(v) we have obtained a priori bounds for [IBVP] on u : D ∞ → R, the partial derivative u x : D ∞ → R, together with h : (0, ∞) → R. Consequently, both global existence and uniqueness for [IBVP], can be anticipated by adopting an iterative approach to accommodate the quasi-linear terms (see, for example [9]), whilst following, in principle, Cannon and Hill [3]. We next develop the analysis of [IBVP] by considering the structure of the solution as t → 0 + and correspondingly as t → ∞.

Coordinate Expansions
We begin this section by analysing the structure of the solution to [IBVP] as t → 0 + . After which we consider the structure of the solution to [IBVP] as t → ∞.

Coordinate Expansions as t → ∞
We now consider the structure of the solution to [IBVP] as t → ∞. From Theorem 3.11, it follows that, as t → ∞. Thus we write, On substituting from (4.28) and (4.29) into [IBVP] we obtain the leading order problem forū and h, as We can eliminateh(t) from (4.30)-(4.33) to obtain, with λ ∈ C to be determined. After substituting from (4.38) into (4.34)-(4.36) we obtain the linear eigenvalue problem, with eigenvalue λ ∈ C. It is now convenient to introduce the function ψ : [0, h s ] → R defined by, together with the prescribed function ∆ : [0, h s ] → R given by, Substituting from (4.42) into (4.39)-(4.41) we obtain the equivalent eigenvalue problem, Henceforth, we will refer to the generalised Sturm-Liouville eigenvalue problem given by (4.44)-(4.46) as [S-L]. We note that [S-L] has the following properties: (S1) The eigenvalues of [S-L] are all real and may be written as, with λ r → ∞ as r → ∞.
We observe that (H4) and (5.22) are in agreement with (4.65). An implicit form for H(t) is obtained from (5.20) and (5.21) as Finally, having determined h 0 (t), we obtain from (5.11), with A(t) given in terms of h 0 (t) in (5.17). It is worth observing from (5.24) that, which, via (5.18) and (5.19), is monotonic decreasing from unity to k −1 with increasing t.

Numerical Solution to [IBVP]
In this section we consider numerical solutions to [IBVP] for comparison with the theory of the previous sections. For [IBVP] with constant diffusivity, the method of fundamental solutions was employed in [12] to provide numerical approximation of the solution to [IBVP] 3 . A useful feature of this method is that a node can be placed at (x, t) = (0, 0) in the domain to encapsulate the initialboundary conditions (3.5) and (3.6). However, for non-constant diffusivity, the partial differential equation (3.1) is quasi-linear, and consequently, the nonlinearity precludes numerical methods based on fundamental solutions. Hence, in the current situation, we employ an explicit finite-difference method to provide numerical approximations to [IBVP]. We note that although this method is simple to apply, in this setting there are several limitations, primarily due to the representation of conditions (3.5) and (3.6). We first transform [IBVP] to a rectangular domain, by introducing

It then follows from (3.1)-(3.4) that [IBVP] becomes
for X ∈ [0, 1]. Due to the degeneracy of (6.2) at t = 0, we set the initial conditions for the numerical method at t = δ > 0, with δ sufficiently small so that we can use the asymptotic forms for h and v in (6.6) and (6.7) at t = δ, respectively. We refer to the initial-boundary value problem given by (6.2)-(6.7) as [IBVP * ].
To implement the finite-difference method, we employed a uniform spatial grid with N x grid points to represent the interval [0, 1] so that the i th spatial grid-points X i = (i − 1)dX with dX = 1/(N x − 1). The temporal grid points t j , used to represent [δ, T ], were not uniformly spaced, with the time step dt chosen sufficiently small at each step to accommodate numerical stability conditions on the discrete evolution equations for v i,j ≈ v(X i , t j ) and h j ≈ h(t j ), namely .
We note that this local stability condition limiting dt relaxes, as h increases. The discretization of (6.5), (6.4), (6.2) and (6.3) respectively, in the order that h j+1 and v i,j+1 are computed, is given by: v Nx,j = 1; Using this numerical scheme we approximated the solution to [IBVP * ] in the cases when D(v) = 1 + cv(1 − v), with c = 3, 2, 1, 0, −1, −2 and −3. The parameter k(> 1) and ε(> 0) were chosen from k = 2, 10 and 20 with ε = 0.  Figures 4(a-c), we observe that as ε increases, with D and k fixed, the solutions to [IBVP], converge to the steady state u s at a slower exponential rate. This is further illustrated in Figures 5(a,b) and Figures 5(c,d), which graph h against t for each of the cases illustrated in Figures 4(a-f). Additionally, in Figures 4(d-f) we observe that as k increases, with D and ε fixed, firstly, that the steady state curvature changes from monotone to non-monotone as k increases past k = 2 and secondly, that as k increases, the time taken for u to approach the steady state u s decreases.

Comparison with Experiments
Detailed experiments on the formation of wax deposited layers in straight circular pipes transporting heated oil, when subject to external wall cooling, have been performed and reported by Hoffman and Amundsen [7]. In this section we make qualitative and trend comparisons with the theory presented in this paper and the experimental results in [7]. We observe that the experiments reported in [7] have fixed paraffinic wax and oil properties in all experiments, while the bulk oil temperature T o , the coolant temperature, T c and the oil flow Nusselt number N u are varied in turn, with a series of experiments performed in each case and wax layer evolving profile and equilibrium thickness measured. First, we relate the variations in T o , T c and N u , respectively, to the dimensionless parameters in the model, namely ε and k. We observe from (2.17) that, first k > 1 for wax layer formation to occur, and thereafter, (a) Increasing T o decreases k whilst keeping ε fixed.
(b) Decreasing T c increases k and ε together.
(c) Increasing N u decreases k whilst leaving ε fixed.
Thus, via Section 6, we can refer, for comparison, to the evolution of the dimensionless wax layer thickness with dimensionless time in Figures 5(a-d). Figure 5(a,c) correspond to cases (a) and (c) above. Conversely Figures 5(a,d) correspond to case (b) above. We first consider Figure 11 in [7]. This graphs the experimental wax layer equilibrium height against the wall temperature T c , at two different values of N u . Each graph has a critical value of T c , above which a wax layer does not form. This is consistent with the theory, corresponding to the critical value k = 1. As T c decreases from the critical value, the equilibrium wax layer thickness increases; this is born out for the theory Figures 5(a,c), where Figure 5(c) has larger k and ε values than Figure 5(a).
We now consider Figure 15 in [7], which graphs the evolution of the wax layer thickness with time, for a number of increasing oil temperatures T o . These profiles show remarkably similar structure to those theoretical profiles in Figures 5(a-d). The experimental graphs show a lowering of the profile with increasing T o . These can be compared with Figures 5(a,c), which show lowering profiles with decreasing k. This is consistent via (a). Note that the experiments show that a wax layer will not form for T o sufficiently large. This corresponds to k decreasing to k = 1 in the theory.
Finally we consider Figure 18 in [7]. This shows the evolution of the wax layer thickness with time, for a number of Nusselt numbers, with N u increasing. We see that the wax layer profiles are lowering with increasing Nusselt number. This corresponds to case (c) in the theory, and we compare with Figures 5(b,d). We see that, with ε fixed, the corresponding profiles are lowering as k decreases, in line with the experiments. To conclude, we note the striking similarity between the experimental profiles in Figure 18 of [7] and the model profiles given in Figure 5.

Discussion
In this paper we have developed and analysed in detail the simple thermal model for the development of a wax layer on the interior wall of a circular pipe transporting heated oil containing dissolved paraffinic wax, which was introduced in [12]. This approach is gaining considerable traction compared to the traditional mechanical and material diffusion theories; it is able to describe features associated with wax layer formation which have been absent from, or even contrary to, the outcomes from the mechanical and material diffusion theories. This view point is vindicated in a number of recent detailed reviews of the wax layer formation process, see, for example [10], [11], [8] and [17]. The current paper extends the theory developed in [12] to allow for the dependence of solid wax thermal conductivity on local temperature, which is a significant feature of solidified paraffinic wax, and which should be included as a principle component in the thermal model. This inclusion modifies the free boundary problem formulated in [12], which now becomes, in general, a nonlinear free boundary problem. Despite the introduction of this fundamental nonlinearity, we have still been able to develop a comprehensive theory for this improved model. The variable conductivity model was developed in Section 2, and formulated as a nonlinear free boundary problem [IBVP] in Section 3. A detailed and rigorous qualitative theory for [IBVP] has been developed in Section 3, whilst the small and large time asymptotic structure to the solution to [IBVP] was given in Section 4. In Section 5 we developed the asymptotic solution to [IBVP] in the situation when the heat conductivity time scale in solid wax is much faster than the time scale associated with the inter-facial formation of solid wax. These substantial analytical developments have been supplemented by an efficient and readily implemented numerical scheme for [IBVP] in Section 6. Finally, a qualitative comparison between the present theory and the experimental results in [7] was briefly considered in Section 7. These initial comparisons are very encouraging for the thermal mechanism model introduced in [12] and developed here. A more detailed experiment with model comparison to consider quantitative agreement appears to a worthwhile and appropriate endeavour to undertake.