Mathematical symbols

Partial derivative

Jacobian matrix, vector a related to vector

*I* Identity matrix

[]^{-1} Inverse matrix

[]^{t} Transposed matrix

| | Scalar quantity absolute value

δ_{ξ} , δ_{η} Centered difference operator

Δ_{ξ} ,Δ_{η} Forward difference operator

∇_{ξ} ,∇_{η} Backward difference operator

∞ Free-flow properties (non disturbed flow)

*x* Matrix multiplying signal, when not written in the same line

Latin characters

*a* Sound speed (m/s)

Jacobian matrix associated with flux vector E

Jacobian matrix associated with flux vector F

*c* Typical velocity

_{ Cv } Specific heat at constant volume

_{ Cp } Pressure coefficient

*D* Jacobian matrix associated with conservative variables vector

*e* Specific energy (by volume)

_{ ei } Specific internal energy (by volume)

*E* Flux vector at x-direction, (Beam and Warming formulation)

*F*Flux vector at h-direction, (Beam and Warming formulation)

*h* _{0} Stagnation enthalpy

*h* Specific enthalpy

*H* Enthalpy

*J* Jacobian for coordinate transformation

*K*2, *K*4 Non-linear parameters in the artificial dissipation scheme

*L*, *M*, *N* Decomposed “L-U” matrix components

_{ lref } Reference length

*LHS*Left Hand Side for the Euler equations (discretized)

ξ-direction implicit operator

η-direction implicit operator

*O*Order of magnitude

*p*Static pressure (local)

P_{0} Stagnation pressure

*q* Primitive variables vector, [*p*, *u*, *v*, *T*]

q_{d} Dynamic pressure (local)

conservative variables vector (generalized curvilinear coordinates)

*R*Gas universal constant

*Re*Reynolds number (dimensionless)

RHS 2-D Euler equations Right Hand Side (discretized)

_{ Rn } Numerical residue

*R* _{ξ} RHS operator at ξ-direction

R_{η} RHS operator at η-direction

*S* _{ξ} Artificial dissipation at algorithm RHS at ξ-direction

*S* _{η}Artificial dissipation at algorithm RHS at η-direction

*t* Time

*T* Temperature

*T* _{0}Stagnation temperature

*T* _{ξ} Transformation that diagonalizes

Tη Transformation that diagonalizes

u x-direction velocity (Cartesian component)

U ξ-direction velocity (Contravariant component)

v y-direction velocity (Cartesian component)

V η-direction velocity (Contravariant component)

x, y, z Cartesian coordinates

Greek characters

γ Specific heat ratio

εe, εi Artificial dissipation coefficients (linear components)

Δt Time-step

∧ξ Eigenvalues associated with matrix

∧_{η} Eigenvalues associated with matrix

dimensionless time

ξ Curvilinear general coordinates at longitudinal direction

ξ_{ t } , ξ_{ x } , ξ_{ y } Metrical terms at x-direction

η Curvilinear general coordinates at normal direction

η_{ t } , η_{ x } , η_{ y } Metrical terms at η-direction

ρ Fluid density

ρ_{0} Stagnation density

θ Inlet airflow angle

INTRODUCTION

This work is a contribution based on a computational code developed as a numerical tool for simulating flows at low speeds in configurations of complex geometry. Typical examples of situations that can be solved in the context of this work include internal flows in channels and other industrial devices ^{1}.

The primary goal in the long term is to obtain the ability to simulate with great accuracy the flow in such devices, so as to enable the use of numerical simulation in the complex geometry optimization for lowering flow energy use. Numerical code developed to solve general problems, as for internal flows in complex geometries have been shown effective in predicting the physical results, previously available only for experimental basis as laboratories and wind tunnels. There is great interest in having access to a numerical tool able to predict the behavior of the internal flows and problems in aerodynamics area, for aeronautical purposes or not, for understanding the phenomena and flow physical properties that occur inside the channels and devices.

Specific flow fields of interest were modeled by two-dimensional Euler equations. However, the focus of the work presented here is the use of mesh technique for multiple blocks (of meshes) together with a numerical algorithm for all speed regimes (from low speed to supersonic speed) ^{1}.

The results herein considers the existence of an air jet in the flow, air supply real conditions (with specified speed, m/s or mass flow, kg/s). That occurs in an industrial device which can have various purposes, such as movement and separation of grains, sugarcane, among other agricultural products, as well as the possibility of drying by using hot air.

