Numerical Simulation of Magnetohydrodynamic Flow of Power-law Fluid in Channel – An Investigation Under Three Different Flow Conditions

Numerical Simulation of Magnetohydrodynamic Flow of Power-law Fluid in Channel – An Investigation Under Three Different Flow Conditions

Hardik Manjeshwar, Ranjith Maniyeri

Biophysics Laboratory, Department of Mechanical Engineering, National Institute of Technology Karnataka,

Surathkal, Karnataka.

Abstract – Characterized by their non-Newtonian rheology, power-law fluids arise in numerous engineering and biological applications where fluid viscosity varies with shear rate. Among various flow configurations, pressure-driven, shear-driven and oscillatory flows are of particular importance due to their relevance in confined geometries. In this work, pressure-driven, shear-driven, and oscillatory flows of power-law fluid in straight channel are simulated in two-dimensional computational model developed using MATLAB programming language and it is further extended to analyze the impact of external magnetic field oriented perpendicular to flow direction. The governing continuity and Navier-Stokes equations are discretized using finite volume method and solved with the help of semi-implicit fractional step algorithm under prescribed inlet velocity and pressure boundary conditions. The results reveal that change in Reynolds number, power-law index, and Hartmann number significantly affect the flow characteristics. Pressure-driven flow shows a reduction in boundary layer thickness with an increase in Reynolds number from 10 to 50 and 100 before fully developed flow condition is reached. Shear thickening fluid shows thicker boundary layer as compared to Newtonian and shear thinning fluids. Application of magnetic field inhibits fluid motion, with fluid attaining a more plug-like velocity profile. Shear-driven flow with top wall moving shows no change in velocity with a change in Reynolds number and power-law index, with all fluids displaying a linear velocity profile. Shear thickening fluid shows highest velocity in presence of magnetic field due to its resistive nature. In case of oscillatory flow, fluctuations are introduced at every time instance in the flow field. Higher fluctuation levels are observed at low Reynolds number, while they diminish with increasing Reynolds number with shear thinning fluids responding actively to the imposed disturbance as compared to other fluids. Under magnetic field, fluctuations are negligible and the flow stabilizes.

Keywords – finite volume method, magnetohydrodynamic flow, oscillatory flow, power-law fluid, shear-driven flow.

