On implementation of a semi-analytic strategy to develop an analytical solution of a steady-state isothermal tube drawing model

In this paper, we consider an isothermal glass tube drawing model consisting of three coupled nonlinear partial differential equations. The steady-state solution of this model is required in order to investigate its stability. With the given initial and boundary conditions, it is not possible to determine an analytical solution of this model. The difficulty lies in determining the constants of integrations while solving the second order ordinary differential equation analytically appearing in the steady-state model. To overcome this difficulty, we present a numerical based approach for the first time to develop an analytical solution of the steady-state isothermal tube drawing model. We use a numerical technique called shooting method to convert the boundary value problem into a set of initial value problems. Once the model has been converted into a system of differential equations with initial values, an integrating technique is implemented to develop the analytical solution. The computed analytical solution is then compared with the numerical solution to better understand the accuracy of obtained solution with necessary discussions.

In this paper, we consider an isothermal glass tube drawing model consisting of three coupled nonlinear partial differential equations. The steady-state solution of this model is required in order to investigate its stability. With the given initial and boundary conditions, it is not possible to determine an analytical solution of this model. The difficulty lies in determining the constants of integrations while solving the second order ordinary differential equation analytically appearing in the steadystate model. To overcome this difficulty, we present a numerical based approach for the first time to develop an analytical solution of the steady-state isothermal tube drawing model. We use a numerical technique called shooting method to convert the boundary value problem into a set of initial value problems. Once the model has been converted into a system of differential equations with initial values, an integrating technique is implemented to develop the analytical solution. The computed analytical solution is then compared with the numerical solution to better understand the accuracy of obtained solution with necessary discussions.
In glass industry, tubes are drawn through various manufacturing processes used to achieve the continuous production of glass tubes having correct wall thickness and diameter. The most commonly used are the Danner process and the Vello process [1][2][3][4][5] having great importance in glass fabricating industry and are still in use today. In Danner process, glass is melted in a furnace to the stage where it is soft and pliable. Molten glass is then let to fall with low feeding speed v 0 on the surface of a cylindrical device called mandrel kept in a temperature controlled tank called oven. Mandrel is slightly inclined and hollow such that the air can be blown through it. By continuous rotation of mandrel about its axis of symmetry, the molten glass falling downward creates a smooth layer around the mandrel. It cools down gradually and takes the shape of a thick-walled hollow glass tube with all desired properties at just below the end of mandrel. The length of hot-forming zone is taken as L. It is then pulled out by a drawing machine with a drawing speed v L > v 0 . This ratio v L /v 0 > 1 is called the draw ratio. Keeping a constant temperature in the hot-forming zone leads to develop an isothermal tube drawing model. The drawn tube is then conveyed straight by rollers to further process of cutting, finishing, polishing and packaging at the end of the spinline. This manufacturing process is explained and illustrated in 3 . For the Vello process, we refer to [1][2][3][4][5] .
The geometry of the drawn tube is illustrated in Fig. 1. The length of hot forming zone is denoted by L. Inner and outer radii, inside pressure, feeding and take up speeds, axis of symmetry of the tube during the production process are shown in Figure. All other parameters involved along with their numerical values are illustrated in Table 1.
The shaping parameters such as the wall thickness and cross-sectional area (or diameter) are the main characterizations of the drawn tube. In either of the manufacturing processes, the required shape of the tube can be maintained by the stream of the air gently blown through the mandrel. Insufficient quantity of the air blown through the mandrel can disturb the desired shape of the glass tube. As a result, this insufficient balance collapse the walls of the glass tube. Moreover, the geometry of the glass tube can be controlled by the parameters involved herein such as the glass temperature, the composition of the raw material, the pressure of the blowing air in the mandrel and the rate of draw. We have added a Fig. 1 to better understand the geometry of the drawn glass tube.
A variety of mathematical models both for the isothermal and non-isothermal tube drawing processes, with different levels of descriptions and needs, are available in the literature e.g., see 1-4 and 6-11 . In 1-3,7 numerical solutions of both the isothermal and non-isothermal models have been found and used to optimally control the geometry of the glass tube. In recent years, a variety of mathematical models representing physical and real world problems have been investigated for numerical solutions and stability analysis (for example see [12][13][14][15][16][17][18][19][20][21][22][23][24] ). Thus, the concept of finding numerically accurate and exact solutions of real world problems has attracted attention from all over the world.
In this paper, we have considered an isothermal glass tube drawing model consisting of three coupled nonlinear partial differential equations of first and second order. The steady-state numerical solution of this model is required in order to investigate its stability. Recently, we have analyzed the stability of an isothermal tube drawing model and performed a complete mathematical analysis in 3 by incorporating a steady state numerical solution. To better understand the physical model and to examine the accuracy of the obtained numerical solution, our objective here is to develop an analytical solution of the steady-state isothermal model. To start with, it is not possible to determine an analytical solution of this type of model with the given initial and boundary conditions. When we solve the second order ordinary differential equation appearing in the steady-state model, the difficulty lies in determining the constants of integrations. To overcome this concern, we utilize a numerical technique called shooting method to convert the boundary value problem into a set of initial value problems. An integrating technique is implemented to develop an analytical solution of the converted initial value problems. www.nature.com/scientificreports/ Similar approach has been suggested in 10,22 to construct an analytical solution of a steady-state melt-spinning model. The semi-analytical approach developed in this paper has an advantage over other numerical methods that it may be applied directly to both linear and nonlinear systems of equations without the need for discretization, linearization, or perturbation. The organization of the paper is as follows: In "Governing equations" section, we give a brief description of the mathematical model of an isothermal tube drawing process. A strategy converting a boundary value problem into a set of initial value problems is described in "Numerical approach" section. Analytical solution is developed in "Analytical solution" section and compared with the numerical solution in "Accuracy comparison" section. Conclusion is given in "Conclusions" section. We have included some leftover parts of the modeling of tube drawing process in the "Appendix". In this section, we have included one example concerning the application of the analytical solution obtained through the implemented semi-analytic technique.

