Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0117 Acta Polytechnica 61(SI):117–121, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague ON THE SENSITIVITY OF THE NONLINEAR TERM IN THE OUTFLOW BOUNDARY CONDITION Tomáš Neustupa∗, Ondřej Winter Czech Technical University in Prague, Faculty of Mechanical Engineering, Karlovo nám. 13, 121 35 Praha 2, Czech Republic ∗ corresponding author: tneu@centrum.cz Abstract. This paper studies the artificial outflow boundary condition for the Navier-Stokes system. This type of condition is widely used and it is therefore very important to study its influence on a numerical solution of the corresponding boundary-value problem. We particularly focus on the role of the coefficient in front of the nonlinear term in the boundary condition on the outflow. The influence of this term is examined numerically, comparing the obtained results in a close neighbourhood of the outflow. The numerical experiment is carried out for a fluid flow through the channel with so called sudden extension. Presented numerical results are obtained by means of the OpenFOAM toolbox. They confirm that the kinetic energy of the flow in the channel can be controlled by means of the proposed boundary condition. Keywords: Navier-Stokes equations, natural outflow boundary condition, finite volume method. 1. Introduction In computational fluid dynamics, the boundaries where the velocity is not known in advance are usually denote by open/artificial boundaries. This situation occurs in mathematical models of many types of fluid flow (e.g. the flow of blood, the flow in various blade machines, etc., see e.g.[1–5]. For instance, the accu- racy of the dynamics of micropolar fluid depends on the boundary conditions [6–8]. In these cases, the ve- locity profile is rarely available in advance on the whole boundary of the flow field, the pressure is available in some special cases when it is measured or computed with the aid of a reduced model. Furthermore, the necessity of setting an appropriate boundary condi- tion on an artificial part of the boundary becomes also important when the computational domain is obtained by truncating the length of the domain in order to reduce the computational cost. One of the boundary conditions addressing the prob- lem of the open boundary is the so called “do–nothing” boundary condition used e.g.by Heywood, Rannacher and Turek in [9] (see also the so called natural bound- ary condition [10]), i.e., −ν ∂u ∂n + p n = 0, (1) where u = (u,v) denotes the velocity of the fluid, p is the kinematic pressure, i.e., the pressure divided by constant fluid density, ν is the kinematic viscosity (which is assumed to be a positive constant) and n is the unit outer normal vector to the boundary of the considered domain. However, the condition (1) does not enable one to control the amount of kinetic energy in the domain if a backward flow appears on the “open boundary” (which is the part, where the velocity profile is not given). Bruneau and Fabri proposed in [11] the boundary condition −ν ∂u ∂n + p n − 1 2 (u · n)− u = 0, (2) as a natural modification of (1). This extension comes naturally when the symmetric part of the convective term is integrated by parts. The superscript “−” de- notes the negative part (i.e. (u · n)− = −u · n if u · n < 0, otherwise (u · n)− = 0). Thus, the in- equality (u · n)− > 0 is satisfied only in the case of a “backward flow” on the open boundary. This modification enables one to prove the existence of a weak solution, but only if the inflow velocity profile is bounded with respect to the viscosity see e.g. [2]. The same condition is also used in [12] on a part of the boundary. Here, the authors prove the existence of a weak solution, under the stronger assumption, i.e. that the inflow velocity is zero. Neustupa in [13] proposed a modification of (2), i.e., −ν ∂u ∂n + p n − 1 + ξ 2 (u · n)− u = h, (3) where ξ is a constant dimensionless parameter and h is an arbitrary vector function. This boundary condi- tion, in comparison to (2), contains a correction in the nonlinear part. The correction enables the author to derive necessary a priori estimates of a solution in the case of an arbitrarily large inflow. These results show that the coefficient in front of the nonlinear part of the used boundary condition plays an important role in theoretical considerations, particularly in the existen- tial theory. If it is less then 12 (which corresponds to ξ < 0 in condition (3)) then the existence of the weak solution is an open problem. If it is exactly 12 (which 117 https://doi.org/10.14311/AP.2021.61.0117 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en Tomáš Neustupa, Ondřej Winter Acta Polytechnica D i D e li le Ω Γi Γo Γw+ Γw− x y Figure 1. Sketch of the computational domain Ω. corresponds to ξ = 0 in condition (3)) then the exis- tence of the weak solution can be proved, assuming a certain restriction of the size of the inflow velocity. If the coefficient is greater than 12 (which corresponds to ξ > 0 in condition (3)) then the existence of the weak solution can be proved for arbitrary large inflow, see [13]. In the field of numerical simulations, the problem with the original “do–nothing” boundary condition is well known and the modification (2) is widely used. Our goal is to test numerically the behaviour of the flow in the neighbourhood of the outflow part of the boundary, in dependence on the dimensionless param- eter ξ in front of the nonlinear part of the boundary condition. We are especially interested in compari- son of numerical results, obtained in the three cases, i.e. when the dimensionless parameter is less that, equal to, or greater than zero. 2. Mathematical Model The stationary flow of a viscous incompressible Newto- nian fluid is described by the Navier–Stokes equation and the equation of continuity, i.e., (u ·∇)u + ∇p = ν∆u + f, in Ω, (4) div u = 0, in Ω. (5) Here, f is a specific volume force. Since the parts of the boundary of Ω are of different types, we impose different boundary conditions. Concretely, we assume that the boundary of Ω consists of four curves: ∂Ω = Γi ∪ Γo ∪ Γw+ ∪ Γw−, see Figure 1. The curve Γi represents the inlet (i.e. the part of boundary where the fluid enters the considered domain Ω) and Γo is the outlet (where the fluid leaves Ω). The curves Γw+, Γw− are the non–permeable fixed walls of the channel. We assume that the whole boundary ∂Ω is Lipschitz–continuous. We consider Dirichlet’s boundary conditions on Γi, Γw+ and Γw−, i.e., u = g on Γi, (6) u = 0 on Γw+. (7) u = 0 on Γw−. (8) The curve Γo represents an artificially chosen part of the boundary, and the velocity profile on Γo is therefore not known in advance. We consider this concrete “artificial” boundary condition on Γo: −ν ∂u ∂n + p n − 1 + ξ 2 (u · n)− u = 0 on Γo. (9) 3. Numerical Approximation The numerical solution of the governing system is based on a collocated finite-volume method imple- mented in the freely available CFD toolbox Open- FOAM [14]. The solver uses segregated approach and pressure–velocity coupling is done with aid of the Semi-Implicit Method for Pressure-Linked Equa- tions (SIMPLE) algorithm, see e.g. [15]. The con- vective term appearing in Navier–Stokes equation is discretized using the limited piece-wise linear recon- struction and the viscous term is approximated using a central scheme, see [16] or [17] for more details on the spatial discretization. 3.1. Boundary Conditions The boundary conditions are implemented in the usual way, used in the finite volume framework, see e.g. [16]. Since the SIMPLE algorithm uses the elliptic equation for pressure we need to prescribe boundary condition for the pressure on the whole boundary of Ω. The numerical implementation of the boundary conditions given by (6, 7, 8, and 9) is realized as follows: • Γi: The velocity profile u = (u(y), 0) is prescribed, i.e., the fully developed parabolic profile is given as u(y) = 3U 2 − 6U D2i y2 , (10) where Di is inlet diameter and U denotes the mean value of the magnitude of the velocity at the inlet. The homogeneous Neumann boundary condition for pressure is used; • Γw+ ∪ Γw−: no-slip boundary condition, i.e., u = 0 and the homogeneous Neumann boundary condition for pressure is used; • Γo: the artificial boundary condition (3) is imple- mented in the finite volume framework in the fol- lowing way, i.e., ∂u ∂n ∣∣∣ b = 1 ν [ (pb −p0)n − 1 + ξ 2 (up · n)−up ] , (11) where p0 is the referential value of the pressure, constant on Γo. Indices b and p denote the value on the boundary and the internal value at the nearest degree of freedom placed along the normal direc- tion to the boundary, respectively. The pressure is realized so that� Γo pdΓ − � Γo p0 dΓ = 0. (12) This means that we are prescribing only one piece of information on the pressure on the whole line segment Γo. 118 vol. 61 Special Issue/2021 On the Sensitivity of the Nonlinear Term. . . −2 0 2 0 0.5 1 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 −2 0 2 0 0.5 1 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 Figure 2. u along the boundary Γo, Re = 10, ξ < 0,ξ = 0. Long - result obtained on prolonged domain (le = 50Di). 4. Numerical Results The influence of the coefficient ξ is studied in case of a flow through the channel with the so called sudden extension, see Figure 1. The computations were done for ξ = −0.6,−0.5,−0.3̄, 0, 0.3̄, 0.5, 0.6 corresponding to 15, 1 4, 1 3, 1 2, 2 3, 3 4, 4 5 of the coefficient in front of the nonlinear term in (3). Furthermore, the computa- tions were done for three different Reynolds numbers, Re = UDi/ν, Re = 10, 100, 500. For an optimal occurrence of the backward flow in the neighbour- hood of the outlet, we used different ratios of the sudden extension De/Di (outlet diameter / inlet di- ameter) for each Reynolds number, namely 4, 2, 1.5 for Re = 10, 100, 500, respectively. The following data were considered: Di = 1 m, U = 1.5 m/s, li = 1.5 m, le = 0.5 m. Figures 2 and 3 show profiles of the velocity u on Γo for Re = 10, ξ ≤ 0 and ξ ≥ 0, respectively. There are no visible differences in dependence on the varying ξ. Figure 4 shows details of the contours of u for Re = 10 for ξ = −0.6, 0, 0.6. One can see that ξ = 0 and ξ = −0.6 gives almost the same results, but there is a significant difference (shift) between the contours obtained with ξ = 0 and ξ = 0.6. Figure 5 shows profiles of u on Γo for Re = 100, ξ ≤ 0. Small differences can be observed for different ξ. However, for Re = 100 we were not able to obtain any solution for ξ > 0 due to lack of convergence. More cases were computed (not published in this work) with different extension ratios for ξ ∈ 〈0, 0.3̄〉 and we were not able to find a general sharp borderline for value of ξ. The lack of convergence seems to be dependent on the magnitude of the velocity occurring in the region of Γo where a backward flow occurs. Figure 6 shows a detail of the contours of u for Re = 100 −2 0 2 0 0.5 1 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 −2 0 2 0 0.5 1 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 Figure 3. u along the boundary Γo, Re = 10, ξ > 0,ξ = 0. Long - result obtained on prolonged domain (le = 50Di). Figure 4. Detail of the contours of u, in region D, see Figure 1, Re = 10. Black, red, blue, and green lines indicate results for prolonged domain (le = 50Di), ξ = 0, ξ = −0.6 and ξ = 0.6, respectively. for ξ = 0,−0.6. One can see that there is a good correspondence between the contours for ξ = 0 in a “longer” domain and the results for ξ = −0.6 are slightly shifted. This can be explained by the fact that as ξ → 1, the boundary condition (3) converges to the do-nothing boundary condition (1) which does not take into account any backward flow on Γo. Figure 7 shows profiles of u at Γo for Re = 500, ξ ≤ 0. No virtual differences can be observed for different ξ. Figure 8 shows a detail of the contours of u for Re = 500 for ξ = 0,−0.6. Similar conclusion can be made for ξ > 0 as for the case of Re = 100. 119 Tomáš Neustupa, Ondřej Winter Acta Polytechnica −1 0 1 0 0.5 1 1.5 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 −1 0 1 0 0.5 1 1.5 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 Figure 5. u along the boundary Γo, Re = 100, ξ < 0,ξ = 0. Long - result obtained on prolonged domain (le = 50Di). Figure 6. Detail of the contours of u in region D, see Figure 1, Re = 100. Black, red, and blue lines indicate results for prolonged domain (le = 50Di), ξ = 0, ξ = −0.6, respectively. 5. Conclusion Numerical investigation of the boundary condition (3), allowing to control the amount of kinetic energy in the domain was done. The boundary condition was tested for different magnitudes of the inflow velocity with respect to the viscosity, namely Re = 10, 100, 500. From the obtained results one can conclude that if the inflow velocity is sufficiently small then ξ can be chosen so that ξ > 0 (ξ ∈ 〈−0.6, 0.6〉 in our simula- tions). With an increasing inflow velocity, lack of the convergence for ξ > 0 can be observed. From other numerical results (not presented in this work), it is possible to conclude that the convergence for ξ > 0 strongly depends on the magnitude of the possible reverse velocity on the outflow and also on the size of the backward flow area. Studying the dependence of the backward flow area, as a subset of Γo, on multiple parameters, e.g. inlet velocity, Reynolds number, and extension ratio we were not able to find any sharp general border for convergence criteria on coefficient ξ. However, one can observe that the region in which −0.5 0 0.5 0 0.5 1 1.5 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 −0.5 0 0.5 0 0.5 1 1.5 y Long ξ = −0.3̄ ξ = −0.5 ξ = −0.6 ξ = 0 Figure 7. u along the boundary Γo, Re = 500, ξ < 0,ξ = 0. Long - result obtained on prolonged domain (le = 50Di). Figure 8. Detail of the contours of u in region D, see Figure 1, Re = 500. Black, red, and blue lines indicate results for prolonged domain (le = 50Di), ξ = 0, ξ = −0.6, respectively. we are able to find a numerical solution is approxi- mately Re ∈ (0, 50〉 for ξ ∈ 〈−0.6, 0.6〉. Hence, the used numerical method converges for small inlet data both for ξ ∈ 〈−0.6, 0.6〉. For larger inlet data, the method converges only for ξ ∈ 〈−0.6, 0〉, where ξ = 0 corresponds to boundary condition introduced in [11]. The used numerical method does not confirm the the- oretical conclusion presented in [13], i.e. that for ξ > 0 the existence of a weak solution can be proved for an arbitrary large inflow, in the sense that the numerical method does not converge in the situations described above. There are still relatively many open problems, con- nected with the boundary condition (3), both in the field of qualitative analysis and numerical computa- tions. For example, the condition (3) can be reformu- lated in terms of the pressure (see e.g. [18]) which can be possibly more suitable for finite volume framework and SIMPLE algorithm. Furthermore, it also seems desirable to test another approach, e.g. based on the finite element method, due to its closer relation with the weak formulation. 120 vol. 61 Special Issue/2021 On the Sensitivity of the Nonlinear Term. . . Acknowledgements The research was supported by European Regional De- velopment Fund-Project “Center for Advanced Applied Science” No. CZ.02.1.01/0.0/0.0/16_019/0000778. Authors acknowledge support from the EU Opera- tional Programme Research, Development and Education, and from the Center of Advanced Aerospace Technology (CZ.02.1.01/0.0/0.0/16_019/0000826), Faculty of Mechan- ical Engineering, Czech Technical University in Prague and by the Grant No. 19-04477S of the Czech Science Foundation. References [1] M. Feistauer, T. Neustupa. On Non-Stationary Viscous Incompressible Flow Through a Cascade of Profiles. Mathematical Methods in the Applied Sciences 29(16):1907–1941, 2006. doi:10.1002/mma.755. [2] M. Feistauer, T. Neustupa. On Some Aspects of Analysis of Incompressible Flow through Cascades of Profiles. In Operator Theoretical Methods and Applications to Mathematical Physics, vol. 147, pp. 257–276. Birkhäuser Basel, 2004. doi:10.1007/978-3-0348-7926-2_29. [3] M. Feistauer. Mathematical Methods in Fluid Dynamics. Longman Scientific & Technical, 1993. [4] S.I. Abdelsalam. M.M Bhatti. The Study of Non-Newtonian Nanofluid with Hall and Ion Slip Effects on Peristaltically Induced Motion in a Non-uniform Channel. RSC Advances 8(15):7904–7915, 2018. doi:10.1039/C7RA13188G. [5] Y. Abd elmaboud, S.I. Abdelsalam, K.S. Mekheimer Couple Stress Fluid Flow in a Rotating Channel with Peristalsis. Journal of Hydrodynamics 30:307–316, 2018. doi:10.1007/s42241-018-0037-2. [6] S. S. Motsa, I. L. Animasaun Bivariate Spectral Quasi-linearisation Exploration of Heat Transfer in the Boundary Layer Flow of Micropolar Fluid with Strongly Concentrated Particles over a Surface at Absolute Zero due to Impulsive. International Journal of Computing Science and Mathematics 9:455–473, 2018. doi:10.1504/IJCSM.2018.10016504. [7] O. K. Koriko, I. L. Animasaun New Similarity Solution of Micropolar Fluid Flow Problem over an Uhspr in the Presence of Quartic Kind of Autocatalytic Chemical Reaction Frontiers in Heat and Mass Transfer 8, 2017. doi:10.5098/hmt.8.26. [8] O. K. Koriko, A.J. Omowaye, I. L. Animasaun, M. E. Bamisaye Melting Heat Transfer and Induced-Magnetic Field Effects on the Micropolar Fluid Flow towards Stagnation Point: Boundary Layer Analysis International Journal of Engineering Research in Africa 29:10–20, 2017. doi: 10.4028/www.scientific.net/JERA.29.10. [9] J. G. Heywood, R. Rannacher, S. Turek. Artificial Boundaries and Flux and Pressure Conditions for the Incompressible Navier-Stokes Equations. International Journal for Numerical Methods in Fluids 22(5):325–352, 1996. doi:10.1002/(sici)1097- 0363(19960315)22:5<325::aid-fld307>3.0.co;2-y. [10] R. Glowinski. Numerical Methods for Nonlinear Variational Problems. Springer, 1984. [11] C.-H. Bruneau, P. Fabrie. New Efficient Boundary Conditions for Incompressible Navier-Stokes Equations: A Well-Posedness Result. Mathematical Modelling and Numerical Analysis 30(7):815–840, 1996. doi:10.1051/m2an/1996300708151. [12] M. Braack, P. B. Mucha. Directional Do-Nothing Condition for the Navier-Stokes Equations. Journal of Computational Mathematics 32(5):507–521, 2014. doi:10.4208/jcm.1405-m4347. [13] T. Neustupa. A Steady Flow through a Plane Cascade of Profiles with an Arbitrarily Large Inflow: The Mathematical Model, Existence of a Weak Solution. Applied Mathematics and Computation 272:687–691, 2016. doi:10.1016/j.amc.2015.05.066. [14] H. G. Weller, G. Tabor, H. Jasak, C. Fureby. A Tensorial Approach to Computational Continuum Mechanics Using Object-Oriented Techniques. Computers in Physics 12(6):620–631, 1998. doi:10.1063/1.168744. [15] J. H. Ferziger, M. Perić. Computational Methods for Fluid Dynamics. Springer, 3rd edn., 2002. [16] F. Moukalled, L. Mangani, M. Darwish. The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM and Matlab. Springer, 2016. [17] O. Winter, P. Sváček. On Numerical Simulation of Flexibly Supported Airfoil in Interaction with Incompressible Fluid Flow using Laminar-Turbulence Transition Model. Computers & Mathematics with Applications 2020. doi:10.1016/j.camwa.2019.12.022. [18] T. Bodnár, P. Fraunié. Artificial Far-Field Pressure Boundary Conditions for Wall-Bounded Stratified Flows. In Topical Problems of Fluid Mechanics 2018. Institute of Thermomechanics, AS CR, v.v.i., 2018. doi:10.14311/tpfm.2018.002. 121 http://dx.doi.org/10.1002/mma.755 http://dx.doi.org/10.1007/978-3-0348-7926-2_29 http://dx.doi.org/10.1039/C7RA13188G http://dx.doi.org/10.1007/s42241-018-0037-2 http://dx.doi.org/10.1504/IJCSM.2018.10016504 http://dx.doi.org/10.5098/hmt.8.26 http://dx.doi.org/ 10.4028/www.scientific.net/JERA.29.10 http://dx.doi.org/ 10.4028/www.scientific.net/JERA.29.10 http://dx.doi.org/10.1002/(sici)1097-0363(19960315)22:5<325::aid-fld307>3.0.co;2-y http://dx.doi.org/10.1002/(sici)1097-0363(19960315)22:5<325::aid-fld307>3.0.co;2-y http://dx.doi.org/10.1051/m2an/1996300708151 http://dx.doi.org/10.4208/jcm.1405-m4347 http://dx.doi.org/10.1016/j.amc.2015.05.066 http://dx.doi.org/10.1063/1.168744 http://dx.doi.org/10.1016/j.camwa.2019.12.022 http://dx.doi.org/10.14311/tpfm.2018.002 Acta Polytechnica 61(SI):117–121, 2021 1 Introduction 2 Mathematical Model 3 Numerical Approximation 3.1 Boundary Conditions 4 Numerical Results 5 Conclusion Acknowledgements References