1. INTRODUCTION

   Pressure-driven, shear-driven and oscillatory flows of power-law fluid plays a vital role in modelling non-Newtonian transport phenomena in various engineering and biological systems. In pressure-driven flow, pressure gradient is the driving mechanism of fluid motion. In shear-driven flow, fluid motion is induced by the movement of channel boundaries, such as in Couette flow, where velocity distribution depends nonlinearly on the shear rate due to fluid shear thinning or shear thickening

   nature. Oscillatory flow, on the other hand, involves time- periodic motion typically driven by sinusoidal velocity variations, leading to complex dynamic behaviour. These flow types exhibit distinct velocity and stress profiles due to power- law model, offering valuable insights for applications in chemical mixing devices, turbomachinery, rotating-tube heat exchangers, coating processes, viscometers, oil and gas drilling, control of blood flow during surgical operations, manufacturing and processing of foods and paper, modelling of respiratory functions in lungs, modelling of chemical/blood dispensing in biochemistry/clinical labs, cooling of electronic device, etc. Understanding these flows enables more accurate prediction and control of non-Newtonian behaviour in confined geometries.

   Various advanced numerical techniques are used in prior investigations to evaluate flow characteristics in scenarios where analytical solutions are impractical and experimental setups are costly or challenging to implement. Chaudhuri and Sahoo [1] used finite difference method (FDM) with an explicit time marching scheme to examine non-Newtonian fluid flow in rectangular channel under both pressure-driven and shear-driven conditions. They looked at how different fluid parameters affect velocity and shear stress in pseudoplastic and dilatant fluids. Garang et al. [2] investigated the time dependent flow of non- Newtonian fluid in a channel with varying magnetic field keeping the lower wall fixed and upper wall in motion at constant velocity. FDM with forward time central space (FTCS) scheme was used to solve the governing equations. The variation of temperature and velocity of fluid were analyzed by changing various fluid parameters. Ershadnia et al. [3] clarified how shear thinning effect affects the dynamic behavior of pressure loss (ratio) in rotating yield-power-law fluids and how it interacts with inertial forces. Using a coupled simulated annealing (CSA) approach to improve the least square support vector machine (LSSVM), they created a unique and computationally effective machine learning based prediction model. Chaudhuri et al. [4] investigated heat transfer in a third-grade fluid that is non- Newtonian, which flows between two parallel plates due to pressure gradient and shear forces. They applied semi-analytical approach to solve the non-linear governing equations using least-square method (LSM) and examined how Brinkman number and heat flux ratio affect Nusselt number and temperature of fluid. Liu and Liu [5] conducted a study on linear stability of power-law fluid in plane Couette flow by employing eigenvalue analysis through Chebyshev collocation method,

   energy method, and non-modal stability theory. According to their findings, viscosity gradients generated by shear thinning enhanced the stability of general plane shear flows, whereas viscosity disturbances contributed to destabilization. Makinde and Franks [6] used FDM with Runge-Kutta Fehlberg integration solver to study the collective impact of time- dependent effects, variable properties, and magnetic fields on reactive third-grade hydromagnetic Couette flow and associated heat transfer. In their study, Rath and Mahapatra [7] conducted numerical investigation into pulsatile flow of Oldroyd-B viscoelastic fluid within a parallel plate channel utilizing open- source RheoTool toolbox integrated with OpenFOAM framework. They employed fully implicit finite volume method (FVM) to solve the governing equations, which was based on time-marching pressure-correction algorithm. For coupling pressure and velocity, semi-implicit method for pressure-linked equations-consistent (SIMPLEC) algorithm was used. Their study highlighted the interplay between pulsation parameters and fluids complex rheology on flow characteristics. It was found that as fluids viscoelastic properties increased, there was a notable rise in axial velocity and polymeric shear stress near the wall, whereas axial velocity decreased with higher pulsation amplitude or frequency of time averaged pulsatile flow. Bukhari et al. [8] statistically studied the behavior of thermally generated Casson fluid undergoing pulsatile flow in a constricted channel under magnetic field effects. FDM was utilized for numerical treatment of governing equations. This study looked at how axial velocity and dimensionless temperature profiles were affected by magnetic field parameter, pulsatile flow parameter, Casson fluid parameter, Darcy parameter, Prandtl number, and thermal radiation parameter.

   While non-Newtonian fluid flows have been extensively studied under pressure gradients, magnetic fields and complex geometries, the influence of steady pressure-driven, shear- driven and oscillatory motions on power-law fluids in straight channel remains insufficiently explored. Existing research is either focused on Newtonian fluids or more complex rheological models, overlooking the simpler yet widely applicable power- law model. This study aims to address this gap by investigating the flow characteristics of power-law fluids in straight channel under pressure-driven, shear-driven and oscillatory effects for three power-law index (z) values 0.75, 1.00 and 1.50 at three different Reynolds number (Re) values 10, 50 and 100. Continuity and Navier-Stokes equations governing the fluid flow are solved on staggered Cartesian grid system using FVM. We develop a code using MATLAB programming language to perform numerical simulations to demonstrate the flow. The velocity variations of each case are captured and analyzed to understand the fluid flow behaviour.