Governing equations
In the literature, different types of models for the drawing processes with different level of demands and descriptions are available. A considerable amount of work has been carried out by different researchers 1-4,6,8,9 and 11,25-28 to model the tube drawing process. In this section, we briefly explain the mathematical model for an isothermal tube drawing process.
To model the tube drawing process, we consider an incompressible Newtonian flow of a molten glass between two free surfaces r = r 1 (z, t) and r = r 2 (z, t) where r 1 (z, t) and r 2 (z, t) respectively denote the inner and outer radii of the glass tube, and assume that temperature remains constant throughout the forming zone. Glass tube during production process is illustrated in Fig. 1. In the draw-down zone, the surface tension force and the inertial force acting upon the molten glass are insignificant and hence can be neglected. This kind of flow is governed by the equations The Eqs. (1a)-(1c) are taken from 4 and are known as the standard equations showing the axi-symmetric stokes flow. The first equation is the continuity equation and the last two equations are the momentum equations respectively in r and z directions. We denote the derivatives by subscripts r and z where z denotes the distance along the axis of the glass tube and r measures the distance perpendicular to it. The velocity of the molten glass is defined to be v = (u, v) where u and v are the components of velocity v along z and r direction respectively. The pressure, density and the acceleration due to gravity are denoted by p, ρ and g respectively. Now, at the free surfaces r = r 1 (z, t) and r = r 2 (z, t) , it is necessary to specify the stress conditions and the kinematic conditions.
On the inner and outer surfaces of the glass tube, the stress conditions are given as: where n i and n o are the unit normals on the surfaces r = r 1 and r = r 2 of the tube respectively defined as and p s is the inside pressure applied on the surface r = r 1 of the glass tube and τ is the stress tensor given as The kinematic conditions are described as www.nature.com/scientificreports/ and thus the stress conditions (1d)-(1e) can be expanded to give where µ denotes the viscosity of the molten glass which remains constant at the given temperature.
To take benefit of the small parameters present in the problem, it is now convenient to convert the Eq. (1) into dimensionless form. For the dimensional quantities, the appropriate scaling is defined as: W is the width of the glass tube that has a very small value than the typical length of hotforming zone L, U is the typical drawing speed and µ 0 is the typical melt glass viscosity which has a constant value for the isothermal case. After dropping the bar notation, the system of governing equations in dimensionless form is given as follows: where St = ρgL 2 µ 0 U is known as Stokes number.
The stress conditions (1h)-(1k) read as The kinematics conditions (1f)-(1g) become The above model can be simplified by means of an asymptotic expansion, in which the inverse aspect ratio ε is used as scaling parameter. Assuming the glass flow as a thin layer flow and ignoring the large aspect ratio of the flow, one derives the simplified equations to model the isothermal tube drawing process. In this derivation, the surface tension and inertial forces acting upon the molten glass have also been ignored due to their insignificant www.nature.com/scientificreports/ contributions. For details, we refer to 1-6 and 9,27 . However, we have included the remaining considerations of the modeling process in "Appendix A" and reach at the following model equations representing the isothermal process: equipped with the initial conditions and the boundary conditions where A 0 and R 0 are the cross-sectional area and average radius of the glass tube at the time of entering the hotforming zone, respectively. Acceleration due to gravity and density of the molten glass are denoted by g and ρ , respectively. Average radius R of the tube is defined as Equations in system (3) give us cross-sectional area A, velocity v, and average radius R of the glass tube. Width W of the tube can be determined by the equation Since the temperature throughout the forming zone of the tubing process remains constant, the viscosity µ of the melt glass also remains constant.
Dimensionless form. We introduce the following dimensionless quantities into the model (3) to get the dimensionless equations where we have removed the asterisk notation and are the dimensionless parameters.
The initial and boundary conditions are respectively transformed to The steady state form of the model (4) is written as www.nature.com/scientificreports/ subject to the conditions: The system (5) is a system of coupled nonlinear ordinary differential equations. We are interested in developing an analytical solution of the nonlinear system (5). However, we have already proved the existence and uniqueness of the solution of steady state isothermal tube drawing model (5) in 3 . For details, we refer the reader to 3 .

