Turbulence transition in pulsating pipe flow: shear-thinning effects

Hola Universitat Polit`ecnica de Catalunya Treball Final de Grau en Enginyeria F´ısica Turbulence transition in pulsatile pipe flow: shear-thinning effects ZARM - University of Bremen June 2021 Author: Laura Dot Casadevall Director: Prof. Dr. Marc ´ Avila Codirector: Dr. ´ Alvaro Meseguer Supervisor: Daniel Mor´on Acknowledgements I would like to express my gratefulness to Prof. Dr. Marc ´ Avila for allowing me to join the Center of Applied Space Technology and Microgravity (ZARM) at University of Bremen and encouraged me to keep researching on pipe flow and turbulence transition topic. Moreover, I would really like to thank Daniel Mor´on, who has always helped me with everything I needed, no matter what and when, and he has supervised my work better than ever. I would like to extend my gratefulness to the other research group members for sharing their points of view and different questions about this work topic, as this has uncovered other perspectives to address the problem. Finally, I also want to thank my co-director Dr. ´ Alvaro Meseguer from Universitat Polit`ecnica de Catalunya (UPC), given that his lessons on Numerical and Computational Methods, plus Physics of Fluids inspired me to choose this bachelor’s thesis topic and he helped me to contact Prof. Dr. Marc ´ Avila. 1 Notation dictionary ttime u=(u,v,w) velocity field ρdensity ppressure τshear stress tensor νkinematic viscosity µviscosity ˆµ∞infinite shear-rate viscosity ˆµ0zero shear-rate viscosity ˆ λfluid’s time constant ncshear-thinning index ashift term Re Reynolds number Wo Womersley number wpulsation angular frequency Apulsation amplitude Rpipe radius Dpipe diameter Nnumber of radial nodes αaxial wavenumber βazimuthal wavenumber Genery gain 2 Contents Abstract 6 1 Introduction 7 1.1 Motivation and state of the art . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Objectives.................................. 7 2 Description of the problem 8 2.1 Non-NewtonianFluids ........................... 8 2.1.1 Carreau-Yasuda model . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Governingequations ............................ 10 3 Base Flow 12 3.1 SteadyBaseFlow.............................. 12 3.1.1 Chebyshev discretization . . . . . . . . . . . . . . . . . . . . . . 13 3.1.2 Iterative Newton’s method . . . . . . . . . . . . . . . . . . . . . 14 3.2 UnsteadyBaseFlow ............................ 15 4 Transient Growth Analysis (TGA) 17 4.1 Disturbance and linearized Navier-Stokes equations . . . . . . . . . . . 17 4.2 Theoretical fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2.1 Linearstability ........................... 19 4.2.2 Transient Growth assessment . . . . . . . . . . . . . . . . . . . 19 4.3 Improved Projection Scheme (IPS) . . . . . . . . . . . . . . . . . . . . 20 4.3.1 IPS for Newtonian fluids . . . . . . . . . . . . . . . . . . . . . . 20 4.3.2 IPS for non-Newtonian fluids . . . . . . . . . . . . . . . . . . . 22 5 Convergence Analysis 23 6 Results and discussion 27 6.1 Baseflowanalysis.............................. 27 6.2 TGAresults................................. 32 7 Conclusions 37 References 38 3 List of Figures 1 Fluid behaviour according to (a) shear stress, (b) viscosity vs shear rate relation ................................... 8 2 Carreau viscosity profile . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 Steady base flow validation . . . . . . . . . . . . . . . . . . . . . . . . . 13 4 Newton’s interative method flowchart. . . . . . . . . . . . . . . . . . . 15 5 Unsteady base flow validation. . . . . . . . . . . . . . . . . . . . . . . . 16 6 Velocity (left) and viscosity (right) profiles for Newtonian fluid in two different instants (t=0 and t=0.5T) for Wm=15, Re=2000 and A=1. . 23 7 Velocity and its corresponding derivative profiles for t=0.5T Newtonian fluid with Wm=15, Re=2000 and A=1. . . . . . . . . . . . . . . . . . . 24 8 Convergence error for λ=0.5 and µ∞=0.0616, using equispaced nodes in finite elements with stencil 7. . . . . . . . . . . . . . . . . . . . . . . . 24 9 Minimum radial nodes (Nmin) required to fulfill viscosity convergence for different λand µ∞values. ....................... 25 10 Viscosity profiles for different λwith each respective Nmin required for a good shape definition. . . . . . . . . . . . . . . . . . . . . . . . . . . 26 11 Blood steady laminar velocity and viscosity profiles. . . . . . . . . . . . 27 12 Base flow velocity and viscosity profiles for fixed µ∞=0.2, nc=0.35, a=2, Re=2000 and changing λ.......................... 28 13 Base flow velocity and viscosity profiles for fixed λ=5, nc=0.35, a=2, Re=2000 and changing µ∞......................... 28 14 Velocity and viscosity profiles for λ=5, Re=2000, Wo=15 and changing µ∞for different times. Red dashed lines correspond to Newtonian case, blue and green to non-Newtonian case with µ∞=0.2 and µ∞=0.0616 repectively. ................................. 29 15 Velocity and viscosity profiles for µ∞=0.0616, Re=2000, Wo=15 and changing λfor different times. Red dashed lines correspond to Newtonian case; blue, green and pink to non-Newtonian case with λ=5, λ=20 and λ=200repectively............................ 30 16 Velocity and viscosity profiles for λ=5, Re=2000, µ∞=0.2 and changing Wo for different times. Red dashed lines correspond to Newtonian case; blue, green and pink to non-Newtonian case with Wo=15, Wo=12.5 and Wo=10repectively.............................. 31 17 m=0,1,2... Pulsatile base flow Descending Aorta velocity profiles for different instants. λ=242.2631, µ∞=0.0616 and Wo=15. Red dashed and blue lines correspond to Newtonian fluid and Desc. Aorta respectively. 31 18 m=0,1,2... Pulsatile base flow Left Carotid Artery velocity profiles for different instants. λ=722.4477, µ∞=0.0616 and Wo=5. Red dashed and pink lines correspond to Newtonian fluid and Left Carotid Artery respectively.................................. 32 19 m=0,1,2... Pulsatile base flow Cerebral Capillaries velocity profiles for different instants. λ=1308.6350, µ∞=0.0616 and Wo=2. Red dashed and green lines correspond to Newtonian fluid and Cerebral Capillaries respectively. ................................ 32 4 20 Optimal (a) Newtonian and (b) non-Newtonian perturbations at maximal growth (t=1T) on top of non-Newtonian base flow with λ=0.06, Wo=15, Re=2000, A=1, a=2, nc=0.35 and µ∞=0.2. . . . . . . . . . . 33 21 Energy plot for (a) Newtonian and (b) non-Newtonian perturbations on top of non-Newtonian base flow with λ=0.06, Wo=15, Re=2000, A=1, a=2, nc=0.35 and µ∞=0.2. ........................ 33 22 Upper right plot: Optimal perturbations at maximal energy growth for (a) λ=5, t=0.9T and (b) λ=20, t=0.8T. Upper left plot: contour of vorticity on a r-θcross-section. Bottom plot: time series of the energy growth of the optimal (a) helical (α,β)=(3,1) and (b) classic (α,β)=(0,1) disturbances at (Re,A,Wo)=(2000,1,15). . . . . . . . . . . . . . . . . . 34 23 Maximal energy growth with respect (a) Reynolds number for different λ, (b) fluid’s time constant for different Wo and (c) fluid’s time constant for different Wo and µ∞. ......................... 35 24 (a) Energy growth of the optimal stream-wise disturbance at (Re,A,Wo)=(2000,1,15) for the desc. aorta λ=242.2631 and µ∞=0.0616. (b) Time series of kinetic energy contribution of each velocity component for the optimal disturbance.................................. 36 25 Optimal perturbation at maximal energy growth (tf=1.7T) for the descending aorta λ=242.2631, µ∞=0.0616. Upper left plot: vorticity contour on a r-θcross-section. Bottom plot: time series of the energy growth fortwoperiods................................ 36 List of Tables 1 Carreau and Carreau-Yasuda blood parameters . . . . . . . . . . . . . 9 2 Steady base flow validation errors . . . . . . . . . . . . . . . . . . . . . 13 3 Unsteady base flow validation errors . . . . . . . . . . . . . . . . . . . . 16 4 Cardiovascular sections parameters . . . . . . . . . . . . . . . . . . . . 27 5 Newtonian and non-Newtonian perturbations TGA results. . . . . . . . 33 6 Desc. Aorta parameters and TGA results for the most dangerous perturbation to trigger turbulence transition. . . . . . . . . . . . . . . . . 35 5 Abstract In the last years the major cause of death worldwide has turned out to be because of cardiovascular diseases; for instance, owing to the presence of blood turbulence in our vessels. The problem has already been addressed from simplified points of view that ignore many blood characteristics; consequently, one of the main goals of this project is to take into account the blood shear-thinning effects in an infinite cylindrical pipe and study turbulence transition in such scenario. With this aim, the first thing to do is a non-Newtonian viscosity models research in order to find the best one to model blood behaviour. After this, an exhaustive base flow study is carried out, for steady and unsteady scenarios, so that finally the Transient Growth Analysis (TGA) provides us with the most dangerous perturbation that could trigger turbulence transition more easily in our veins. Besides, a parametric study has been performed to take into account a wide range of different possibilities and to be able to highlight more general conclusions. 6 1 Introduction 1.1 Motivation and state of the art In 1883 [16], Reynolds identified the basic problem in the context of flow in a pipe: when and how do high-speed flows undergo transition from laminar state to more complicated scenarios such as puffs, slugs and turbulence? To this day, the problem of the control of fluid turbulence has been under research [18]; often, but not always, to delay its occurrence or mitigate its effects. In addition, it has many practical applications ranging from aeronautics [17] to pipeline engineering [2]. Understanding the phenomenon of turbulence transition initiated by the growth of infinitesimal perturbations is a necessary prerequisite to subsequently find effective transition-delay strategies. This bachelor’s thesis is mainly focused on a more biological background. Cardiovascular diseases represent one, if not the major, cause of death worldwide [20]. With its ever increasing danger in mind, the scientific community is interested on determining its causes for later to be able to prevent them. One of these causes appears to be the presence of turbulence in our cardiovascular system [1]. The only solid understanding regarding transition to turbulence in our arteries or veins comes from the analogy between cardiovascular flow and the elementary Newtonian fluid flow in rigid circular cross-section pipes. It comes with no surprise that, with such a simplified model, many characteristics of blood flow in our flexible vessels are ignored. Most recently, the effects that a pulsating driving flow (instead of a steady driving) has on the transition scenario has been uncovered for a Newtonian fluid [6], [11]. As further steps, it would be interesting to check which are the effects that other previously ignored characteristics of cardiovascular flow have on turbulence transition, specifically shear-thinning effects. 1.2 Objectives As said before, the main intention for this project is to take a step forward in the blood transition to turbulence problem and consider this fluid as non-Newtonian with shear-thinning effects. That is why the main goals of this thesis are: - First of all, perform a literature research of existing shear thinning fluids models and turbulence transition in order to find the best model to assess this project’s problem. - Secondly, implement a Matlab code to obtain steady/unsteady laminar flows; as well as, implement a further script to perform Transient Growth Analysis. - Afterwards, an analysis of the flow behaviour and structures in the scope of the linear Navier-Stokes equations is expected to be done in order to find the most dangerous perturbation in each case. - Finally, the concluding goal would be to perform a parametric study for different shear-thinning conditions, not forgetting the corresponding discussion and conclusions. 7 2 Description of the problem In order to model pipe flow for non-Newtonian shear thinning fluids an infinite circular cross-section pipe of radius ˆ Ris considered, where an incompressible fluid flows, the blood (constant fluid density). This last assumption considerably simplifies the equations that must be solved and furthermore, we are dealing with a viscous fluid. Dimensional parameters are expressed with hat (ˆ), so the three-dimensional velocity field vector at a given time tis defined as ~ ˆu(t)=(ˆu,ˆv, ˆw)(t), being ˆu, ˆvand ˆwthe radial, azimuthal and axial velocity components respectively. Therefore, in this whole project cylindrical coordinates are used owing to problem’s geometry, with the corresponding derivatives expressed as subscripts. For instance: ˆuris the radial derivative of ˆucomponent and ˆvθis the azimuthal derivative of ˆv. 2.1 Non-Newtonian Fluids The well known scientist Isaac Newton owes his fame to various discoveries and developed theories in mathematics and physics. Among these, he perceived that fluids behave with constant viscosity that could only changes with temperature [19]. However, it was uncovered that stress applied to a fluid plays an important role. Non-Newtonian fluids change their viscosity or flow behaviour under stress, which means that their viscosity is dependent on shear rate; so they are able to get thicker and act more like a solid, or on the other hand, the fluid can get runnier than before. Hence, not all non-Newtonian fluids behave the same way when stress is applied. To begin with, Newtonian fluids follows a linear relation between shear stress (τ) and shear rate (˙γ): τ=µ·˙γ(2.1) so viscosity (µ) is a constant value, no matter how fast the fluid is forced to flow through a pipe or channel. In contrast, non-Newtonian fluids display a non-linear relation between shear stress and shear rate, as it can be seen in the following figure. (a) (b) Figure 1: Fluid behaviour according to (a) shear stress, (b) viscosity vs shear rate relation So a fluid is said to be shear thickening when its viscosity increases as the shear rate does and is said to be shear thinning if its viscosity decreases for higher ˙γ. However, this project is mainly focused on shear-thinning fluids, since we want to analyse blood 8 Figure 4: Newton’s interative method flowchart. 3.2 Unsteady Base Flow Unsteady base flow corresponds to the scenario when the laminar flow changes over time; specifically, in this project a pulsatile model has been chosen such that the imposed pressure gradient constraint meets the following flow rate relation, Q=um-Acos(wt), being ’um’ the mean velocity; ’A’ and ’w’ the amplitude and angular frequency of the pulsation respectively. Transition to turbulence in this case also depends on a new parameter called Womersley number Wo, which represents the relation between pulsation frequency (f) and viscous effects: Wo =D 2r2πf ν(3.9) So now, since the problem is time dependent, an additional term corresponding to the laminar velocity time derivative must be added to the equation 3.1. Re∂wb ∂t =ReGp+Lwb+CNL ·NL (3.10) Time discretization has been chosen to be carried out by a second order backward differentiation formula (BDF2), which is an implicit method for numerical integration. For a given function and time, this multistep method approximates the respective derivative using previously computed values. Hence, the corresponding BDF2 formula for the problem y’=f(t,y) would be: yn+2 −4 3yn+1 +1 3yn=2 3hf(tn+2, yn+2) (3.11) where h denotes the integration step size. Bearing this in mind, as well as recalling that ’Gp’ and ’wb’ are the unknown variables, the unsteady flow Jacobian reads as: J=0 2 ·r·wei −2 3dtRe Re −2 3dt[L+CNL ·(µp∂ ∂r +∂wb ∂r ∂µp ∂wb)](3.12) 15 Finally, the pulsatile Matlab code validation has also been successfully done for the Newtonian case, such that the computed profile has been validated with Sexl- Womersley exact analytical solution. For instance, the following plots correspond to a Newtonian fluid with parameters: A=1, Wo=15 and Re=2000 in two different times. (a) t=mT, m=0,1,2... (b) t=(2m+1)T/2, m=0,1,2,... Figure 5: Unsteady base flow validation. Once again, the mean square error corresponding to both different times has been computed; so taking into account that m=0,1,2,... the respective mean square errors for N=80 are: t=mT t=(2m+1)T/2 MSE (m/s)21.46 ·10−69.95 ·10−7 Table 3: Unsteady base flow validation errors 16 4 Transient Growth Analysis (TGA) In this section, the transient evolution of perturbations on top of the previous base flows is considered. So first of all, the corresponding linearization of the Navier-Stokes equations will be seen, followed by a brief description of the theoretical fundamentals and as a last section, also a short and concise explanation of the method used to implement TGA. 4.1 Disturbance and linearized Navier-Stokes equations ”Hydrodynamic stability theory is concerned with the response of a laminar flow to a disturbance of small or moderate amplitude” [10]. Generally, disturbances are considered infinitesimal so that further simplifications could be justified; for instance, the governing linear equation of the disturbance evolution is yearned. However, nonlinear effects turn out to be significant once the disturbance velocities grow above the base flow; at this point, the linear equations are not completely precise to predict the disturbance evolution. Despite the fact that the linear equations have a limited region of validity they are of great importance for physical growth mechanisms detection, as much as dominant disturbance types identification. Thus, the base flow is initially perturbed by an infinitesimal disturbance (u,p) such that: U=~ub+uand P=pb+p, where variables ’~ub’ and ’pb’ correspond to the base flow. The governing time evolution equation of the disturbance is obtained by introducing these perturbations to the momentum equation, keeping the lowest order (linear) terms in and setting quadratic terms to zero: ∂u ∂t + (~ub∇)u+ (u∇)~ub=−∇p+1 Re∇(~τub+u−~τub) (4.1a) ∇u= 0 (4.1b) Regarding the solution of 4.1, two main questions have to be considered. The first one states: is there a solution that grows boundless? If so, the base flow ~ub is linearly unstable. The second question stands for: given a linearly stable base flow, is there any bounded solution that manifest large transient growth before the inevitable decay? In these cases, it can be considered that the base flow is linearly stable; however, susceptibility to nonlinear instability is present, as well as some local regions of convective instability. Hence, when answering both questions, it is wished to find the perturbation that produces optimal growth. As noted before, an infinitesimal disturbance is applied; therefore, the term standing for the shear stress tensor ~ τ0=~τub+u−~τubhas to be linearized around the base flow (~ub,pb). To this purpose, detailed analysis of each term is performed, starting with viscosity term µ(~ub+u), which is linearized by a Taylor expansion. It should be called up that ~τ(~u)=µ(~u) ˙γ(~u), besides taking into account that ˙γis a linear term so: ˙γ(~ub+u) = ˙γ(~ub) + ˙γ(u). µ(~ub+u) = µ( ˙γb+ ˙γ(u)) = µb+1 2 ∂µ ∂˙γij b : ˙γij(u) = µb+∂µ ∂˙γrz b∂w ∂r +∂u ∂z (4.2) 17 Introducing 4.2 in ~ τ0, as well as neglecting infinitesimal terms: ~ τ0=~τub+u−~τub=µ(~ub+u) [ ˙γ(~ub) + ˙γ(u)] −µb˙γb=µb˙γ(u) + ∂µ ∂˙γrz b∂w ∂r +∂u ∂z ˙γb (4.3) ∂µ ∂˙γrz b = (1 −µ∞)(nc−1)λa1 + λa ∂wb ∂r  anc−1 a−1 ∂wb ∂r  a−1∂wb ∂r |∂wb ∂r |(4.4) Once the disturbance shear stress tensor has been worked out, next step would be applying the divergence at each term: ∇(µb˙γ(u)) = µb4u+∂µb ∂r ,0,0˙γ(u) = µb4u+  2ur∂µb ∂r [vr+1 r(uθ−v)]∂µb ∂r [uz+wr]∂µb ∂r  (4.5) ∇∂µ ∂˙γrz b (wr+uz) ˙γb=∇       0 0 ∂µ ∂˙γrz b (wr+uz)∂wb ∂r 000 ∂µ ∂˙γrz b (wr+uz)∂wb ∂r 0 0       = =       ∂ ∂z ∂µ ∂˙γrz b (wr+uz)∂wb ∂r 0 ∂ ∂r ∂µ ∂˙γrz b (wr+uz)∂wb ∂r      (4.6) In accordance, shear stress tensor divergence would write as: ∇~ τ0=µb4u+   2ur∂µb ∂r +∂µ ∂˙γrz b(wrz +uzz)∂wb ∂r [vr+1 r(uθ−v)]∂µb ∂r [uz+wr]∂µb ∂r + [A1] ∂wb ∂r +∂µ ∂˙γrz b(wr+uz)∂2wb ∂r2   (4.7) A1 = ∂ ∂r ∂µ ∂˙γrz b(wr+uz) + ∂µ ∂˙γrz b (wrr +uzr) (4.8) Finally, after having developed the other terms from 4.1a, plus splitting in velocity components, the linearized Navier-Stokes equations (LNSE) for an infinitesimal disturbance are obtained: ut+wbuz=−pr+1 Re µb4u−u r2−2 r2vθ+ 2ur ∂µb ∂r +∂µ ∂˙γrz b(wrz +uzz)∂wb ∂r  (4.9) vt+wbvz=−1 rpθ+1 Re µb4v−v r2+2 r2uθ+vr+1 r[uθ−v]∂µb ∂r (4.10) 18 wt+u∂wb ∂r +wbwz=−pz+1 Re µb4w+ (wr+uz)∂µb ∂r +A1∂wb ∂r +∂µ ∂˙γrz b(wr+uz)∂2wb ∂r2 (4.11) Equations 4.1b, 4.9, 4.10 and 4.11 are linear, periodic in θand unbounded in z (due to the cylindrical domain); so moreover, as introduced in Meseguer et al. paper [4], all solutions to these equations can be expressed as superposition of complex Fourier modes of the form u(r, θ, z, t) = ei(βθ+αz)u(r, t); p(r, θ, z, t) = ei(βθ+αz)p(r, t) (4.12) where β∈Zand α∈Rare the azimuthal and axial wavenumbers respectively. Introducing 4.12 into 4.9, 4.10 and 4.11 gives simplified equations that are easier to implement in Matlab. 4.2 Theoretical fundamentals 4.2.1 Linear stability As distinctively explained in [5], a linear evolution operator Λ(t) could be defined to represent the disturbance evolution forward in time: u(t)=Λu(0) (4.13) Through this forward evolution operator it is easy to access the linear stability of the base flow. Two different scenarios could be distinguished: for a steady base flow, LNSE are autonomous with eigenmode solutions of the form u(x,t)=exp(λjt)~uj(x)+C; for an unsteady periodic base flow (with period T), LNSE becomes a Floquet problem with solutions of the form u(x,t)=exp(λjt)~uj(x,t)+C, where ~uj(x,t+T)=~uj(x,t). In either case, an eigenvalue problem in terms of the operator Λ(T) is obtained once time t is settled to a value tf. Λ(tf)~uj=ηj~uj, ηj=exp(λjtf) (4.14) The dominant eigenvalues of Λ(tf), ηjof the largest modulus, are used to determine the classical linear stability of the base flow. – If any ηjcomply with |ηj|>1, there exist exponentially growing solutions and the base flow is linearly unstable. – If every ηjcomply with |ηj|<1, all solutions decay to zero and the base flow is linearly stable. –|ηj|=1 is a sign of bifurcation point 4.2.2 Transient Growth assessment In particular, what is about to be studied is the second main question introduced in the linearized equations section (4.1): given a linearly stable base flow, is there any bounded solution that manifest large transient growth before the inevitable decay? 19 Perturbations will exhibit substantial transient response due to regions of localized convective instability, as [5] corroborate. The asymmetry of the convective terms in the Navier-Stokes equations leads to non-orthogonal eigenmodes; the aforementioned affects the form of the dynamics of interest, which in fact won’t be an exponential function of time multiplying a fixed modal shape. Hereafter, transient growth analysis is going to be the way to quantify such dynamics. Transient growth is assessed as the normalized perturbation energy at time τto its initial energy. Considering normalized initial perturbations, a new parameter known as energy growth (G), the one that is optimized in this method, is defined as the ratio between the total perturbation kinetic energy at time τand t0: G(τ) = E(τ) E(0) =||u(τ)||2= (u(τ),u(τ)) (4.15) Where the standard L2inner product is present: (u,v) = ZΩ uvdυ (4.16) Introducing 4.13 into 4.15 leads to: E(τ) E(0) = (Λ(τ)u(0),Λ(τ)u(0)) = (u(0),Λ∗(τ)Λ(τ)u(0)) (4.17) Being Λ∗(τ) the adjoint evolution operator of Λ(τ). Henceforth, the main question should be remodeled: which is the disturbance that leads to maximal growth? For the stability problem, the greatest importance falls on the dominant eigenvalues, since they stand for the ”worst case scenario”, so the eigenfunctions corresponding to the dominant eigenvalues (λj) of Λ∗(τ)Λ(τ) will describe the largest possible growth. Gmax(τ) = maxE(τ) E(0) =maxλj(4.18) 4.3 Improved Projection Scheme (IPS) In this subsection, the numerical method implemented to compute transient growth is explicitly detailed. This specific scheme has been chosen since its efficiency and outstanding temporal behaviour for sufficiently large time integration are a key feature to capture three-dimensional simulations of unsteady flows. One of the major problems when solving the incompressible Navier-Stokes equations concerns the coupling of pressure and velocity to satisfy the incompressibility constraint. Indeed, the pressure gradient varies for time-dependent flows while preserving its boundary value fixed, held at its initial value, during time integration. To assess this, a pressure predictor is introduced to provide a suitable pressure field with a divergence free velocity. As it is also highlighted in [12], the temporal scheme’s accuracy order used is kept for the pressure and velocity when using this algorithm. 4.3.1 IPS for Newtonian fluids Considering an incompressible Newtonian fluid, the respective governing Navier-Stokes equations: ∂u ∂t + (u∇)u=−∇p+1 Re4u(4.19) 20 Distinguishing the linear and non-linear terms, the former equation could be expressed as: ∂u ∂t +NL(u) = −∇p+1 ReL(u) (4.20) ∇u= 0 (4.21) Applying divergence on the Navier-Stokes, plus taking into account that velocity is divergence free (∇u=0), the Poisson equation is derived: ∇NSE → ∇(∂u ∂t ) + ∇NL(u) = −4p+1 Re∇L(u)→ ∇2p=−∇NL(u) (4.22) However, a consistent boundary condition is required to solve the Poisson equation, which can be derived from 4.20 and 4.21: ∂p ∂ˆn= ˆn−∂u ∂t −NL(u) + 1 ReL(u)(4.23) As remarkably mentioned in [12], the time accuracy of the global solution is directly dependent on the treatment of these pressure boundary conditions. So now, it should be remembered that spatial approximation is accomplished by Chebyshev discretization and the temporal scheme obey an Adams-Bashforth and BDF combination of secondorder accuracy, such that equation 4.20 reads as: 3un+1 −4un+un−1 2∆t+ 2NL(un)−NL(un−1) = −∇pn+1 +1 Re4(un+1) (4.24) By such means, the main steps to follow in order to get the perturbation components proceed as: (i) Pressure pn+1 prediction: 4pn+1 =−∇2NL(un)−NL(un−1)(4.25) ∂pn+1 ∂ˆnwall = ˆn−3un+1 bc + 4un−un−1 2∆t−2NL(un) + NL(un−1) + 1 Re(2L(un)−L(un−1)) (4.26) Where ubc stands for the velocity boundary conditions. (ii) Velocity u∗prediction from the momentum equation, including the new pressure field and with boundary condition u∗=un+1 bc . 3u∗−4un+un−1 2∆t+ 2NL(un)−NL(un−1) = −∇pn+1 +1 Re4(u∗) (4.27) (iii) Intermediate step. An intermediate variable φ, defined as φ=2∆t 3(p∗(n+1) −pn+1), is introduced when the final divergence-free velocity field is explicitly evaluated, such that: 3un+1 −3u∗ 2∆t=−∇(p∗(n+1) −pn+1) (4.28) So taking the divergence of the former equation, a new Poisson equation for φ arise 4φ=∇u∗(4.29) 21 which can be solved together with the consistent Neumann boundary condition, as a means to get φ: ∂φ ∂ˆn= 0 (4.30) (iv) As a last step, update the corrected velocity field: un+1 =u∗−∇φ(4.31) Finally, this scheme is applied to the linear Navier Stokes equations for a Newtonian fluid and their corresponding adjoint to integrate forward and backward in time iteratively, until a Krylov space is built and from which is possible to obtain the optimal perturbation for each specific Wm, Re, λand A scenario. It has to be noted that TGA code was already implemented in [6] study; new contribution appears in the next section where the non-Newtonian terms will be taken into account to compute Transient Growth Analysis. 4.3.2 IPS for non-Newtonian fluids In this section, IPS for an incompressible non-Newtonian fluid is introduced. For this purpose, few line code modifications have been done in order to insert the non- Newtonian contribution to the model. Therefore, the equations leading this scenario are nothing more than the linearized Navier Stokes equations for a non-Newtonian fluid presented in section 4.1, corresponding to equations 4.9, 4.10 and 4.11, which could be expressed in the compact form (similar to eq. 4.20): ∂u ∂t +NL(u) = −∇p+1 Re [µbL(u) + NN(u)] (4.32) What is more, NN term can be treated just like NL. Hence, NL’(u)=NL(u)- 1 Re NN(u)and: ∂u ∂t +NL0(u) = −∇p+1 ReµbL(u) (4.33) Applying again divergence on the former equation, a different expression than the Poisson equation is obtained; since as quickly seen, the appearance of the viscosity (µb) prevents the linear term from vanish. ∇(∂u ∂t )+∇NL0(u) = −4p+1 Re∇(µbL(u)) → ∇2p=−∇NL(u)+ 1 ReL(u)∇µb(4.34) The consistent boundary condition required to solve the previous equations reads as: ∂p ∂ˆn= ˆn−∂u ∂t −NL0(u) + 1 ReµbL(u)(4.35) Hereafter, the following discretization treatment keeps the same as the IPS Newtonian case, as well as the four main steps. 22 5 Convergence Analysis As well argued in section 3.1.1, the discrete rendering of the problem’s geometry is a crucial key when solving CFD tasks in view of the fact that convergence and solution accuracy depend on it. During the base flow results analysis, it has been surprising to discover that the viscosity profile is not at all well defined, wiggles appear throughout the whole domain. Hence, the main goal of the following analysis is to determine the minimum number of radial nodes (Nmin) needed to satisfactory capture the viscosity behaviour for different increasing λvalues, which is the dimensionless fluid’s time constant. Besides, in this convergence analysis Chebyshev nodes are left out for finite differentiation using an equispaced grid of stencil 7. With the aim of finding Nmin, a fixed velocity profile, specifically corresponding to t=0.5T Sexl-Womersley base flow for Wo=15, Re=2000 and A=1, has been chosen to be the one for viscosity computation. It should be noted that any other velocity profile could be valid to carry out this convergence analysis; however, the interest falls on the profiles with more complex shapes. In other words, the interest falls on the ones with maximums/minimums near the pipe wall that would trigger the viscosity rapidly to 1 (because the velocity derivative would be zero), giving way to the appearance of sudden peaks in the viscosity profile. dWb dr = 0 →µ=µ∞+ (1 −µ∞)1 + λa·dWb dr a(nc−1)/a →µ= 1 (5.1) Hence, this is the reason why t=0.5T Sexl-Womersley profile is chosen and not t=0 for example. Figure 6 shows the appearance of the mentioned peak due to the minimum in the velocity profile between 0.8 <r<1. Moreover, the peak is not well defined and the purpose of this study falls on finding the minimum radial nodes Nmin such that this crest is well delineated. Figure 7 displays the chosen velocity profile with its corresponding radial derivative, which remains the same throughout the analysis. Figure 6: Velocity (left) and viscosity (right) profiles for Newtonian fluid in two different instants (t=0 and t=0.5T) for Wm=15, Re=2000 and A=1. 23 Figure 7: Velocity and its corresponding derivative profiles for t=0.5T Newtonian fluid with Wm=15, Re=2000 and A=1. So, starting from a fixed low λvalue and N=80, the viscosity profile (µb) is computed, from which the term C (used to assess convergence) can be determined: C(N) = Zr µ2 b(N)dr (5.2) Subsequently, while increasing the radial nodes it is possible to compare the new C value with the previous one, such that this Nrising procedure is continued until a constraint is fulfilled. In this specific case, the required restriction was for the last four C values to have an error between them smaller than 10−8. Once this is reached, it could be said that the Nmin value is found; for instance, a plot of the error tendency can be seen below, so that error decreases as increasing N. Figure 8: Convergence error for λ=0.5 and µ∞=0.0616, using equispaced nodes in finite elements with stencil 7. Afterwards, the same methodology is carried out but for a slightly higher λvalue; such that finally, convergence analysis has been performed for λ∈[0.5,25] with small dλ=0.5 increments. Moreover, the current study has also been completed for two different µ∞values; µ∞=0.0616 corresponding to the blood shear viscosity at infinite shear rate and µ∞=0.2, three times (approximately) the former one. 24 (a) t=mT, m=0,1,2... (b) t=(4m+1)T/4, m=0,1,2... (c) t=(2m+1)T/2, m=0,1,2... (d) t=(4m+3)T/4, m=0,1,2... Figure 16: Velocity and viscosity profiles for λ=5, Re=2000, µ∞=0.2 and changing Wo for different times. Red dashed lines correspond to Newtonian case; blue, green and pink to non-Newtonian case with Wo=15, Wo=12.5 and Wo=10 repectively. Finally, unsteady base flow profiles for the same cardiovascular sections presented at the beginning of this section (descending aorta, left carotid and cerebral capillaries) have also been computed. Different fluid time constants and Womersley numbers describe the blood problem, the only fixed parameter in this case is the infinite shear-rate viscosity. Even if we could expect a softer profile for the descending aorta, because it has the lowest time constant, what is actually seen is the narrowest behaviour near the wall. This is due to Womersley number, which is higher, and its influence on the profile shape is stronger than λ’s one. On the other side, cerebral capillaries are the ones with highest λbut lowest Wo, so the profile is the softest among all. (a) t=mT (b) t=(4m+1)T/4 (c) t=(2m+1)T/2 (d) t=(4m+3)T/4 Figure 17: m=0,1,2... Pulsatile base flow Descending Aorta velocity profiles for different instants. λ=242.2631, µ∞=0.0616 and Wo=15. Red dashed and blue lines correspond to Newtonian fluid and Desc. Aorta respectively. 31 (a) t=mT (b) t=(4m+1)T/4 (c) t=(2m+1)T/2 (d) t=(4m+3)T/4 Figure 18: m=0,1,2... Pulsatile base flow Left Carotid Artery velocity profiles for different instants. λ=722.4477, µ∞=0.0616 and Wo=5. Red dashed and pink lines correspond to Newtonian fluid and Left Carotid Artery respectively. (a) t=mT (b) t=(4m+1)T/4 (c) t=(2m+1)T/2 (d) t=(4m+3)T/4 Figure 19: m=0,1,2... Pulsatile base flow Cerebral Capillaries velocity profiles for different instants. λ=1308.6350, µ∞=0.0616 and Wo=2. Red dashed and green lines correspond to Newtonian fluid and Cerebral Capillaries respectively. 6.2 TGA results So far, how to approach the problem and a deep review of what is happening with pipe base flow for a non-Newtonian fluid, as the blood, have been seen. Following, in this section, the TGA results are analysed and discussed. However, before directly addressing the results, a brief refresh could be useful. The Transient Growth Analysis method performs a non-modal stability analysis making use of the linearized Navier Stokes equations and assuming small enough flow perturbations; linearization has been performed for Newtonian and non-Newtonian disturbances (as seen in section 4.3). Thus, in this framework the perturbations grow through linear mechanisms on top of the base flow profile. So, the most dangerous initial velocity perturbation is computed using TGA in terms of energy growth, out of a wide range of azimuthal and axial wavenumbers, is to say, out of all initial perturbation shapes, initial ’t0’ and final ’tf’ pulsation phases; such that the larger this energy gain is, the easier the transition to turbulence will be. Newtonian and non-Newtonian perturbations on top of the same base flow, for small fluid’s time constants, have been analyzed; what is captured is basically the same optimal disturbance. The helical shape is preserved for both scenarios, as well as the initial and final phase to obtain maximal growth. The only term that differs concerns the energy growth, such that non-Newtonian perturbation is slightly more amplified than the Newtonian one (see table 5). However, as the results are quite the same, proved with the following figures (20 and 21), and it is known that for higher 32 λthe viscosity problem is not completely well defined and it would require too many discretization nodes, hereafter only Newtonian perturbations on top of non-Newtonian base flows are considered. For the following table concerning the Newtonian and non-Newtonian perturbations TGA results, remember that αand βstand for axial and azimuthal wavenumbers respectively. Gmax α β t0/T tf/T Newtonian 1064.6160 3 1 0.3 1 non- Newtonian 1104.7449 3 1 0.3 1 Table 5: Newtonian and non-Newtonian perturbations TGA results. (a) (b) Figure 20: Optimal (a) Newtonian and (b) non-Newtonian perturbations at maximal growth (t=1T) on top of non-Newtonian base flow with λ=0.06, Wo=15, Re=2000, A=1, a=2, nc=0.35 and µ∞=0.2. (a) (b) Figure 21: Energy plot for (a) Newtonian and (b) non-Newtonian perturbations on top of non-Newtonian base flow with λ=0.06, Wo=15, Re=2000, A=1, a=2, nc=0.35 and µ∞=0.2. Hence, the dependence of the Newtonian perturbation energy gain with respect of the different parameters that define the flow could be analyzed to highlight the most important features. Starting with Reynolds number, one could verify that for higher Re the maximum gain Gmax is also remarkably higher, fact that has already been proven at previous works as Duo Xu et al. [6] for Newtonian base flows. Thinking about this result plus remembering that Reynolds number is defined by Re=uD/ ν, a logical explanation can be found: higher Re implies lower ν, is to say, less viscous effects, such 33 that the perturbation has more freedom and less resistance to develop and amplify itself (see figure 23 (a)). Concerning the Womersley number, the behaviour is now the opposite: for higher Wo (Wo>10), less amplification. This can be related with the optimal perturbation shape; it has been proven that lower Wo corresponds to helical perturbation shapes (for all λvalues), while higher Wo lead to stream-wise geometries, the classical one, also for all λvalues. Therefore, it could be deduced that in a certain point, as Womersley number and λincrease, a change in the perturbation shape (from helical to streamwise) will be perceived. After a proper parametric analysis it has been found when this shape alteration occurs, for two different scenarios: for µ∞=0.0616 the perturbation reshaping takes place at Wo=10 and λ∈[5,20]; and for µ∞=0.2 the perturbation modulation takes place at Wo=12.5, also in the fluid’s time constant range λ∈[5,20]. So, as well as in [6], for sufficiently large Wo the classic disturbance dominates and the optimal gain Gof steady pipe flow is recovered. A graphic representation is attached in the following figures, for µ∞=0.0616, Wo=10 and Re=2000 scenario. (a) λ=5 (b) λ=20 Figure 22: Upper right plot: Optimal perturbations at maximal energy growth for (a) λ=5, t=0.9T and (b) λ=20, t=0.8T. Upper left plot: contour of vorticity on a r-θ cross-section. Bottom plot: time series of the energy growth of the optimal (a) helical (α,β)=(3,1) and (b) classic (α,β)=(0,1) disturbances at (Re,A,Wo)=(2000,1,15). For the parameter values of figure 22 (a), the optimal disturbance has a helical structure with (α,β)=(3,1) and is localized at the outer half of the pipe. The optimal point to disturb is during the deceleration phase at t0=0.5T, as well as the maximum amplification is reached also during the deceleration phase at t0=0.9T. In this case, the kinetic energy of the optimal helical perturbation is mostly distributed in the stream-wise and azimuthal components (which indicates a strong three-dimensional effect), which self-amplify rapidly during the deceleration phase and start decaying from t=0.9T on. Initially the disturbance spirals clock-wisely towards the pipe center; as the energy grows, the perturbation switches the spiraling direction. The classic α=0 , β=1 optimal disturbance observed in figure 22 (b) consists of stream-wise vortices, which means the energy is subsequently transferred to the streamwise velocity components, while the cross-stream ones decay monotonically. Overall the classic perturbation’s behavior appears to be rather insensitive to the change in flow profile throughout the cycle and the decay is very slow. In fact, this tendency was already noticed by Duo Xu et al. [6]. 34 The aforementioned result regarding the change of geometry for different µ∞may lead to deduce that the energy gain depends on the infinite shear rate viscosity too. Figure 23 (c) clearly shows the corresponding relationship between Gand µ∞; for lower viscosity values, less energy gain, so less danger to trigger transition to turbulence. Moreover, this result agrees with the one above, such that lower viscosity values correspond to lower energy amplifications; this could entail a sooner perturbation geometry change. (a) Wo=15, µ∞=0.2 (b) Re=2000, µ∞=0.2 (c) Re=2000 Figure 23: Maximal energy growth with respect (a) Reynolds number for different λ, (b) fluid’s time constant for different Wo and (c) fluid’s time constant for different Wo and µ∞. Finally, moving to the blood scenario one should pay more attention on the Descending Aorta cardiovascular section, since it is the one where turbulence transition is most critical regarding deaths. It has been found that the most dangerous perturbation, which can trigger easily transition to turbulence, consists of stream-wise vortices. This outcome completely agrees with the result previously found; perturbation geometry change from helical to classic disturbation for µ∞=0.0616 (blood scenario) occurs at Wo=10, and actually the descending aorta works for higher Womersley, as seen in the following table. λ µ∞Wo Gmax α β t0/T tf/T Desc. Aorta 242.2631 0.0616 15 201.2604 0 1 0 1.7 Table 6: Desc. Aorta parameters and TGA results for the most dangerous perturbation to trigger turbulence transition. Figure 24 (b) shows that the energy is transferred to the stream-wise velocity component, while the cross-stream components decay quickly and monotonically. The optimal point to disturb is during the acceleration phase, specifically just at the beginning of the period, and the maximum energy amplification is reached at the deceleration phase, after 1.7 periods from the initial phase: tf=1.7T. Finally, figure 25 displays the 3D perturbation shape at its maximum energy amplification point and contour vorticity on a r-θcross-section. 35 (a) (b) Figure 24: (a) Energy growth of the optimal stream-wise disturbance at (Re,A,Wo)=(2000,1,15) for the desc. aorta λ=242.2631 and µ∞=0.0616. (b) Time series of kinetic energy contribution of each velocity component for the optimal disturbance. Figure 25: Optimal perturbation at maximal energy growth (tf=1.7T) for the descending aorta λ=242.2631, µ∞=0.0616. Upper left plot: vorticity contour on a r-θ cross-section. Bottom plot: time series of the energy growth for two periods. 36 7 Conclusions To put the whole project in a nutshell, I could gladly state that the main goals have been successfully fulfilled. For instance, the state of the art study led us to use the Carreau-Yasuda model to approximate the viscosity behaviour; this model has been proven to be exceptionally outstanding for the steady driven case. Surprisingly, after the pulsatile flow study, Carreau-Yasuda law has turned out to not really be a good model to approximate the blood flow problem. This is because the corresponding shear-thinning blood parameters (ˆµ0, ˆµ∞,nc,a,ˆ λ) are determined from steady driven experiments (not for time dependent flows) where shear-thinning is applied until the flow stabilizes; after that, the results are fitted in a curve which meets the model. So, to get a better approximation for the blood flow problem, one should take into account new parameter values which fit the pulsatile scenario and other blood properties, as it could be the viscoelasticity. Refering to the practical part of the project, base flow Matlab codes have been implemented with flying colours, validating them with previous works as the one from Xu et al. [6] for Newtonian fluids and with Lopez-Carranza paper [15] for steady non-Newtonian case. What is more, the unsteady non-Newtonian case has also been covered and it has been discovered that the dimensionless fluid’s time constant λparameter has a huge impact on the laminar flow profile, which later will also be crucial for Transient Growth Analysis. Concerning the unsteady pulsatile case, the new parameter, Womersley number, takes possession and it is basically the one which rules the velocity profile and consequently, the viscosity profile, which has let us notice that Chebyshev discretization grid is not the best one to treat the problem, but using finite elements. Finally, Transient Growth study has uncovered that Newtonian and non-Newtonian perturbations introduce quite the same danger concerning turbulence transition (for small lambdas) in non-Newtonian base flows, such that the energy gain behaviour with respect Reynolds and Womersley numbers is the same as found by Duo Xu et al. [6]. For higher Reynolds, higher energy amplification and for higher Womersley, lower gain. Last but no least, despite the model’s limitations, the most important goal of this project was to find the most dangerous perturbation for turbulence transition in the descending aorta, in view of the fact that it is the cardiovascular section where turbulence is most critical regarding deaths. The optimal perturbation has turned out to be a stream-wise disturbance with maximum amplification after 1.7 periods since its introduction in the base flow. I truly think this research line is of deep importance in order to better understand what is going on in our vessels and also to be aware and try to prevent misfortunes. However further work and model improvement are still left to be done as a means to properly describe the problem and to be able to suitably modelize real blood behaviour; starting from choosing an upgraded non-Newtonian viscosity model, until an appropriate discretization scheme that can actually capture the precise blood behaviour. 37 References [1] Gotlieb AI. Cunningham KS. The role of shear stress in the pathogenesis of atherosclerosis. Lab Invest., 2005. [2] Michael Eckert. Pipe flow: a gateway to turbulence. Springer, 2020. [3] A. M. Robertson et al. Rheological models for blood. Springer, 2009. [4] A. Meseguer et al. Linearized pipe flow to Reynolds number 107. Journal of Computational Physics, 2003. [5] D. Barkley et al. Direct optimal growth analysis for timesteppers. International Journal for numerical methods in fluids, 2008. [6] Duo Xu et al. Non-modal transient growth of disturbances in pulsatile and oscillatory pipe flow. Journal of Fluid Mechanics, 2020. [7] J. S. Horner et al. Investigation of blood rheology under steady and unidirectional large amplitude oscillatory shear. The Society of Rheology, 2018. [8] M. Armstrong et al. Evaluating rheological models for human blood using steady state, transient, and oscillatory shear predictions. Springer, 2018. [9] N. Bessonov et al. Methods of Blood Flow Modelling. EDP Sciences, 2016. [10] P.J. Schmid et al. Stability and Transition in Shear Flows. Springer Science and Business Media New York, 2001. [11] R. Trip et al. An experimental study of transitional pulsatile pipe flow. American Institute of Physics, 2012. [12] S. Hugues et al. An improved projection scheme applied to pseudospectral methods for the incompressible Navier-Stokes equations. Int. J. Numer. Meth. Fluids, 1998. [13] S. Tabakova et al. Application of the Carreau Viscosity Model to the Oscillatory Flow in Blood Vessels. AIP Conf. Proc., 2015. [14] S. Tabakova et al. Oscillatory Carreau flows in straight channels. Royal Society Open Science, 2020. [15] S.N. L´opez-Carranza et al. Pipe flow of shear-thinning fluids. C.R. Mecanique, 2012. [16] Derek Jackson and Brian Launder. Osborne Reynolds and the Publication of His Papers on Turbulent Flow. Annu. Rev. Fluid Mech., 2007. [17] M. A. Leschziner and D. Drikakis. Turbulence modelling and turbulent-flow computation in aeronautics. The Aeronautical Journal, 2002. [18] John Lumley and Peter Blossey. Control of Turbulence. Annu. Rev. Fluid Mech., 1998. [19] I. Newton. Newton’s Principia : the Mathematical Principles of Natural Philosophy. New-York :Daniel Adee, 1642-1727. 38 [20] World Health Organization. The top 10 causes of death. [21] A. Sequeira and J. Janela. An overview of some mathematical models of blood rheology. Springer, 2007. [22] Lloyd N. Trefethen. Spectral methods in Matlab. SIAM, 2000. 39