2. METHODOLOGY

   The non-Newtonian nature of fluid is modelled using Ostwald-de Waele power-law formulation, wherein the dependence of shear stress on shear strain rate is mathematically formulated as follows,

   Numerical model is developed based on the assumptions of unsteady incompressible laminar flow.

   1. Pressure-driven flow

      The dimensionless schematic of computational domain of pressure-driven straight channel configuration with length L and height H and associated boundary conditions are presented in Fig. 1.

      At the inlet of this geometry, magnitude of streamwise velocity (u) is taken as unity and transverse velocity (v) is kept zero. Neumann boundary condition for pressure (p) is assumed at the inlet. The velocity of fluid relative to walls is taken to be zero and flow continuity is maintained at the outlet.

   2. Shear-driven flow

      For shear-driven flow configuration, same computational domain is considered; however, flow is generated by motion of one of the channel walls instead of an imposed inlet velocity condition. The lower wall is maintained stationary, while upper wall moves with unit dimensionless velocity. To maintain the continuity of flow, periodic boundary conditions are used at inlet and outlet as shown in Fig. 2.

   3. Oscillatory flow

      For oscillatory flow case, same straight channel geometry is retained, but the driving mechanism is time dependent. To replicate the behaviour of sinusoidal velocity input sin(),

      an external sinusoidal pressure gradient of 12 sin is given

      with zero velocity at the walls and periodic boundary conditions

      at inlet and outlet as shown in Fig. 3. Here, A is amplitude of oscillation, is angular frequency, and t is time.

      Fig. 1. Pressure-driven straight channel geometry with boundary conditions.

      Fig. 2. Shear-driven straight channel geometry with boundary conditions.

      = () (1)

      where is shear stress of fluid, is shear strain rate of fluid, and is consistency factor. Accordingly, apparent viscosity () of power-law fluid can be expressed as,

      (2)

      1

      = ()

      Fig. 3. Oscillatory flow straight channel geometry with boundary conditions.

   4. Governing equations: non-dimensional form

   To model the MHD flow problem, continuity and Navier- Stokes equations are formulated in non-dimensional vector form as follows,

   = 0 (3)

   + ( ) = + 1 { [ + ()]} 2

   (4)

   where u is vector form of velocity, Re is Reynolds number, Ha

   is Hartmann number, and

   is a dimensionless term given by,

   1

   = {2 [(

   )2

   + (

   2

   ) ] + (

   +

   )2 } 2

   (5)

   The equations (3) and (4) are solved numerically over the entire computational domain.