Numerical approach
There are different techniques available in the literature to solve boundary value problems analytically. The obtained general solutions contains arbitrary constant which can be determined by using boundary conditions. However, there exist some problems (e.g. problem (5) and see 10,22 ) where it is not possible to find integration constants with the given conditions. It is interesting to mention that integration constants can be found if boundary value problem (BVP) is transformed to an initial value problem (IVP).
Shooting method is a technique that transforms a BVP of the form into an IVP where s is a guess to be determined in such a way that it shoots φ(b). We apply RK4 method to the second order differential Eq. (6) to determine approximation s i , i = 1, 2, . . . for φ(b) . Let the first guess s 1 be taken as RK4 method with s 1 gives us an approximation β 1 for φ(b) . If the absolute error |β 1 − φ(b)| is less than some pre-assigned tolerance, we stop, otherwise we refine our guess by considering and compute β 2 by RK4 method. If the error |β 2 − φ(b)| agrees to pre-assigned tolerance, we stop. In the case of disagreement, we make further guesses using the secant formula where is the error of approximation. RK4 method is continued each time with a new guess s i , i = 3, 4, . . . computed from (7) until β i agrees with the value φ(b) to pre-set tolerance.
Once we able to compute a most suitable guess s i → s , we are able to transform the BVP to IVP, and hence an analytical solution of the problem of type given in (6) can be found.

Analytical solution
With the initial conditions (5d), the Eq. (5a) is solved to get Equation (5b) can now be put in the form along with the conditions St, www.nature.com/scientificreports/ With the usual techniques, it is not an easy task to develop an analytical solution of the BVP (9). We use a non-conventional approach to build an analytical solution of the given nonlinear model. We convert the BVP (9) into two IVPs. For this purpose, we set and Then the IVPs corresponding to BVP (9)  Equation (16) along with c 1 and c 2 respectively given in (13) and (17) represents an analytical solution of Eq. (5b).
Next we put the steady state solutions for A and v in Eq. (5c) to get where u = R 2 . Equation (18) is a well-known Riccati's equation 29 whose solution is supposed to be of the form www.nature.com/scientificreports/ where one solution u 1 , guessed by hit and trial procedure, is given as and y is any nonzero unknown function of z.
Differentiating (19) with respect to z and then putting the values of u and du dz in Eq. (18), we reach at which is a first order linear ordinary differential equation in y. After few steps of calculations, the solution of this differential equation is found to be where c 3 is a constant of integration, and A(z), B(z) are given as Using (19) and the condition (5d), we get where and Substituting c 3 in (20), we obtain Therefore, where and hence is an analytical solution of the Eq. (5c) along with c 1 and c 2 given respectively in (13) and (17). Thus, the solutions given in (8), (16) and (22) together constitute an analytical solution of the steady-state model (5). .

