{'fraction_non_alphanumeric': 0.030845720663588858, 'fraction_numerical': 0.006595302115569987, 'mean_word_length': 4.680691642651297, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast computing of velocity field for flows in industrial burners and pumps

In this work we present a technique of fast numerical computation for solutions of Navier-Stokes equations in the case of flows of industrial interest. At first the partial differential equations are translated into a set of nonlinear ordinary differential equations using the geometrical shape of the domain where the flow is developing, then these ODEs are numerically resolved using a set of computations distributed among the available processors. We present some results from simulations on a parallel hardware architecture using native multithreads software and simulating a shared-memory or a distributed-memory environment.


Introduction
Design, development and engineering of industrial power burners have strong mathematical requests, such as numerical resolution of differential problems involving Navier-Stokes equations for velocity and pressure fields, geometrical design of combustion heads for a correct shape and optimal efficiency of flame, geometrical design of ventilation fans and computation of correct air inflows for optimal combustion, computation of flows in the internal vanes of oil pumps for estimating the correct mass rate and dimensions of gears. Rapid prototyping for an accurate design of the correct geometries involves a numerical simulation of the gas or oil flows in these burner's components. The necessity of an high graphic resolution for the visualization of the flows requires the knowledge of a large amount of streamlines. Therefore the numerical computation is very onerous about amount of memory and power of processors of the used hardware environment. A technique for reducing the total time of computation can be very useful. In this paper we present a method based on reduction of Navier-Stokes system to a set of not coupled ordinary differential equations. Each of these equations describes the fluid velocity field along each flow streamline. The advantage of this method, compared to an usual one such as finite differences scheme (see e.g. [1]), is a faster computation of the numerical resolution and an easy parallelization in a multicore environment, using a distribution of blocks of streamlines among the available processors.

The mathematical model
In this section we briefly present the mathematical model developed for simbolic and numerical resolution of the differential problems about some flows in burner's components. The main set of equations is the Navier-Stokes system describing the dynamics of the fluids (see e.g. [2]): where v is the velocity field, p is the pressure, ρ the density and µ the dynamic viscosity of the gas or oil. Sometimes there are other equations to be considered, such as continuity equation, but here we would describe a general simplified case. A geometric characterization of the flow is offered by the streamlines, that is the continuous lines on which the tangent at any point is parallel to the velocity vector. At every instant of time, the velocity field in the domain of the flow is the set of vectors which give the velocities of the particles at every point (eulerian description of a flow). The computation of this vectorial field is done by numerical discretizations of (1), e.g. using finite temporal steps, which are usually expensive for CPUs and central memory of an hardware environment.
The fluid motion can be described by the set of all the trajectories followed by the particles and the speed along these geometric lines in a given interval of time (lagrangian description of a flow). The computation of the trajectories usally requires the knowledge of the velocity field, therefore numerical methods are expensive too (see [3] for a review on computational fluid dynamics techniques). We present a method for computing the velocity field by reduction of Navier-Stokes equations to ordinary differential form.
is a time value between initial t i and final instant t f , and r φ ∈ [r 1 , r M ], the curve φ : s −→ Φ(r φ , s, t 0 ) is a C k ([s 0 , s 1 ], R 3 ) parameterized invertible representation (see [4]) of a streamline at time t 0 . The r parameter identifies a single streamline, and in numerical computation it takes integer values in interval 1, M , where M is the total number of streamlines. Then temporal evolution of the flow is described by the map t −→ Φ(r, s, t), for (r, s) ∈ [r 1 , r M ] × [s 0 , s 1 ]. At any t 0 ∈ [t i , t f ], the map (r, s) −→ Φ(r, s, t 0 ) is a snapshot of all the streamlines at time t 0 . We can consider the velocity field of this snapshot as a set of vectors depending only from (r, s) and not from t, so that for the instant t 0 +dt the field is computed using the new snapshot s −→ Φ(r, s, t 0 + dt). For fixed t 0 and r φ , let (x i ) 1≤i≤3 be a cartesian coordinates system in R 3 and s the parameter of the representation s −→ φ(s) = Φ(r φ , s, t 0 ). In general, under usual hypothesis of regularity, the relation between a streamline and its tangent vector field is whereφ i = dφ i ds . Now we define the parameterized velocity field along the streamline as so that u = u(r φ , s, t 0 ) = u(s). Using chain rule we haveu i = 3 m=1 v imφm and, from (2), the following relation holds: in Navier-Stokes system, so it can be substituted by v iui (φ i ) −1 . Now let g = g(x) be the inverse of φ and q = p • φ; using standard computations we obtain Along the streamline s −→ φ(s) we have ∂ t v = ∂ t (u • g) = 0, so Navier-Stokes equations can be translated into the following ordinary differential equations for the functions (u i ) 1≤i≤3 : where ρ and µ should be considered as function of the parameter s and the derivatives of g should be expressed throught those of φ. These equations are again nonlinear, and a possible symbolic resolution required non standard analytical methods (as Lie symmetries, see [5]), but now all the derivatives are ordinary and the equations are no more coupled.

Example: internal gears pump
As example of technological interest, we consider an internal gears pump for oil, usually formed by a central internal rotor (pignon) which moves an external one (crown).
The oil is contained in the geometrical zones (vanes) between the two rotors. The knowledge of the flow in these vanes is important for a comprehension of some critical phenomena as cavitation (bubbles incoming) or loss of lubrication. We would study the flow in the critical zone of oil compression, where the rotors are very near and the vane has a small volume. Here we consider, for simplicity, that streamlines have the same geometrical shape of the rotors. In the case, e.g., of a classical teeth gears pump, without considering a dimensional factor, a single profile (left side of the tooth) can be described by the following representation where π + arctan(s 0 ) ≤ s ≤ 3 2 π, for an opportune value s 0 . The evolution in time can be handled, e.g., by a rotation of angle θ = ωt for an interval of time small enough to have a costant shape for the volume between the two rotors. Then in this interval the streamlines have a new parametrization given by

The computational algorithm
For computing a solution of (8), also in case of more complex geometries than previous or in case of turbolent flow, we have tested an algorithm based on smooth interpolation of streamlines which are computed by a first rough numerical resolution of original Navier-Stokes equations. For example, in the oil suction channel (see fig. 4) there is a geometric angle where the flow can become turbulent and problems of cavitation (formation of air or vapours bubbles due to local gradient of pressure) can happen.
We describe now the algorithm. For a chosen set of time values, a suitable number of streamlines are computed using a finite difference scheme on a grid covering the geometrical domain of interest (see [2]). Then each streamline is divided into couples of geometrical consecutive points, i.e. {P 1 , P 2 }, Every couple is interpolated with a cubic polynomial p = p(t), t ∈ [0, 1] (t is not time), imposing the following four analytical conditions: where 1 ≤ k ≤ N − 1 and v is the numerical solution (velocity field) of the original Navier-Stokes system (see [6] for details). In this way we obtain a set of class C 1 new streamlines. Note that cubics are the smallest degree one variable polinomials useful to have continuous differentiation along the interpoled lines. Each new streamline can be represented using, e.g., a unique parameter s ∈ [0, 1] with the re-parametrization s = (m+t) is a cubic between two generic points, the curve φ = φ(s) has derivativeφ = N df dt (13) and two consecutive cubics have the same derivative in their common point. The first order derivativeφ is necessary for the ODE form (8) of the Navier-Stokes equations. It can be shown (see [6]) that, if M is the number of streamlines which we want to interpolate, all the 4 × (N − 1) × M coefficients of the cubics can be computed by N − 1 matrix-vector products Ab where A is the 4M × 4M matrix Note that A is a sparse matrix with density number (see [7]) at most 1 M . Also, A is a block diagonal matrix where each block is constant, so the execution of a product like Ab can be optimized using standard techniques (see [7]).

About numerical resolution
In this section we discuss the numerical resolution of the nonlinear differential equations (8). For physical reasons, as typical in eulerian point of view, we might resolve an initial value problem for every streamline and for every instant of discretized time, imposing two known values u 0 and u 1 for u at s = s 0 and s = s 1 and two known valuesu 0 andu 1 foru at s = s 0 and s = s 1 respectively. The amount of total computations for this method depends on the number of streamlines, on the number of time values considered, and on the spatial step size (that is, degree of discretization for the independent variable interval [s 0 , s 1 ]) of the algorithm used for the resolution of a single equation. In addition, for complex geometries or turbulent flows, there is the roughresolution of the classical Navier-Stokes system and then the interpolations of the streamlines so found. Using an hardware environment for numerical computations, we have tried to limit the total time spent on this work and the amount of total memory required for handling the data, so that one can use typical workstations usually devoted to technical and scientific computing. In particular, we have used a system with 2 dual-core processors (3.2 GHz as clock frequency) and 4 GByte of central memory. We have prefered the software Mathematica, version 5.2, for the general numerics and its internal function NDSolve for the numerical resolution of differential equations.
The function NDSolve (see [8] for details) has the option Method for specify the desired numerical scheme and the option StartingStepSize for specify the initial grid step for the independent variable. This step is then automatically adapted by an internal procedure if necessary. We have chosen for Method the Adams-Moulton or Backward Differentiation schema for their reliability when used with nonlinearity as uu type (see e.g. [9] for technical informations and [10] for some examples). The application of NDSolve required the specification of the initial values u andu respectively for velocity and acceleration along a streamline. The numerical resolution of associated algebraic systems is then executed using advanced methods based on Lapack technology (see [11]). Version 5.2 of the software Mathematica provides support for multicore computations, spreading a single massive operation into a suitable number of parallel threads among the available processors (see [12]). In particular the numerical linear algebra operations seem to have the major benefits from multicore support. Therefore we have tested these features for the matrixvector products for the computation of all the coefficients of the cubics in the case of complex shape streamlines, and for the numerical resolution of the nonlinear differential equations, where the use of implicit discretization schema involves matrix computation and resolution of many systems of algebraic equations. Also, Mathematica has the function BlockMatrix of the package MatrixManipulation for construction and manipulation of block diagonal matrix A from block T, and the multiplication of sparse-dense matrices like Ab is highly optimized (see [13]).
In Mathematica all the calculations are executed by a central logical unit called kernel, which distributes parallel threads of a single computation among the processors or cores available. So we have measured the performances of a single kernel spreading computations on 4 cores (2 physical processors), simulating a shared-memory system. As many kernels can run indipendently on the same machine, we have also tested the performances of 2 kernels running computations on 2 physical processors, but in this case using suitable procedures we have previously equi-distributed, among the two units, the data to handle. The two kernels have independently addressed their own memory allocation, simulating a distributed-memory system. For comparison, we have executed the same computations on a single onecore processor (3.2 GHz) machine with 4 GByte of memory. Table 1: Speedup for computing cubics coefficients, 100 couples of points for each of M streamlines. Right column shows speedup results in the case of 2 kernels with 1 dual processor for each of them (distributed programs), while central column shows results in the case of 1 kernel with 2 dual processors (multicore program). Multicore single program presents better performances when amount of data increases.

About performances
The Mathematica internal function AbsoluteTiming gives the total elapsed time spent by the kernel for the execution of a single task. Therefore this time amount is the sum of the time spent by CPUs and the time spent by program (and operating system) for handling data among threads. For example, the command AbsoluteTiming[A*b;] gives the total elapsed time of execution of all the necessary operations for the calculation of the product c=A*b, without front-end visualization of the result (note the use of ; in the expression), but with product vector c stored into memory. In this way, if a task is spreaded by the kernel among parallel threads, the output given by this built-in function is the effective total time spent for completion of calculations, and it's not the sum of the single times spent by single threads.
For an evaluation of the performances of the multicore technology, we consider the notion of speedup S = T 0 T , that is the ratio between the elapsed time spent by a single processor session and the elapsed time spent by a multicore session for the execution of the same task. Now we discuss the results obtained for interpolation of streamlines in complex geometry ( fig. 4) and for resolution of corresponding flow differential equations. In Table 1 we report the performances registered for the computation of M streamlines interpolated by cubics, in the case of 100 couples of points for each of them. Note the good performance of the multicore technology on single kernel case when one increases the number of streamlines, while Table 2: Speedup for computing velocity field, N = 100 pairs of points for each of M streamlines.
Right column shows speedup results in the case of 2 kernels with 1 dual processor for each of them (distributed programs), while central column shows results in the case of 1 kernel with 2 dual processors (multicore program). Multicore single program has better performances than distributed parallel programs. the other case offers a slower speedup due to time spent in spreading and collecting data among the two kernels. The efficiency shown by the single kernel method when the size of the numerical parallel task increases is an experimental proof of the so called Gustafson generalization of Amdahl Law (see e.g. [14]).
When solving differential equations, if N is the number of couples of points used for interpolation, the number of ODEs for every velocity component is N too, because the cubics coefficients are not the same for two consecutive pairs. So the computed numerical value of the solution at the end point of the k-th pair is the initial value for computing the solution of the (k + 1)th pair. In this way, using a built-in function like NDSolve, the numerical computation of the solution along a streamline becomes an essentialy serial task. For evaluating the effect of multicore Mathematica technology in the case of vectorial operations, we have tabulated all the coefficients and initial values in a N dimensional array inputs and applied NDSolve to it using a vectorial single instruction like Map[modifiedNDSolve,inputs], being Map a built-in optimized iterator and modifiedNDSolve a suitable modified version of NDSolve which accepts as input the cubics coefficients. This algorithm is mapped to the array of all the M streamlines, assigning at each step as initial values for u andu the end values computed in the previous step. In particular, for N = 1 we have the case of a priori assigned streamlines family, as that of the flow in gears pump. Table 2 reports the performances for the numerical resolution of ODEs. The speedups are smaller than those of the previous case, probably for the greater percentage of not parallelizable operations required by the used algorithm, but again the performances are better for the multithreads single program. Note that the virtual distributed-memory environment fails when the size of the problem reaches a critical value, while the virtual shared-memory one can complete the task, probably for a better management of memory and centrally cached informations. Note also that the elapsed time necessary for a full calculation of velocity field by a single kernel on a single processor has a range, depending on M , variable from few seconds to tens of minutes, therefore for great sizes or large physical domains the speedup of a multicore computation could help in a significant way.

Conclusions
From the discussion about the used algorithms and registered performances, we can deduce the following conclusions on the parallel resolution, in a multicore environment, of blocks of ordinary differential equations: • the multicore technology gives better speedups in the case of a shared-memory environment with a single central multithreads program which executes all the computations; • on a multicore machine, from the last line of table (6) we deduce that shared-memory single kernel resolution can execute the computations on a greater amount of data than distributed-memory multi-kernel resolution, which requires a greater allocation of memory; probably this fact depends on the use, made by shared-memory method, of a central caching of temporary data instead of a distributed one; • mathematical computations on large amount of scientific or technical data can take advantage from the use of numerical software specifically designed for a central multicore architecture, as already shown in previous very particular experiences (see e.g. [15]); • the streamlines algorithm and a multicore computation can help to speedup the search, in a suitable set of possible flows, of particular solutions of physical and technical interest; for example, in the case of a small gradient of pressure we have numerically tested, from a set of possible pressure fields, the existence of a solution of Navier-Stokes equations along a given turbulence streamline like that of Fig. (5), which has importance for the prediction of cavitation in pump.