3. RESULTS AND DISCUSSION

   The continuity and momentum equations as mentioned in section II are solved in a two-dimensional (2D) computational domain of 5 × 1. A uniform Eulerian grid size of 201 × 41 in x and y directions respectively is used for simulation with a time step of 0.001. A finite volume framework is adopted for discretization of governing equations with staggered grid configuration, and pressure-velocity coupling is treated through semi-implicit fractional step algorithm [9]. Incomplete Cholesky conjugate gradient (ICCG) method is also employed in the simulation [10]. The resulting discretized equations are solved numerically by constructing a MATLAB code. Numerical simulations are conducted for Reynolds number, Re = 10, 50, and 100, and power-law index, z = 1.00, 0.75, and 1.50.

   1. Validation

      Initially, the code is validated for plane Poiseuille flow of fluid governed by power-law rheology. Results of plane Poiseuille flow are compared with the numerical results obtained by [11] using lattice-Boltzmann method (LBM). The physical model describes pressure-driven flow between two parallel plates along streamwise direction. Under steady state conditions, flow becomes unidirectional along transverse direction. Simulations are carried out for power-law indices of 0.50, 1.00, and 1.50, with no-slip conditions at upper and lower walls and periodic boundary conditions along x-direction. It can be observed from Fig. 4 that velocity profiles match closely with the profiles obtained by [11], thus validating the developed code.

   2. Pressure-driven flow

      Fig. 5 shows the u-velocity contour and vorticity contour plots for power-law fluid flow in straight channel with selected values of z and Re. It can be observed from Fig. 5 that for Newtonian fluid (z = 1.00) at low Reynolds number (Re = 10) viscous forces are dominant, and the fluid experiences significant resistance to motion near the wall leading to thicker boundary layer. At high Reynolds number (Re = 50 and 100) inertial forces dominate over viscosity causing the fluid to accelerate more rapidly within boundary layer which leads to steeper velocity gradients near the wall and thinner boundary layer before it reaches fully developed flow condition. The maximum velocity in case of Newtonian fluid reaches 1.50 for fully developed flow as expected in plain Poiseuille flow. Same effect can be observed in case of shear thinning (z = 0.75) and shear thickening fluids (z = 1.50) as shown. Shear thinning fluid

      Fig. 4. Velocity profile comparison of Poiseuille flow for non-Newtonian power-law fluid with results obtained by [11].

      exhibits decrease in viscosity as shear rate increases. At the wall, fluid experiences a very high shear rate due to velocity gradient. Because the viscosity decreases with increased shear rate, fluid can flow more easily near the wall, reducing resistance to motion and reduction in boundary layer thickness. Shear thinning fluid shows maximum velocity of 1.43 which is lower than Newtonian fluid. Shear thickening fluid on the other hand shows increase in viscosity with shear rate. This increased viscosity slows down the flow near wall causing a more gradual transition from no-slip condition to free stream velocity. As a result, momentum diffusion is slower leading to thicker boundary layer. Shear thickening fluid gives maximum velocity of 1.60 which is higher than Newtonian fluid.

      Fig. 5. uvelocity contour plots for pressure-driven power-law fluid flow through straight channel for z = 1.00, 0.75, 1.50; Re = 10, 50, 100; and Ha = 0.

      The distribution of vorticity field can be understood from Fig. 6. For Newtonian fluid, at low Re, vorticity contour shows positive values along the upper wall and negative values along the lower wall due to onset of rotational flow in boundary layer region. As Re is increased, due to inertia, stronger positive and negative vorticity values are observed that extendfurther downstream. In case of shear thinning fluid, due to reduced viscosity near high shear regions, vorticity is more confined towards the walls and contours appear sharper. Shear thickening

      Fig. 6. Vorticity contour plots for pressure-driven power-law fluid flow through straight channel for z = 1.00, 0.75, 1.50; Re = 10, 50, 100; and Ha = 0.

      fluid with an increase in viscosity at high shear region allows smoother transition of vorticity from wall towards channel centerline and has more spread as compared to Newtonian fluid.

      Now, if an external transverse magnetic field is applied across the channel, Lorentz force is induced in flow field that acts counter to fluid motion due to MHD effect. This force hinders the movement of fluid, thereby reducing the magnitude of streamwise velocity. The u-velocity contours and vorticity contours are shown in Fig. 7 and Fig. 8 respectively showcasing this effect. Due to the presence of magnetic field velocity profile attains more plug-like structure. So, transition from core velocity to zero velocity at wall has to occur in a thin region near walls. This region is called Hartmann boundary layer. Strengthening magnetic field will reduce the thickness of this boundary layer.

      Fig. 7. uvelocity contour plots for pressure-driven power-law fluid flow through straight channel for z = 1.00, 0.75, 1.50; Re = 10, 50, 100; and Ha =

      20.

      Fig. 8. Vorticity contour plots for pressure-driven power-law fluid flow through straight channel for z = 1.00, 0.75, 1.50; Re = 10, 50, 100; and Ha =

      20.

   3. Shear-driven flow

      Fig. 9 shows the u-velocity profiles for top wall driven shear flow of power law-fluid in the channel for power-law indices 0.75, 1.00 and 1.50 and Reynolds number 10, 50, and 100 at Ha

      = 0. It can be observed that the velocity distribution is independent of power-law index and Reynolds number under fully developed condition as indicated by straight overlapping velocity profiles for all fluids. This is due to the fact that in shear- driven flow, the flow is governed solely by imposed wall motion in the absence of any external pressure gradient. Thus, shear stress remains constant across the channel width for Newtonian as well as non-Newtonian fluids. This makes the shear strain rate

      Fig. 9. uvelocity profile comparison for shear-driven flow of power-law fluid with z = 1.00, 0.75, 1.50 at Re = 10, 50, 100 and Ha = 0.

      Fig. 10. uvelocity profile comparison for shear-driven flow of power-law fluid with z = 1.00, 0.75, 1.50 at Re = 10, 50, 100 and Ha = 20.

      constant as well. As a result, the velocity gradient does not vary in transverse direction, leading to a linear velocity distribution between the two walls.

      Fig. 10 shows the u-velocity profiles for top wall driven shear flow in presence of magnetic field applied perpendicular to streamwise direction with Ha = 20 for power-law indices 0.75, 1.00 and 1.50 and Reynolds number 10, 50, and 100. Due to Lorentz force, there is a deviation in the magnitude of velocity of all three fluids from straight line velocity profile. At a given y-coordinate, shear thickening fluid has velocity greater than that of Newtonian and shear thinning fluid because of higher momentum diffusion due to increased viscosity making the magnetic field less effective.

   4. Oscillatory flow

   To understand the effect of sinusoidal velocity input on flow characteristics of power-law fluid, it is important to observe how the imposed velocity changes with time. Fig. 11 shows the sinusoidal waveform with different time instants, namely ,

   + , + 2, + 3, and + 4 where dt is the time step. These time levels help illustrate the gradual variation of input velocity during one cycle. The sinusoidal amplitude A is taken as unity with an angular frequency = 2. The velocity profiles of power-law fluid are compared at each time instant and the effect of Reynolds number and Hartmann number are illustrated further.

   The comparison of u-velocity profiles for oscillatory flow of power-law fluid with power-law indices 1.00, 0.75, and 1.50 at Reynolds number values of 10, 50, and 100 are shown in Fig. 12, Fig. 13, and Fig. 14 for different time instances of sinusoidal velocity input. At Re = 10, noticeable differences are observed

   Fig. 11. Sinusoidal velocity input waveform.

   among the velocity profiles of three fluids. The velocity of fluid completely reverses as it reaches the negative peak from first time instant to second time instant + as shown in Fig. 12 (a) and Fig. 12 (b). As time advances from time instant + to + 2 and + 3 to reach the positive peak, the fluid velocity tends to attain positive values as shown in Fig. 12 (c) and Fig. 12 (d). At time instant + 4, fluid again attains the velocity profile obtained at time instant and

   (a)

   (b)

   (c)

   (d)

   Fig. 12. uvelocity profile comparison for oscillatory flow of power-law fluid with z = 1.00, 0.75, 1.50 at Ha = 0 and Re = 10 at (a) = (b) = + ,

   (c) = + 2, (d) = + 3.

   the cycle repeats. Shear thinning fluid responds quickly to the change in input velocity due to reduced viscosity but shear thickening fluid shows the least variation on account of increased viscosity. As the Reynolds number is increased to Re

   = 50, differences between velocity profiles become considerably smaller as illustrated in Fig. 13. The influence of fluid rheology weakens and profiles corresponding to different power-law indices begin to merge. This happens because the fluid does not get enough time to maintain its state from previous time instant, as the flow is reversed in next time instant. A further increase in Reynolds number to Re = 100 results in an almost complete overlap of velocity profiles for all three fluids as demonstrated in Fig. 14. This indicates that the flow is controlled by inertial forces and under these conditions, Newtonian and non-Newtonian fluids exhibit very similar responses.

   The influence of transverse magnetic field on oscillatory flow behaviour is investigated by setting Ha = 20. Lorentz force generated suppresses the oscillatory movement induced by sinusoidal inlet velocity. As a result, overall velocity magnitude decreases. Fig. 15, Fig. 16, and Fig. 17 compare the u-velocity profiles of shear thinning, Newtonian, and shear thickening fluids at different time instants for Ha = 20 and Re = 10, 50, and 100.

   (a)

   (b)

   (c)

   (d)

   Fig. 13. uvelocity profile comparison for oscillatory flow of power-law fluid with z = 1.00, 0.75, 1.50 at Ha = 0 and Re = 50 at (a) = (b) = + ,

   (c) = + 2, (d) = + 3.

   The presence of magnetic field introduces damping effect throughout the flow domain, which opposes oscillatory motion and weakens velocity variation. Consequently, velocity profiles become more confined and exhibit smaller amplitudes compared to non-magnetic case. Similar to the non-magnetic case, influence of inertia reduces the effect of apparent viscosity with an increase in Reynolds number to Re = 50 and 100 as shown in Fig. 16 and Fig. 17 respectively.

   (a)

   (b)

   (c)

   (d)

   Fig. 14. uvelocity profile comparison for oscillatory flow of power-law fluid with z = 1.00, 0.75, 1.50 at Ha = 0 and Re = 100 at (a) = (b) = +

   , (c) = + 2, (d) = + 3.

   (a)

   (b)

   (c)

   (d)

   Fig. 15. uvelocity profile comparison for oscillatory flow of power-law fluid with z = 1.00, 0.75, 1.50 at Ha = 20 and Re = 10 at (a) = (b) = +

   , (c) = + 2, (d) = + 3.

   (a)

   (b)

   (c)

   (d)

   Fig. 16. uvelocity profile comparison for oscillatory flow of power-law fluid with z = 1.00, 0.75, 1.50 at Ha = 20 and Re = 50 at (a) = (b) = +

   , (c) = + 2, (d) = + 3.

   (a)

   (b)

   (c)

   (d)

   Fig. 17. uvelocity profil comparison for oscillatory flow of power-law fluid with z = 1.00, 0.75, 1.50 at Ha = 20 and Re = 100 at (a) = (b) = +

   , (c) = + 2, (d) = + 3.