Considerable effort has been made in order to find a numerical methodology able to overcome the lack (at the time the results were obtained) of a single way (keeping the same methodology) to solve efficiently and accurately since incompressible flow at low speeds, up to supersonic flows for which there is a different physical formulation. Considering such lack, ^{2} presented alternatives and comparison capabilities and its respective limitations for various commercial codes available for flow simulations. Multiblock type formulations, applied to structured meshes are also presented by ^{3}.

The numerical method is an extension of an approximate factorization algorithm for compressible flows of ^{4} and ^{5}, in which the working variables are pressure, temperature and speed (Cartesian components). The numerical solution is obtained for steady state conditions with the aid of a convergence acceleration process, derived from ^{6}, without impairing the quality of the result.

In a general way, what is available in non-commercial CFD codes are methods that solve compressible or incompressible flows, individually, not being common to resolving the two conditions by the same method. For incompressible flows, density (ρ, kg/m^{3}) is considered a constant or as a function of temperature *(T,* K), while for compressible flows, density (ρ, kg/m^{3}) varies with both pressure *(p, Pa)* and temperature *(T, K).*

Incompressible flows, for example, use pressure correction in a finite volume formulation, ^{7}, while for compressible flows, classical solution involves the work of ^{(4, 5, 8)} and ^{9} in which the fluid mechanics equations, in the conservative form, are solved simultaneously. The work of ^{(10, 11)} and ^{12} show that incompressible schemes can be extended to the compressible regime speed range, by using pressure-based methods.

GOVERNING EQUATIONS MODEL

The Euler equations, for Fluid Mechanics, written in conservative way, the dependent (or primitive) variable vector (q *= f {p, u, v,* T}) and the conserved variable vector (, generalized curvilinear coordinates) as function of primitive properties, are described respectively by Equations (1), (2) and (3):

The flux vectors *E*(ξdirection) and *F* (η direction), according to formulation given by ^{5}, re-written in terms of the set of dependent variables are presented in Equations (4) and (5):

Equations for *U* and *V,* contravariant velocity components, written as a function of the Cartesian components of speed and metric terms are:

To easier algebraic manipulation, equations represented in vectors , E and F are multiplied by R (J/kg.K). Furthermore, in this nomenclature, is the ratio of specific heats and R is the gas constant.

Equations (1) through (7) are already in dimensionless forms, and the reference values are based on the flow or jet stagnation properties, in the case of internal flows. Further details on the non-dimensional procedure can be found in ^{1}.

The density variation effect is indistinguishable from errors arising from numerical truncation due to the algorithm spatial precision order, thus generating physically dubious results, as also discussed in ^{13}. It is noteworthy that the method presented herein is always a second order in space, and for Mach numbers around 0.1 the flow is typically incompressible. As reported by ^{14}, for Mach numbers lower than 0.3, typical algorithms for compressible flow begin to present the difficulties and errors herein indicated.

Numerical convergence, residues and solution

Numerical derivative for the equation systems given by the theoretical formulation (governing equations) is described in ^{1}. Implicit Euler is the numerical-step scheme in pseudo-time for convergence of typical steady-state solution.

In the previous equation, it is considered , where Δt is the time step and . The time instant in which is to be assessed is *n*+1.

Algorithm final format applied to an all-peed regime, with convergence acceleration, after algebraic manipulation and development of the equation system, can be expressed as:

Eq. (10) will be used for the present work, applying the multi-block methodology solution. The convergence acceleration is achieved because the left side involves the resolution of two tri-diagonals matrixes of diagonal blocks, rather than two full block tri-diagonals as would be required in the original algorithm.

The numerical order for the solution and the manner in which the processes are executed by the computational algorithm, when it is implemented, are briefly presented next. The idea is to show a general view in the solution process, for internal flow in the channels and devices, considering the multi-block procedure indicated in this work.

1°) Reading of computational mesh data;

2°) Reading the individual conditions in the neighborhood, for each mesh at the geometric configuration, recognizing each boundary condition type in neighboring blocks;

3°) Metrics terms and Jacobian matrix calculation;

4°) Initial condition setting in the computational domain, and on its borders;

5°) Right Hand Side, Eq. (10), the explicit calculation in the whole block,

6°) RHS explicit multiplication by, ;

7°) Solving Tri-diagonal matrix (of diagonal blocks) in ξ direction, ;

8°) Explicit multiplication, ;

9°) Solving Tri-diagonal matrix (of diagonal blocks) in η direction, ;

10°) Explicit multiplication, , ;