Accuracy comparison
The developed analytical solution is plotted and compared with the numerical solution in Figs. 2, 3, 4. The parametric values used in the model are given in Table 1. The solutions obtained by RK4 method are reliable and competitive with the analytical solutions. We observe that the errors between numerical and analytical solutions for the steady state variables A, v and R are negligibly small (see Table 2). We also determine W, the width of the glass tube, by using the relation A = 2πRW . The steady state analytical and numerical solutions for W and the corresponding error are also shown in the Fig. 5. It is observed from Table 2 that error between the numerical and analytical solution for A, v, and R is maximum at z = 0.2 , z = 1.0 , and z = 0.1 respectively. Moreover, the error for the solution W is maximum at z = 0.1.

Conclusions
In this work, we have presented a numerical based approach to find analytical solution of a boundary value problem where one cannot find constants of integration with the given boundary conditions. The approach is to convert a boundary value problem into initial value problems and then to develop an analytical solution of the resulting problems with the usual methods. To explain the approach, we have developed an analytical solution of a steady-state model of the isothermal tube drawing process with the help of shooting method. The obtained analytical solution is almost in agreement with the numerical solution that justifies our approach. www.nature.com/scientificreports/ Solving the Eq. (25) for G(z, t) and p 0 , we get By using the kinematic boundary condition (23d), the Eq. (24) at boundary r = r 1 is written as Now using (26), we obtain Similarly the Eq. (24) at boundary r = r 1 gives us Equation (29a) gives us the inner radius r 1 and the Eq. (29b) gives the outer radius r 2 of the glass tube. We do some manipulations with the Eq. (29) to obtain a differential equation for mean radius R in a dimensional form that is, where R in terms of r 1 and r 2 is the mean radius of the glass tube and is given by the relation and A = 2πrW is the cross-sectional area of the glass tube where W is the width of the glass tube defined by the relation Equations (29) yield us the continuity equation In dimensional form, the Eq. (32) becomes as Since u 0 and µ are not depending or r. Therefore, the leading-order r-momentum equation (23b) yields which mean that the leading-order pressure is not a function of r.
We consider the z-momentum equation to get equations in closed form for the tube drawing process.
Multiplying the Eq. (34) by r and then integrating from r = r 1 to r = r 2 , we get The boundary conditions of order ε 2 for the normal stress are . Differentiating p 0 with respect to z and using the conditions (37), the momentum equation (35), after certain simplification, becomes or in a dimensional form For simplicity, we dropped the zero subscripts of the leading order quantities present in the derived Eqs. (30), (33) and (38). Moreover, we need to attach these derived equations with the initial and boundary conditions.
The initial conditions read as: The boundary conditions are: The derived model Eqs. (30), (33), (38) along with the initial conditions (39a) and the boundary conditions (39b) are highly nonlinear, coupled and describe an isothermal tube drawing process. To make it more readable, we replace u by v. Then the model for the isothermal tube drawing process in a dimensional form is given as: with the initial conditions and the boundary conditions where L is the length of the hot-forming zone, p s is the inside pressure, µ is the input viscosity, v 0 is the feeding and v L is the drawing speed.
The system (40) gives us the velocity v, the cross-sectional area A and the average radius R of the glass tube. The width, denoted by W, of the glass tube can be determined by the relation Throughout the forming zone of the tubing process, the temperature remains constant and hence the viscosity µ of the molten glass is also remains constant. (42) y(x) = 2 3 (x + 1) 3/2 + 4 3 .
(43) 2y ′ y ′′ = 1, y(0) = 0, y ′ (0) = w 0 ≈ 1,  www.nature.com/scientificreports/ Case-II: when the exact solution is impossible with the given BCs. In this case, it is not always possible to determine an analytical solution of a system of equations with the given initial and boundary conditions. As an example, we refer the reader to 22 . absolute error x Figure 8. Error between the derived semi-analytic solution and that obtained through the shooting method using discretization step size h = 0.5.