4. CONCLUSION

In this work, finite volume based numerical method is developed and validated for simulating pressure-driven, shear- driven, and oscillatory flow of power-law fluids in straight channel. The results show that Reynolds number and power-law index significantly influence velocity distribution and vorticity generation.

In pressure-driven flow, increasing Re reduces boundary layer thickness before getting fully developed and intensifies velocity gradients near the walls. Shear thinning fluids exhibit lower apparent viscosity in high shear regions, resulting in thinner boundary layer and more concentrated vorticity, whereas shear thickening fluids display thicker boundary layer and smoother vorticity distribution. The introduction of transverse magnetic field generates Lorentz forces opposing fluid motion, leading to suppression of velocity magnitude and vorticity. Velocity profiles become increasingly plug-like with the formation of Hartmann boundary layers near the walls.

In shear-driven flow without magnetic field, all fluids produce identical linear velocity profiles, indicating that fully developed velocity distribution is independent of Re and z. This happens due to constant shear stress across the channel width for all fluids. With magnetic field, shear thickening fluid shows highest velocity at a particular y-coordinate followed by Newtonian and shear thinning fluid due to inherent rheological behaviour.

In oscillatory flow without magnetic field, sinusoidal velocity input imposes fluctuations in the flow field at every time instance. At low Re, the velocity fluctuations are apparent, whereas at high Re the fluctuations diminish due to rapid inversion in velocity magnitude. Shear thinning fluids respond more rapidly to oscillatory forcing owing to their reduced effective viscosity, while shear thickening fluids exhibit greater resistance to flow reversal and velocity fluctuations. Application of magnetic field dampens oscillatory motion, reduces velocity amplitudes, and promotes flow stabilization throughout the oscillation cycle.