11°) Update of the current loop solution, qn+1 = qn+Δ qn;

12°) Applying of the specified boundary condition, for each mesh (or block);

13°) Numerical convergence check for the solution in the whole flow field (all blocks or meshes);

14°) If no convergence is reached, returning to step 5;

15°) If convergence is reached, finishes the computational loop.

Decomposed "L-U" matrix components are given by: and . Signs "←" in equations at steps 1 up to 15, represent the allocation of what is on the right side over what is on the left side. The process of "LU" decomposition, known as the Thomas algorithm is used to "reverse" the matrixes.

In (^{1}), it is available the complete methodology developed by the authors, detailing the terms and considerations presented here for the computational design and iterative resolution of the equations (1) to (7), considering the finite difference method. In that other work, the multi-block method developed and used to obtain the results presented here takes the form of exchange of information between blocks exemplified in more detail ways.

Computational results were obtained using the non linear artificial dissipation model since this offers a more careful consideration of the terms (^{15}).

Boundary conditions, constant parameters and relationships

The following values and relationships are used in the equations described, also in the implementation of the proposed computational algorithm and solution of the problems considered in this paper:

Boundary conditions considered include:

✓ Inlet

✓ Outlet

✓ Wall

✓ Free-flow

✓ Inlet jet flow

At the inlet, properties are obtained through the use of one-dimensional (1-D) relationship characteristics, by fixing temperature and speed, and pressure complex geometry conditions to be simulated inside calculated as a result.

At the outlet, the pressure is fixed, and other properties are calculated.

On the wall, zero gradients for pressure and temperature are imposed, and these properties updated as identical to its neighborhood. Speed Cartesian components are calculated through the contravariant components of velocity, U and V, so that the speed in the wall tangent thereto.

For inlet, an air jet positioned at the channel entrance, it is accepted as an initial condition that the whole field has the same properties of stagnation jet properties.

RESULTS AND DISCUSSION

The flow inside ducts, mainly the classical case of internal flow with backward step, was obtained and presented by ^{1} and ^{16} to the methodology used herein. These results are of lower complexity compared to the ones herein presented, but supported the proposed methodology, allowing to try more complex geometry conditions to be simulated inside

industrial devices.

Figure 1 (a) and 1 (b) show ducts and industrial devices of interest in this work and its mesh/grid (computational domain) representation. Other details on how to obtain these meshes are available in ^{1}. A similar device (experimental model scale) was described by ^{17} and used for solid element separation studies in a two-phase flow (liquid/air, solid/crops) through the internal flow of the mixture inside the device indicated.

At the duct inlet side, there is an airjet, featuring a confined internal flow. Numerical results are dimensionless values, based on the stagnation properties ^{1}, given T_{0} = 302.5 *K* (stagnation temperature) and _{ ujet } = 43.80 *m/s* (airjet speed), whose reference value is adimensionalized by reference sound speed. These values are indicated in the Figure legends.

Typical Duct - Airjet at the inlet position

The airflow inside a typical industrial duct (one inlet and 1 outlet), using a simple configuration with 3 (three) multi-block meshes, as shown in Figure 1(a), has its CFD solution in Figure 2 and 3, for the whole flow field behavior and quantities for dimensionless velocity and pressure.

Solutions were obtained for 3 orders of magnitude drop in the numerical residue waste, and the convergence behavior provided in Figure 2 for a typical duct (three blocks).

Results from Figure 3 (a), indicate that the main physical aspect is the solid wall boundary closely above and below the computational domain, that does not allow the velocity vectors to take a dispersive direction. Thus, they end up bouncing on the walls and compressing the flow at the central position. It is also observed, vortex formation just after airjet position, X/L ≈ 4, and two other vortices, much larger, already in X/L ≈ 12.5. The observed behavior is robust and qualitatively correct, consistent, even when reducing the computational domain for X/L =15 (^{(1, 16)}).

In Figure 3 (b), it is possible to notice pressure field floating values in nearby the airjet region, i.e., at the entrance region of the computational domain. It is possible to realize, by the labels (4 digits, non-dimensional values) that such variations are quite small and probably could be eliminated with a redefinition of the range of property values (for example, 2 digits, non-dimensional values).

Figure 3 (c) focuses on the results nearby the entrance region, where the airjet is positioned and where the highest speed values occur. In that region, the overlap of velocity vectors (different directions) and velocity field (intensity given by color scale) allows understanding the flow characteristics inside the typical industrial duct. It also shows, recirculation taking place near the entrance of the computational domain.

