Code Revival: Fluid Motion in a Curved Pipe

1. Introduction

The theoretical treatment of laminar flow in curved pipes was pioneered by Dean in the late 1920s [1,2]. He showed that in addition to a displacement of streamwise flow, a secondary flow also exists, taking the form of two vortices, the so-called Dean vortices. A constant D, the Dean number, was identified as the key parameter:

$$D=4Re\sqrt{{\displaystyle \frac{2a}{L}}},$$

(1)

where $Re$ is the bulk Reynolds number, a is the pipe radius and L is the radius of curvature of the curved pipe. Dean [1] obtained simplified equations of motion for the approximation $a/L\ll 1$, which can be interpreted as $a/L\le 0.01$ [2].

Dean [1] obtained theoretical solutions to the equations for the range $0\le D\le 96$ and following this, McConalogue and Srivastava expanded the solutions numerically to $96\le D\le 605.72$[3]. In this paper we treat a subsequent study by Greenspan and Schubert [4,5], where numerical solutions were extended to $D\le 5000$. This extension was done with a FORTRAN 66 code which is available as an Appendix in [4]. The purpose of our paper is to convert - and thus revive - this code to modern Fortran. The reason for this is twofold: (i) we need the revived code for planned future research and (ii) we want to make this code available to the scientific community for others who may need a fast and efficient code to calculate basic solutions of laminar flow in slightly curved pipes.

On a personal note, I have previous experience writing and running FORTRAN 77 code on an IBM mainframe at JET [6] during the period 1997-1998. Thus, the task at hand, to convert FORTRAN 66 code, was feasible.

The paper is structured as follows: In Section 2, we detail the FORTRAN 66 code revival steps; results from the revived code are presented in Section 3. A brief discussion is presented in Section 4 and finally we conclude in Section 5.

2. Code Revival

There were three aspects to the code revival, and they are treated in the three subsections below:

• Digitise the scanned document
• Convert the FORTRAN 66 code to modern Fortran
• Debugging and testing of the code

2.1. Digitisation

We note that the original FORTRAN 66 code is only available in [4] and not available in the official journal paper [5]. As can be seen, the scanned document [4] has a quite poor quality, making in challenging to perform the optical character recognition (OCR) needed. However, with 12 pages (∼600 lines) of code, it would be very cumbersome to retype the code entirely by hand

We did not attempt to improve the quality of the scanned document before conversion.

Our initial attempts at extracting characters from the scanned document was done using Python with EasyOCR [7] and PyTesseract [8]. However, these tools performed poorly on the document.

Our next attempt was with two online converters [9,10]. The results were better than for the Python tools, but still with a correct percentage of around 50%.

The best conversion, with a success rate of roughly 80%, was the Amazon Textract tool from Amazon Web Services (AWS) [11].

In the end, we used the AWS tool to generate the initial digital file and combined it with the output from the two online converters. And finally, a human correction was needed for around 10% of the code.

2.2. Conversion

Once the code was available in digital form, the main conversion from FORTRAN 66 to modern Fortran was to replace the Hollerith string and data notation [12]. To achieve this goal, we also used ChatGPT-4 [13] to convert the Hollerith commands to modern Fortran.

2.3. Debugging and testing

At this stage we had a complete code available, the final step being to debug and test the code.

We chose to use a laptop with Microsoft Windows 11 [14] installed. Microsoft Visual Studio [15] was installed as the graphical user interface along with the Intel Fortran Compiler (ifx) [16]. In terms of Fortran versions, ifx uses the latest standards from Fortran 2018 and some Fortran 2023 features. In that sense, we consider the updated code to represent modern Fortran.

What remained was to compile (build solution) and systematically resolve the errors. Some sections were updated, for example replacing the reading of and writing on punch cards with files. Another major change was to make an IF statement to loop over all D values analysed in [4,5].

The final result of the steps is the modern Fortran code provided in Appendix A.

3. Results

All post processing in this section was done using GNU Octave [17].

3.1. Comparison to Original Findings

We include two figures with contour plots of $\varphi$ and w, Figure 1 and Figure 2, to facilitate a direct comparison to Figure 3 and 4 in [5]. Here, $\varphi$ is the non-dimensional stream function for the secondary flow and w is the non-dimensional velocity in the streamwise direction. Only the upper half of the pipe is calculated by the Fortran code, we mirror this to the lower half for clarity of exposition. 0 (180) degrees indicate the outboard (inboard) positions, respectively.

The results match the original results very well, both in position and magnitude. Again, the $\varphi$ curves move clockwise up to $D=500$ and then counter-clockwise for larger D. And a certain waviness is observed for the core w curves for $D=5000$, which was speculated to perhaps indicate a beginning transition to turbulence [5].

3.2. Improvements

Having a complete modern Fortran code available enables improvements to e.g. the computational grid resolution. The original resolution was $\mathrm{\Delta }r=0.1$ and $\mathrm{\Delta }\varphi ={\displaystyle \frac{\pi }{18}}$. In the code in Appendix A, these resolutions are named NR and NA, respectively. We updated the code to run with double resolution, both for r and $\varphi$. To distinguish the grid resolutions, we call them "original" and "updated".

We compare original and updated contour plots for $\varphi$ and w in Figure 3 and Figure 4. The two largest values of D are compared; qualitatively, the results are very similar, although the contours of the updated resolution are more smooth than those of the original resolution.

A topic of interest for our future research is how the peak w value ($w_{\max }$) moves as a function of D. This is shown in Figure 5, where both the position and magnitude of $w_{\max }$ are included. The values are compared for the original and updated resolution; we observe that there is a substantial difference in the maximum position, but not in the maximum value of w. For reference, $w_{\max }$ values from [3] are included, which agree with the present values for smaller values of D, but deviate for larger D.

An important effect of the secondary flow is the blockage effect, i.e. that it obstructs the streamwise flow. This is also calculated in the Fortran code and shown in Figure 6. Here, the flux (or flow rate) is normalised to that of a straight pipe.

4. Discussion

It has been discussed whether flow in curved pipes only depends on the Dean number, or if one should separate studies into separate $Re$ and curvature effects. Now it is indeed clear that such a separation is important for turbulence studies, see e.g. [18]. For the case we treat with slightly curved pipes, the transition to turbulence is subcritical, which is the same as for straight pipes.

From experimental investigations with higher curvature [19], it is interesting to note the persistence of structures in turbulent flow similar to those identified in laminar flow. Perhaps the laminar flow structures form a skeleton/seed which turbulence builds on?

5. Conclusions

In this paper, we have presented a case study on how to revive older code (FORTRAN 66) to modern Fortran. This code revival has enabled fast calculations of laminar flow through slightly curved pipes. The basis was [5], with the original code being provided in [4]. We hope that our work serves as a template and inspiration for other scientists who have a need for bringing older code back to life again. Also, the modern Fortran code in Appendix A is now available to the scientific community as a fast method to obtain a physical picture of streamwise and secondary flow in curved pipes.

After the code was operational, we tested increasing the grid resolution and found that the position of the maximum streamwise velocity as a function of D was sensitive to this change.

The code can be modified in other ways to e.g. extract information on the stream function for the secondary flow or to make a more detailed study using smaller steps in D.

We plan to use the code in the future to facilitate interdisciplinary work on fluids and plasmas which was initiated in [20].