ACKNOWLEDGMENT

This research was supported by Science and Engineering Research Board, a statutory body of Department of Science and Technology (DST), Government of India through the funded project CRG/2022/004030.

REFERENCES

1. S. Chaudhuri, and S. Sahoo, Characterization of fully developed pressure-driven, shear-driven and combined pressure and shear driven flow of sisko fluids through rectangular channels, Arabian Journal for Science and Engineering, vol. 45(7), pp. 5925-5947, July 2020.

2. G. A. Garang, J. Okelo, and K. Giterere, Unsteady MHD flow of non- Newtonian fluid between two parallel plates with inclined variable magnetic field, Global Journal of Pure and Applied Mathematics, vol. 18(1), pp. 131-154, 2022.

3. R. Ershadnia et al., Non-Newtonian fluid flow dynamics in rotating annular media: physics-based and data-driven modeling, Journal of Petroleum Science and Engineering, vol. 185, 106641, Februrary 2020.

4. S. Chaudhuri et al., Thermal characteristics of forced convection in combined pressure and sheardriven flow of a nonNewtonian thirdgrade fluid through parallel plates, Heat Transfer, vol. 50(7), pp. 6737-6756, November 2021.

5. R. Liu, and Q. S. Liu, Non-modal instability in plane Couette flow of a power-law fluid, Journal of fluid mechanics, vol. 676, pp. 145-171, June 2011.

6. D. O. Makinde, and O. Franks, Thermal decomposition of unsteady non- Newtonian MHD Couette flow with variable properties, International Journal of Numerical Methods for Heat & Fluid Flow, vol. 25(2), pp. 252- 264, March 2015.

7. S. Rath, and B. Mahapatra, Low Reynolds number pulsatile flow of a viscoelastic fluid through a channel: effects of fluid rheology and pulsation parameters, Journal of Fluids Engineering, vol. 144(2), 021201, February 2022.

8. Z. Bukhari, A. Ali, Z. Abbas, and H. Farooq, The pulsatile flow of thermally developed non-Newtonian Casson fluid in a channel with constricted walls, AIP Advances, vol. 11(2), 025324, February 2021.

9. R. Maniyeri, Y. K. Suh, S. Kang, and M. J. Kim, Numerical study on the propulsion of a bacterial flagellum in a viscous fluid using an immersed boundary method, Computers & fluids, vol. 62, pp. 13-24, June 2012.

10. D. S. Kershaw, The incomplete Cholesky-conjugate gradient method for the iterative solution of systems of linear equations, Journal of computational physics, vol. 26(1), pp. 43-65, January 1978.

11. Z. Chen, and C. Shu, Simplified lattice-Boltzmann method for non Newtonian powerlaw fluid flows, International Journal for Numerical Methods in Fluids, vol. 92(1), pp. 38-54, January 2020.