Industrial device: 1 inlet (airjet) and 3 outlets

The air flow inside an industrial device (1 inlet and 3 outlets), using a near complex configuration with 5 (five) multi-block meshes, as shown in Figure 1(a), has its CFD solution in Figure 4, for the whole flow field behavior and quantities for dimensionless velocity and pressure.

Computational solution corresponds to numerical convergence by dropping 4 (four) orders of magnitude in the numerical residues (^{1}).

In Figure 4 (a), we find out that in the regions indicated as outlets, it turns out that part of the velocity vectors are "coming inside" due to recirculation occurring in the flow. Those recirculation regions appear in other flow regions within the industrial device studied, and indicate that the geometric shape and/or intensity of the airjet are not able to establish a laminar speed regime, where all the layer speeds between walls can be aligned. This is why numerical models and/or scale model experiments allow indicating which are the existing physical conditions (recirculation, among others) and which one would be most suitable for the purpose for which the device is intended (turbulence is desired or not?). If looking for transporting agricultural products (sugarcane, straw or grain), it is desirable to have laminar flow, while in drying processes and/ or separation is more interesting a turbulent flow behavior, such as seen in Figure 4 (a).

Low pressure regions, i.e., high speed observed in Figure 3 (b), indicate the center in recirculation bubbles, which are formed when the flow meets a geometrical singularity of the corner type. Theoretical formulation considered herein, is valid for the Euler equations, so that typical viscous phenomena can't be "captured" in the steady state solution of the computational code used herein. However, there are many situations where the flow separation point is geometrically determined, by the existence, for example, of solid body discontinuities in its geometry. In these cases, even non-viscous formulation may adequately capture the vortex or recirculation, occurring in real internal flows. This is exactly the physical situation that occurs, as indicated in several regions ("corners") of Figure 4 (b).

The results presented in Figure 4, correspond to a simplified version of a more complex device, evaluated by ^{1}, but that can also be assumed as a real operating system with component parts not assembled for a simplified sugarcane separation. Simulation of geometry part of the full complex configuration intended to observe any eventual changes in the behavior of the velocity vectors when there are more options for the flow to develop. This allows the possibility to further adds new blocks to build the complex configuration part-by-part, which is one of the motivating factors of the methodology developed and presented in this paper.

The flow field solution, by using a 2-D Euler formulation, together multiple blocks to a code with accelerated numerical convergence has been presented for simple configurations (industrial duct) and complex ones (industrial device, partially reproduced for all-speed methodology evaluation). These results proved physically consistent and satisfactory for validation of the developed code ^{1}, comparing solutions obtained previously in the research group and in the technical literature. In fact, computational algorithms applicable all speed regimes were evaluated in other papers from the research group (^{(18, 19, 20)}).

Other all speed algorithms, for incompressible and compressible flow regimes, are widely discussed in plenty of applications in the literature since the 90's up to actual days in plenty of conditions as indicated in ^{(20, 21, 22, 23, 24, 25, 26)} and ^{27}, for both structured and unstructured meshes and upwind or cell-centered schemes. But normally they lack in present an industrial application, and it is also a representative contribution of the present paper.

CONCLUSIONS

Flow behavior is difficult to be properly predicted inside industrial devices with complex geometry, thus requiring adequate engineering design. Otherwise, it can result in non desirable airflow behavior as the recirculation zones in outlet regions, for example. Experimental results for inside parameters in a real scale industrial device were not provided by the manufacturer, only some geometrical parameters.

Thus, this paper investigated the qualitative behavior of the airflow inside industrial ducts and devices. This information can give subsidy to in later stages of the research, get the ability to simulate or perform fluid mechanics experiments with solid material in the flow (sugarcane, straw, soybean and others), aiming for better understanding and deeper comprehension of the physical problem under analysis.

Contribution provided is to present a numerical code developed by the research group (non commercial), which considers a centered method, able to solve airflows in complex configurations with the same flexibility given by unstructured grids, but using multiblock for structured meshes. The implemented numerical scheme is applicable to all-speed flow regime and showed to be robust and able to deal with an increasing part in a complex device by using the developed multiblock methodology.

That allows a great number of applications in practical problems of the Brazilian industries, which has no knowledge of the behavior of airflows inside ducts and devices at low-medium speeds. Applications include transport and separation in the agribusiness industry (including drying processes or heat transfer), like the grain shift, sugarcane and other agricultural products as well as agricultural residues (straw, stem, earth, gravel, etc.).