TThis draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Active Flows on Curved Surfaces
M. Rank and A. Voigt † Institut f¨ur Wissenschaftliches Rechnen, TU Dresden, 01062 Dresden, Germany Center for Systems Biology Dresden (CSBD), Pfotenhauerstr. 108, 01307 Dresden, Germany Cluster of Excellence - Physics of Life, TU Dresden, 01062 Dresden, Germany(Received xx; revised xx; accepted xx)
We consider a numerical approach for a covariant generalised Navier-Stokes equation ongeneral surfaces and study the influence of varying Gaussian curvature on anomalousvortex-network active turbulence. This regime is characterised by self-assembly of finite-size vortices into linked chains of anti-ferromagnet order, which percolate through theentire surface. The simulation results reveal an alignment of these chains with minimalcurvature lines of the surface and indicate a dependency of this turbulence regime on thesign and the gradient in local Gaussian curvature. While these results remain qualitativeand their explanations are still incomplete, several of the observed phenomena are inqualitative agreement with experiments on active nematic liquid crystals on toroidalsurfaces and contribute to an understanding of the delicate interplay between geometricalproperties of the surface and characteristics of the flow field, which has the potential tocontrol active flows on surfaces via gradients in the spatial curvature of the surface.
Key words: active turbulence, generalised Navier-Stokes equation, toroidal surface
1. Introduction
To model fluids on curved surfaces is a problem which dates back to Scriven (1960),who derived a covariant Navier-Stokes (NS) equation and established the couplingbetween spatial curvature and fluid flow. The influence of geometric properties on both,equilibrium configurations and the dynamics far from it, has consequences in a hugevariety of problems, ranging from planetary flows (Delplace et al. et al. et al. et al. et al. et al. et al. † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b M. Rank and A. Voigt (Pearce et al. et al. et al. et al. (2013) for a review - which provides the mathematical basis for the vorticity-stream function formulation, not only splits the tangential velocity field into curl-free anddivergence-free components, but might also contain non-trivial harmonic vector fields -vector fields which are curl- and divergence-free. Their richness depends on the topologyof the surface (Jost 1991). As these vector fields cannot be described by the vorticity-stream function formulation, the approach is only applicable for surfaces, where harmonicvector fields are trivial, which are only simply-connected surfaces (Nitschke et al. b ).Reuther & Voigt (2018) and Fries (2018) introduce a numerical approach for a surfaceNS equation in its velocity-pressure formulation, which is applicable on general surfaces.The underlying idea was independently introduced in Nestler et al. (2018), Jankuhn et al. (2018) and Hansbo et al. (2020) and is described in detail in Nestler et al. (2019). It relieson an extension to an embedding thin film. This allows to express the covariant derivativesin terms of partial derivatives along the Euclidian basis. The surface partial differentialequation can thus be solved by considering each component by established methods forscalar quantities and enforcing the normal contributions to be zero, either by penalisationor by Lagrange multipliers. Also this approach has been applied to active flows in thecontext of surface active polar gels (Nitschke et al. et al. (2017), Torres-Sanchez et al. (2020), Sahu et al. (2020) and Lederer et al. (2020).We will here consider a modelling approach for active flows and extend studies for ageneralised Navier-Stokes (GNS) equation on a sphere (Mickelin et al. et al. et al. et al. b ). The model serves as a minimal model for active fluids andhad been successful in reproducing active turbulence flow patterns of swimming bacteria,ATP-driven microtubules and artificial microswimmers, see, e.g. (Bratanov et al. et al. et al. et al. et al. et al. ctive flows on curved surfaces et al. et al. et al. et al. ± , are attracted by regions of like-sign Gaussian curvature. A phenomenonwhich has been explored in detail for disclinations in positional order (Bowick et al. et al. et al. et al. (2018) suggest that activeflows on surfaces can be guided and controlled via gradients in the spatial curvature ofthe surface. They show that pairs of defects unbind and segregate in regions of oppositeGaussian curvature. At least qualitatively these results have been reproduced in Pearce et al. (2019), showing a linear dependency of defect creation and annihilation rates onGaussian curvature. As these rates are directly linked to active turbulence (Giomi 2015),a connection between spatial curvature and active turbulence can be expected. Pearce et al. (2019) consider an active nematodynamics model which is based on a simplifiedsurface Landau-de Gennes energy (Kralj et al. et al. (2017), Nitschke et al. (2018) and Nestler et al. (2020). As already pointedout, the numerical treatment in Pearce et al. (2019) is based on a vorticity-stream functionformulations, which is inappropriate for toroidal surfaces. Appropriate numerical studiesfor active nemotodynamics on general surfaces are still under development. However,several open questions, e.g., the influence of spatial curvature on active turbulence canalready be studies using an appropriately coarse-grained mesoscopic GNS equation.Here we focus on one aspect of the active flow. We only consider the influence of spatialvarying curvature on anomalous vortex-network active turbulence. The simulations inMickelin et al. (2018) on a sphere reveal a global curvature-induced transition froma quasi-stationary burst phase to an anomalous vortex-network turbulent phase anda classical 2D Kolmogorov turbulent phase. The new type of anomalous turbulenceis characterised by the self-assembly of finite-size vortices into linked chains of anti-ferromagnetic order, which percolate through the entire fluid domain. The coherentmotion of this vortex chain network provides an upward energy transfer and thus analternative to the conventional energy cascade in classical 2D hydrodynamic turbulence.We will answer the question if this mechanism is altered by gradients in curvature byconsidering the GNS equation on different tori within a parameter setting which leadsto anomalous vortex-network turbulence on a sphere.The paper is organised as follows: In Section 2 we introduce a covariant formulationof the GNS equation which is applicable to general curved surfaces. In Section 3 wedescribe the numerical approach, including basic validations for NS and GNS equations.The discussion on the influence of curvature on anomalous vertex-network turbulence isdone in Section 4 and conclusions are drawn in Section 5.
2. Mathematical model
In Mickelin et al. (2018) the surface generalised Navier-Stokes (GNS) equation ∂ t u ( x , t ) + ∇ u u ( x , t ) = − grad M p ( x , t ) + div M T , (2.1)div M u ( x , t ) = 0 (2.2) M. Rank and A. Voigt with initial condition u ( x ,
0) = u ( x ) ∈ T x M was proposed on M × (0 , ∞ ) with M asphere. The tangential fluid velocity at point x ∈ M and time t ∈ (0 , ∞ ) is denoted by u ( x , t ) ∈ T M and the surface pressure by p ( x , t ) ∈ R . Mickelin et al. (2018) considers thesurface tension − p ( x , t ). ∂ t is the time derivative, ∇ u the covariant directional derivative,grad M the surface gradient, and div M and div M the surface divergence, for vector- andtensor-fields, respectively. The stress tensor T contains passive and active contributionsand reads T = f ( ∆ M ) (cid:16) grad M u + (grad M u ) T (cid:17) , (2.3) f ( ∆ M ) = Γ − Γ ∆ M + Γ ∆ M , (2.4)where ∆ M = ∆ M ∆ M with a surface Laplacian ∆ M and real parameters Γ , Γ , Γ ∈ R .The constants Γ and Γ are assumed to be positive to ensure asymptotic stability,whereas Γ may have either sign. For Γ < f ( ∆ M ) as an effective viscosity, sufficiently negative Γ may turnthis quantity negative and thus providing a source of energy which makes the modeleffectively active. We can think about ∆ M as a wildcard for any surface Laplacian.Mickelin et al. (2018) have been using the Bochner Laplace operator ∆ M = ∆ B . Otherchoices are the Laplace-deRham operator ∆ dR or the Q-Laplacian ∆ Q . They are relatedby ∆ Q u = ∆ B u + κ u = ∆ dR u + 2 κ u , with κ the Gaussian curvature (Abraham et al. κ is constant, evaluating div M T is not an issue for any of the operators.However, evaluating div M T on more general surfaces, like the torus, turns out to bemost convenient using the Q-Laplacian. This leads to the following computation: div M T = div M f ( ∆ Q ) (grad M u + (grad M u ) T )= div M f ( ∆ Q ) · ∇ Q u = 2 Γ div M ∇ Q u − Γ div M ∆ Q ∇ Q u + 2 Γ div M ∆ Q ∇ Q u = Γ · div M ∇ Q u − Γ · div M ∇ Q ∆ Q u + Γ · div M ∇ Q ∆ Q u = Γ ∆ Q u − Γ ∆ Q u + Γ ∆ Q u , (2.5)where ∆ Q = ∆ Q ∆ Q ∆ Q . Introducing the auxiliary quantities v = ∆ Q u and w = ∆ Q v = ∆ Q u , using eq. (2.5) and the Bochner Laplacian ∆ B , eqs. (2.1) and (2.2) canbe rewritten as a system of second order partial differential equations ∂ t u + ∇ u u = − grad M p + Γ ( ∆ B u + κ u ) − Γ ( ∆ B v + κ v ) + Γ ( ∆ B w + κ w ) , (2.6) v = ∆ B u + κ u , (2.7) w = ∆ B v + κ v , (2.8)div M u = 0 , (2.9)which we will consider as a surface GNS equation on M , with M now a compact smooth2-manifold without boundary. In the case of constant Gaussian curvature κ the modelcoincides with the model presented in Mickelin et al. (2018) with the particular choice of Γ = Γ − Γ κ + 4 Γ κ , Γ = Γ − Γ κ and Γ = Γ . Thereby Γ , Γ and Γ denote theaccording parameters introduced in Mickelin et al. (2018). For Γ = Γ = 0 we obtain asa special case the surface NS equation (Scriven 1960). This limit, which also follows asa thin-film limit of the 3D NS equation with Navier boundary conditions (Miura 2018;Nitschke et al. et al. ctive flows on curved surfaces f ( ∆ Q ), whichaccounts for the intrinsic solvent fluid viscosity and contributions representing the stressesexerted by the active components on the fluid. For Γ < κ exact stationary solutions can beconstructed, which show some aspects of this influence (see Mickelin et al.
3. Numerical approach
Following the general approach (Nestler et al. u ( x , t ) ∈ T M only approximately. This system can thenbe numerically solved using a time discretisation based on a Chorin projection methodand a space discretisation by a regular surface triangulation and a scalar-valued surfacefinite element method (Dziuk & Elliott 2013) applied to each component of the extendedvelocity field, each component of the extended auxiliary variables and the pressure field.The approach extends the surface finite element discretisation for the surface NS equation(Reuther & Voigt 2018) to the surface GNS equation.3.1.
Reformulation
Instead of the tangential fields u ( x , t ) , v ( x , t ) , w ( x , t ) ∈ T M in the surface GNSequation (2.6) - (2.9) we will consider R -valued vector fields ˆ u ( x , t ) , ˆ v ( x , t ) , ˆ w ( x , t ) ∈ R which are only weakly tangential to T M . We thereby approximate the surface GNSequation following the general method proposed by Nestler et al. (2019). We considera neighbourhood U δ of M with coordinate projection π : U δ → M of ˜ x = a (˜ x ) + d (˜ x ) ν ( a (˜ x )) ∈ U δ , with d the signed distance function and ν the surface normal, givenby ˜ x (cid:55)→ x = a (˜ x ) ∈ M . Depending on the curvature of the surface M , this coordinateprojection is injective for δ > U δ such that ˜ u (˜ x ) = u ( x ), ˜ v (˜ x ) = v ( x ),˜ w (˜ x ) = w ( x ) and ˜ p (˜ x ) = p ( x ) are obtained for all ˜ x ∈ U δ . We extend the tangentialdifferential operators byˆ ∇ c ˆ u = π M · ( ∇ ˆ u − ∇ ˆ u · νν ) , (cid:99) div M ˆ u = div (ˆ u ) − ν · ( ∇ ˆ u · ν ) , ˆ ∆ B ˆ u = (cid:100) div M ˆ ∇ c ˆ u , where ∇ , div and div denote the common vector gradient and divergence operatorsin R . For the general form for (cid:100) div M we refer to Nestler et al. (2019)(Eq. E.2). Thepointwise normal projection is given by π M = I − νν T , where I denotes the identitymatrix. Adding α ( ν · ˆ u ) ν , α ( ν · ˆ v ) ν and α ( ν · ˆ w ) ν with penalty parameter α ∈ R bigenough penalises the normal components of the vectors ˆ u , ˆ v and ˆ w , see Nestler et al. (2019). This motivates the choice of the operators above, as we assume the velocity ˆ u to be approximately tangential to the surface. In Nestler et al. (2019) it was shownthat ˆ ∆ B ˆ u ≈ ∆ B u , when using an appropriate penalty term. For convergence studiesand possible dependencies of α on mesh size, we refer to Hansbo et al. (2020). Reuther& Voigt (2018) motivate the replacement of the surface divergence div M u by (cid:99) div M ˆ u .This leads to the following approximation of the surface GNS equation in Cartesian M. Rank and A. Voigt coordinates: ∂ t ˆ u + ˆ u · ˆ ∇ c ˆ u = − ˆ ∇ c p + Γ ˆ v − Γ ˆ w + Γ ( ˆ ∆ B ˆ w + κ ˆ w ) − α ( ν · ˆ u ) ν , (3.1)ˆ v = ˆ ∆ B ˆ u + κ ˆ u − α ( ν · ˆ v ) ν , (3.2)ˆ w = ˆ ∆ B ˆ v + κ ˆ v − α ( ν · ˆ w ) ν , (3.3) (cid:99) div M ˆ u = 0 . (3.4)The above model formulation coincides with the initial one just in the case of ν · ˆ u = ν · ˆ v = ν · ˆ w = 0 and ensures only a weak form of tangency of the according solution ˆ u (Reuther & Voigt 2018; Nestler et al. u , ˆ v , and ˆ w to be tangential to the surface, which is legitimate for appropriate α (Reuther &Voigt 2018; Hansbo et al. ∧ -sign of ˆ u , ˆ v , ˆ w ,ˆ ∇ c , (cid:99) div M and ˆ ∆ B will be omitted in the following.3.2. Time discretisation
The procedure applied for the discretisation of time is based on Chorin (1968) andRannacher (1992) and was successfully applied by Reuther & Voigt (2018) for the surfaceNS equation. For numerical reasons, we linearise the nonlinear covariant directionalderivative and obtain u ∗ · ∇ c u ∗ ≈ u n · ∇ c u ∗ . This yields the following problem with asemi-implicit Euler time scheme:Let u := u ∈ T M be the sufficiently smooth initial velocity field and τ n ∈ R be thetime step in the n -th iteration. For n → n + 1 determine successively(i) u ∗ , v ∗ and w ∗ such that1 τ n ( u ∗ − u n ) + ( u n · ∇ c ) u ∗ − Γ v ∗ + Γ w ∗ − Γ ( ∆ B w ∗ + κ w ∗ ) + α ( ν · u ∗ ) ν = 0 , (3.5) v ∗ − ∆ B u ∗ − κ u ∗ + α ( ν · v ∗ ) ν = 0 , (3.6) w ∗ − ∆ B v ∗ − κ v ∗ + α ( ν · w ∗ ) ν = 0 . (3.7)(ii) p n +1 such that τ n ∆ M p n +1 = div M u ∗ , (3.8)with ∆ M the Laplace-Beltrami operator.(iii) u n +1 such that u n +1 = u ∗ − τ n · ∇ c p n +1 . (3.9)3.3. Space discretisation
A regular triangulation M h = (cid:83) T ∈T T of the smooth surface M is constructed bytriangular elements T ∈ T determined by fixed points that are distributed equally overthe surface. Note that the particular choice of M h in general also affects the normalvector field. If analytically known, we define the unit normal vector field ν h of M h asthe analytic normal of M in each degree of freedom (DOF). Otherwise, in order toachieve convergence, an approximation is needed which is at least one order better thanthe approximation of the surface (Hansbo et al. M h and computing the normalsfrom them, as, e.g., considered in Reuther et al. (2020); Nitschke et al. (2020 a ).The surface differential operators are manipulated to operate on M h instead of M essentially by using the pointwise normal projection onto M h , which is given by π M h = I − ν h ν Th . We apply a scalar-valued surface finite element method for each component ctive flows on curved surfaces V h = (cid:8) ϕ ∈ L ( M h ) : ϕ | T ∈ P for all T ∈ T (cid:9) for trial and test functions. For k = 1 , , u ∗ k , v ∗ k , w ∗ k , u nk , v nk and w nk be the sufficiently smooth k -th component of u ∗ , v ∗ , w ∗ , u n , v n and w n , respectively. We multiply each component of equations (3.5)-(3.9) witha smooth test function ϕ ∈ V h , integrate over the domain M h and apply the divergencetheorem to achieve the weak formulation. This yields the following time and space discreteproblem:Let u = (cid:0) u k (cid:1) k with u k ∈ V h for all k = 1 , , n → n +1determine successively(i) u ∗ k , v ∗ k , w ∗ k ∈ V h such that1 τ n (cid:90) M h u ∗ k · ϕ d A h + (cid:90) M h u n ∇ c u ∗ k · ϕ d A h − Γ (cid:90) M h v ∗ k · ϕ d A h + Γ (cid:90) M h w ∗ k · ϕ d A h + Γ (cid:90) M h ∇ c w ∗ k · ∇ c ϕ d A h − Γ (cid:90) M h κw ∗ k · ϕ d A h + α (cid:88) i =1 (cid:90) M h ν i u ∗ i ν k · ϕ d A h = 1 τ n (cid:90) M h u nk · ϕ d A h , (3.10) (cid:90) M h v ∗ k · ϕ d A h + (cid:90) M h ∇ c u ∗ k · ∇ c ϕ d A h − (cid:90) M h κu ∗ k · ϕ d A h + α (cid:88) i =1 (cid:90) M h ν i v ∗ i ν k · ϕ d A h = 0 , (3.11) (cid:90) M h w ∗ k · ϕ d A h + (cid:90) M h ∇ c v ∗ k · ∇ c ϕ d A h − (cid:90) M h κv ∗ k · ϕ d A h + α (cid:88) i =1 (cid:90) M h ν i w ∗ i ν k · ϕ d A h = 0 (3.12)for every test function ϕ ∈ V h and for all k = 1 , ,
3. Thereby, u n = ( u nk ) k with u nk ∈ V h for all k = 1 , , n .(ii) p n +1 ∈ V h such that τ n (cid:90) M h ∇ c p n +1 · ∇ c ϕ d A h + (cid:90) M h div M u ∗ · ϕ d A h = 0 (3.13)for every test function ϕ ∈ V h .(iii) u n +1 k ∈ V h such that (cid:90) M h u n +1 k · ϕ d A h = (cid:90) M h u ∗ k · ϕ d A h − τ n (cid:90) M h (cid:2) ∇ c p n +1 (cid:3) k · ϕ d A h . (3.14)for every test function ϕ ∈ V h and for all k = 1 , , et al. (cid:96) ) (Van der Vorst 1992). M. Rank and A. Voigt
Figure 1: Parameterisation of torus with inner and outer radius r and R = 1 /r , andcorresponding angles θ and θ . 3.4. Validation
We validate the numerical approach against known results for passive flows on a torusand active flows on a sphere. In the following, all presented numbers and values aredimensionless quantities.3.4.1.
Passive flows on torus
We first consider passive flows on a torus. This has been done in detail by Reuther& Voigt (2018), utilising a surface NS equation which is achieved by Γ = 1 / Re and Γ = Γ = 0. We consider the time step τ n = 0 . α = 3000.The used mesh of the torus surface with outer radius R = 2 . r = 0 . u , we consider the arithmetic mean of the two linear independent harmonic vector fields, u harmθ = 12 (cid:107) x (cid:107) ∂ θ x , u harmθ = 14 (cid:107) x (cid:107) ∂ θ x , with θ , θ ∈ [0 , π ], see Figure 1 for a definition of the parameters. We normalise thetotal kinetic energy of the system E ( t ) = (cid:82) M h (cid:107) u ( t ) (cid:107) d A by dividing by the maximumvalue E max = E (0) and compare it with the benchmark problem in Reuther & Voigt(2018), see Figure 2(a). The agreement gives a first validation of the numerical approach.Figure 2(b) shows the normalised kinetic energy over time for different choices of theinner radius r and according outer radius R = 1 /r to preserve the surface area of thetorus. Already for this case a dependency of the dynamics on geometric properties isrealised. This includes not only a faster decay for increasing inner radius r , but also afaster decay in regions of negative Gaussian curvature, see Figure 2(a).3.4.2. Active flows on sphere
We next compare active flows on a sphere with results in Mickelin et al. (2018). Weconsider parameters leading to anomalous vortex-network turbulence, compare case (d)in Figure 2 of Mickelin et al. (2018). The penalty parameter α = 3 ,
000 is as in the passivecase, but the time step is reduced to τ n = 0 . R = 1. Figure 3(a) shows a snapshot of the normalised ctive flows on curved surfaces (a) (b) Figure 2: Normalised total kinetic energy of the implemented model with parameters Γ = 0 . Γ , Γ = 0, penalty parameter α = 3000, time step τ n = 0 . r = 0 . R = 2 . r = 0 . , . , ..., .
9, with R = 1 /r accordingly, on the entiredomain, against time. The distinction between regions of positive and negative Gaussiancurvature for different r leads to the same behaviour as shown in (a). (a) (b) (c)(d) Figure 3: Snapshots of (a) vorticity and (b) surface tension field at time t = 860. (c)Normalised kinetic energy - total kinetic energy divided by mean kinetic energy afterrelaxation, indicated by vertical grey line at t = 300 - against time in comparison withdata from Mickelin et al. (2018). Considered parameters are Γ = 3 . · − , Γ = − . · − and Γ = 3 . · − . Different line types emerge from different randominitial conditions. (d) Distribution of frequencies found by FFT averaged over threesimulation runs, the strength is given by the absolute value of the according coefficientof the FFT. Black lines indicate the according sample standard deviation.vorticity field. The vorticity is computed by φ = curl M u , with curl M the surface rotation,defined by curl M u = div M ( u × ν ) and u to be interpreted as the tangential velocityfield supplemented with a third component in normal direction set to be zero. φ isa scalar field, the magnitude of the vector pointing in normal direction. Figure 3(b)shows the corresponding normalised surface pressure/surface tension. To be comparablethe considered normalisation follows Mickelin et al. (2018). Anomalous vortex-networkturbulence is characterised in Mickelin et al. (2018) by topological measures of the0 M. Rank and A. Voigt
Figure 4: Convergence properties for vorticity, pressure/tension and kinetic energy. The L error is computed with respect to the solution on the finest mesh and shown withrespect to the number of degrees of freedom (DoFs). The black and the red lines indicateorder 1 and 1 /
2, respectively.vorticity and geometric measures of the high-tension domains. The last shows highlybranched chain-like structures which are clearly visible in Figure 3(b). This provides afurther qualitative validation of the numerical approach. Finite-size vortices self-organiseinto chain complexes of antiferromagnetic order that percolate trough the entire surfaceforming an active dynamic network, which provides an efficient mechanism for upwardenergy transport from smaller to larger scales. The dynamics of this process is validatedagainst Mickelin et al. (2018) in Figure 3(c) by comparing the normalised kinetic energy,the total kinetic energy divided by the mean kinetic energy after relaxation. As the initialdata in Mickelin et al. (2018) is not known, we consider different simulations with differentinitial data. The relaxation time, as well as the amplitude and frequency of the observedoscillations are independent on the initial data and in reasonable agreement with Mickelin et al. (2018). To compare the amplitudes of the energy curves, we have computed the rootmean square RMS = ( T − T (cid:82) T T (cid:0) E ( t ) − ¯ E (cid:1) d t ) with T = 300 and T = 1600, andmean kinetic energy ¯ E . The computed values are RMS = 0 . , .
109 and 0 .
070 for oursimulations and RMS = 0 .
087 for the data of Mickelin et al. (2018) with peak-to-peakamplitude values of 0 . , . , .
333 and 0 . et al. (2018).To demonstrate the appropriate choice of the numerical parameters we consider differ-ent mesh resolutions. Mesh refinement is done by bisection with new vertices projected tothe spherical surface. To circumvent relaxation effects in the initialisation phase we startwith initial conditions, obtained from the previous simulations after t = 300, which areinterpolated to the new meshes. The solutions are almost indistinguishable by eye fromthe one plotted in Figure 3(a) and (b). Figure 4 shows the L -errors in the vorticity φ andsurface pressure / surface tension p , as well es differences in the kinetic energy E after400 time steps with constant time step τ n . The mesh used for the simulations shown inFigure 3(a) and (b) corresponds to the second data point in Figure 4. The results indicateconvergence of order 1 / ctive flows on curved surfaces
4. Results
Measures for topological and geometric quantities
In analogy to Mickelin et al. (2018) we determine the topology of the vorticity fieldsand the geometry of the high-tension domains to identify anomalous turbulence. Wetherefore define the normalised Betti number asBetti φ ( r ) = (cid:104) N φ ( r, t ) − N sphere ( t ) (cid:105)(cid:104) N sphere ( t ) (cid:105) , with N sphere ( t ) denoting the zeroth Betti number of a sphere with radius R = 1 at time t as a reference value. The zeroth Betti number N φ ( r, t ), with r the inner radius ofthe considered torus, measures the number of connected domains with a high absolutevorticity { x ∈ M : φ ( x , t ) > α φ · max x ∈ M φ ( x , t ) or φ ( x , t ) < α φ · min x ∈ M φ ( x , t ) } anda threshold α φ > t . The time average (cid:104)·(cid:105) is taken after the initial relaxationperiod. Intuitively, large values of Betti φ indicate many vortices of comparable circulationor many connected structures of similar size, whereas small values suggest the presence ofa few dominant eddies or large connected structures. Accordingly, the normalised Branchnumber is defined by Branch p ( r ) = (cid:104) A p ( r, t ) − A sphere ( t ) (cid:105)(cid:104) A sphere ( t ) (cid:105) , with A sphere ( t ) serving as a reference value. The Branch number A p ( r, t ) denotes the meanof the ratios ∂A/A , where A denotes the area and ∂A denotes the boundary length ofeach connected component of the regions { x ∈ M : p ( x , t ) > β p ¯ p ( t ) } with ¯ p its meanvalue and β p > ∂A/A is a measure of chainlike structuresin the surface pressure/surface tension fields, a large value signaling a highly branchedstructure, whereas smaller values indicate less branching. N sphere ( t ) and A sphere ( t ) areconsidered to compare with anomalous turbulence in the spherical case in Mickelin et al. (2018). 4.2. Active flows on toroidal surface
We are using a mesh with 48,832 vertices and 97,664 resulting triangle elements, whichwe scale to approximate all tori with inner radius r and according outer radius R = 1 /r .The total simulation time is T = 1600, which turns out to be sufficient to reveal thecharacteristic active dynamics of the system.The separation of the general model dynamics into a relaxation phase and an activechaotic phase can be done in the same way as already described for the sphere with gener-ating and dissolving similar vortex patterns for all choices of inner radii r . Figure 5(a),(b)show snapshots of the normalised vorticity and the nomalised surface pressure/surfacetension of a simulation on a torus. As for the sphere, the relaxation phase also runsuntil about t = 300. Thereafter the active regime begins. Comparing the still imagesof the torus and the sphere, the vorticity seems to have the tendency to form slightlymore complex, longer and maze-like structures on the torus. The tension still has thecharacteristic branched chain-like structures. However, the images give an impression onprefered orientations of the structure. Nevertheless, the patterns and dynamics appearto be qualitatively similar to the ones on a sphere. To investigate the differences between2 M. Rank and A. Voigt (a) (b) (c)(d)
Figure 5: Snapshots of (a) vorticity and (b) surface tension field on torus with radii r = 0 . R = 2 .
0. Other parameters are as in Figure 3. (c) Normalised kinetic energy- total kinetic energy for different choices of radii r = 0 . , . , . R = 1 /r dividedby average total kinetic energy over all tori after relaxation - against time. The verticalgrey line at t = 300 marks the beginning of the active turbulence regime. The resultsfor a sphere are shown for comparison. (d) Distribution of frequencies found by FFT,averaged over three simulation runs, the strength is given by the mean absolute value ofthe according coefficient of the FFT. Black lines indicate the according sample standarddeviation. Inner radius r Table 1: Average normalised kinetic energy, see Figure 5(c), peak-to-peak amplitude androot mean square (RMS) values of normalised kinetic energy curves, averaged over timeand simulation runs for tori with inner radius r and sphere with R = 1.the simulation results on the sphere and the toroidal surfaces for different radii r and R = 1 /r we compute the normalised kinetic energy, see Figure 5(c), and the averageRMS values over three simulation runs, see Table 1. The kinetic energy progressionsof the different tori are similar to each other and to the sphere. Nevertheless, thin toriwith inner radius r = 0 . , . r = 0 . , .
9. This observation is underlined by the peak-to-peak amplitude values.The distribution of frequencies in the FFT of the tori is shown in Figure 5(d). We candetect similar dominant frequencies for all tori with no significant differences between thevarious torus dimensions. Depending on the random initial conditions, distinct frequencydistributions emerge, see Figure 5(d).To explain the observed differences in the energy fluctuations requires a deeper analysisof the data. The normalised Betti and Branch numbers for the whole surface have beencomputed for tori with inner radius r = 0 . , . , ..., . ctive flows on curved surfaces Inner radius r φ ( r ) per area -0.436 -0.161 -0.166 -0.067 -0.055 -0.182 -0.568Branch p ( r ) per area 2.438 1.555 1.291 1.563 1.673 0.820 0.733 Table 2: Normalised Betti and Branch numbers per area for the entire surface averagedover simulation runs for tori with inner radius r . (a) (b) (c) (d) Figure 6: Gaussian curvature κ from blue to red on tori with inner radius a) r = 0 . r = 0 .
5, c) r = 0 .
7, d) r = 0 . R = 1 /r . r does not scalebetween (a) - (d). R = 1 /r with parameters α φ = 0 . β p = 1. They are divided by area and theiraveraged values, over different simulation runs, are given in Table 2. We can observelower normalised Betti numbers for very thin ( r = 0 .
3) and very thick ( r = 0 .
9) tori.The normalised Branch numbers show a different trend. They are largest for very thin( r = 0 .
3) and smallest for thick ( r = 0 . , .
9) tori.However, averaging over the entire surface does not account for the local differences insurface curvature. To extract a dependency of the normalised Betti and Branch numberson local curvature we classify regions of equal Gaussian curvature on all consideredtori. Figure 6 shows the different values of Gaussian curvature κ on selected tori. Whilemoderate values of positive Gaussian curvature are found on the outer part for all tori,strong negative values are only present in the inner part of thick tori ( r = 0 . , . − et al. (2018). The behaviour for the normalisedBranch number is similar, with maximal values for Gaussian curvature between − . r = 0 .
3) and decrease towards very thick tori ( r = 0 . et al. M. Rank and A. Voigt (a) (b)(c) (d)
Figure 7: a) Normalised Betti numbers per area for surface areas with Gaussian curvature( κ − . , κ + 0 . κ ∈ {− . , − . , ..., . } and b) normalised Branch number for surfaceareas with Gaussian curvature ( κ − . , κ + 0 . κ ∈ {− . , − . , ..., . } vs. innertorus radius r , the row at the bottom showing the average value over all consideredtorus geometries. A few intervals of Gaussian curvature are not found on some tori,indicated by grey colour, see Figure 6. The average values with according sample standarddeviation of different simulation runs are shown in c) for the normalised Betti numberand values for one simulation in d) for the normalised Branch number. The black linesshow the according mean values over all runs and tori. The results are obtained withthe parameters α φ = 0 . β p = 1. Results for different parameters are shown in theAppendix and demonstrate the robustness of the results.also compute the enstrophy E ( t ) = (cid:82) M h φ dA for the whole surface and per region ofconstant Gaussian curvature, see Figure 8(a) and (b), respectively. The values are largestfor moderate tori ( r = 0 . , . , .
7) and decay towards thinner ( r = 0 . , .
4) and thicker( r = 0 . , .
9) tori, see also Table 3. In Mickelin et al. (2018) it is argued that also theratio between the mean kinetic energy and the mean enstrophy can be used to identifythe anomalous turbulence regime. Even if this can only be a qualitative measure theratio is shown in Table 3. It shows a minimum for moderate tori ( r = 0 . , . , .
7) andmoderate Gaussian curvature values ( κ ∈ [ − . , . et al. ctive flows on curved surfaces (a) (b) Figure 8: a) Enstrophy per surface area nomalised by mean after relaxation over all tori,shown for inner torus radius r = 0 . , . . . , .
9. (b) Normalised enstrophy for surface areaswith Gaussian curvature ( κ − . , κ + 0 . κ ∈ {− . , − . , ..., . } vs. inner torusradius r , the row at the bottom showing the average value over all considered torusgeometries. Inner radius r (cid:104) E ( t ) (cid:105) / (cid:104)E ( t ) (cid:105) κ -2.25 -1.75 -1.25 -0.75 -0.25 0.25 (cid:104) E ( t ) (cid:105) / (cid:104)E ( t ) (cid:105) Table 3: Mean normalised enstrophy per area over time and simulation runs for toriwith inner radius r . Ratio of mean normalized kinetic energy and mean normalizedenstrophy per area for tori with inner radius r as well as intervals of Gaussian curvature( κ − . , κ + 0 . r and R (see Bowick et al. et al. et al. (2018) below.We first address an observation in the movies in the Electronic Supplement. Theyindicate a preferred direction of the chained vortex structures. They seem to align withthe principle curvature lines and prefer the one with lower absolute value of Gaussiancurvature. For regions of positive Gaussian curvature (outer part) they align horizontally,whereas for negative Gaussian curvature (inner part), at least for strong values, thealignment is vertically. This observation is confirmed in Figure 9(b) showing the averagedirection of elongation of the chained vortex structures. The still image in Figure 9(a) canonly partly confirm this, we therefore refer to the corresponding movies in the ElectronicSupplement. Figure 9(d) and (c), show the corresponding average direction of elongationof high-tension domains and a characteristic still image, respectively. Representativestill images for the other radii are provided in the Appendix. The average directionshave been computed using a method first introduced for quantifying deformations offoam structures by Asipauskas et al. (2003) and was later applied for the elongationof cellular structures (Mueller et al. M. Rank and A. Voigt
Inner radius r Table 4: Elongation eigenvalue averaged over time and simulation runs for tori with innerradius r .the segmented structures with a smooth transition over a small interface. All informationabout the elongation can be represented in terms of the gradient of these phase fields. Theelongation is described as the angle against the lines of constant Gaussian curvature, i.e.the green lines in Figure 1. Thus, an elongation of zero represents a horizontal structurealigned with the green lines and a values of π/ r = 0 . , .
4) and thick tori ( r = 0 . , .
9) this preferred alignment ofthe chained vortex structures becomes evident. The elongation of the tension fields is lessevident. Only for thin tori ( r = 0 . , .
4) a preferred horizontal elongation is observed. Apreferred vertical elongation for thick tori ( r = 0 . , .
9) can not be seen. This might resultfrom less dominant branched structures in regions of strong negative Gaussian curvature,see Figure 7(b),(d). For moderate tori ( r = 0 . , . , .
7) the average elongation directionof the tension fields fluctuates much stronger with no preferred orientation. Similaralignment effects with minimal curvature lines have also been reported for surface liquidcrystals (Segatti et al. et al. et al. r = 0 . , .
4) and thick tori ( r = 0 . , . r = 0 . , . , . et al. et al. r = 0 . , . r = 0 . , .
9) results from the different magnitudes of the vorticity regionsthat are counted by the Betti number. Figure 10(b) shows the average total kinetic energyfor different regions of the tori. We can observe significantly lower values for regions ofstrong negative Gaussian curvature. This results in lower normalised Betti numbers for ctive flows on curved surfaces (a) (b)(c) (d) Figure 9: Average direction of elongation of chained vorticity (b) and branched high-tension structures (b) on different tori over time. The elongation direction is consideredwith respect to the parameterisation. Snapshots of the considered vorticity and surfacetension structures for a torus with inner radius r = 0 . (a) (b) Figure 10: (a) Normalised Betti numbers per area vs. inner torus radius for surface areaswith positive and negative Gaussian curvature obtained with parameter α φ = 0 .
5. Dottedlines show the means over all tori and dashed lines show the according means over all toriwith r ∈ { . , . , . , . , . } . (b) Average total kinetic energy for surface areas withGaussian curvature ( κ − . , κ + 0 . κ ∈ {− . , − . , ..., . } divided by overall meantotal kinetic energy of each torus.8 M. Rank and A. Voigt thick tori ( r = 0 . , . r = 0 . , . Γ κ u term in eq.(2.6), which, if neglecting all other contributions, leads to exponential decay for κ < κ >
0. It has to remain speculative if this simple explanationcan also be used for the highly nonlinear surface GNS equation. Another interestingaspect of Figure 10(b) are the values at κ = 0. The differences for different tori indicatenot only a dependency on the local Gaussian curvature but also on its gradient, whichis largest for thin tori r = 0 . r . Again this behaviour isin accordance to the experiments in Ellis et al. (2018).
5. Conclusions
We consider a surface generalised Navier-Stokes (GNS) equation as a minimal model foractive flows on arbitrarily curved surfaces. This extends work of Mickelin et al. (2018),who considered this equation on a sphere. The numerical approach extends work ofReuther & Voigt (2018) for the surface Navier-Stokes (NS) equation and is based on thegeneral concept to solve surface vector-valued partial differential equations on arbitrarysurfaces by surface finite elements (Nestler et al. et al. et al. r or the local Gaussiancurvature κ , and shows the same weak dependency on curvature as the consideredtopological and geometric measures.While a full understanding of the relation of Gaussian curvature on active turbulencerequires many further investigations, the simulations results indicate a clear dependencyof various aspects on the Gaussian curvature of the surface. Some are in qualitativeaccordance with the experiments on active nematic liquid crystals, which are constrainedto lie on a toroidal surface, see Ellis et al. (2018). These are larger numbers of vorticesin regions of negative Gaussian curvature and a dependency on the gradient in Gaussiancurvature on the kinetic energy. Other effects, such as the alignment of the chainedstructures with minimal curvature lines, ask for experimental validation. Anotherinteresting question for future research is the influence of local Gaussian curvature onthe transitions to classical 2D Kolmogorov turbulence. This transition is addressed inflat space and on the sphere (see Bratanov et al. et al. et al.ctive flows on curved surfaces (a) α φ = 0 . α φ = 0 . α φ = 0 . α φ = 0 . Figure 11: Normalised Betti numbers per area as in Figure 7(a) for differentthresholds to define regions with high absolute vorticity { x ∈ M : φ ( x , t ) > α φ · max x ∈ M φ ( x , t ) or φ ( x , t ) < α φ · min x ∈ M φ ( x , t ) } , with (a) α φ = 0 .
4, (b) α φ = 0 . α φ = 0 . α φ = 0 . et al. b , a ).Acknowledgment: This research was supported by the German Research Foundation(DFG) within the Research Unit 3013. We used computing resources provided by ZIH atTU Dresden and by J¨ulich Supercomputing Centre within project HDR06. We furtheracknowledge the provided data from Mickelin et al. (2018) by O. Mickelin and J. Dunkel,as well as support from M. Nestler, M. Salvalaglio and D. Wenzel concerning thepostprocessing. Appendix A.
Figures 11 and 12 provide the same information as Figure 10 for different thresholds α φ and β p , demonstrating the robustness of the results on these values. Figure 13 and14 show the corresponding still images for the vorticity and tension structures in Figure9 for the other tori. For the corresponding videos and additional visualizations of thesurface velocity using LIC we refer to the Electronic Supplement. REFERENCESAbraham, R., Marsden, J.E. & Ratiu, T.
Manifolds, Tensor Analysis, andApplications . Springer New York.
Alaimo, F., Koehler, C. & Voigt, A.
Sci. Rep. , 5211. M. Rank and A. Voigt (a) β p = 0 . β p = 0 . β p = 1 . β p = 1 . Figure 12: Normalised Branch numbers per area as in Figure 7(b) for different regions { x ∈ M : p ( x , t ) > β p · ¯ p ( t ) } with (a) β p = 0 .
8, (b) β p = 0 .
9, (c) β p = 1 . β p a = 1 . (a) r = 0 . r = 0 . r = 0 . r = 0 . r = 0 . r = 0 . r = 0 . Figure 13: Snapshots of normalised vorticity fields on different torus surfaces in the activeregime.
Apaza, L. & Sandoval, M.
Soft Matter ,9928–9936. Asipauskas, M., Aubouy, M., Glazier, J., Graner, F. & Jiang, Y.
Granular Matter , 71–74. Beresnev, I.A. & Nikolaevskiy, V.N.
Physica D , 1–6. ctive flows on curved surfaces (a) r = 0 . r = 0 . r = 0 . r = 0 . r = 0 . r = 0 . r = 0 . Figure 14: Snapshots of normalised surface tension fields on different torus surfaces inthe active regime.
Bhatia, H., Norgard, G., Pascucci, V. & Bremer, P.-T.
IEEE Trans. Vis. Comp. Graphics , 1386–1404. Bowick, M., Nelson, D.R. & Travesset, A.
Phys. Rev. E , 041102. Bratanov, V., Jenko, F. & Frey, E.
Proc. Nat.Acad. Sci. (USA) , 15048–15053.
Chang, Y.-W., Fragkopoulos, A.A., Marquez, S.M., Kim, H.D., Angelini, T.E. &Fernandez-Nieves, A.
New J. Phys. , 033017. Chorin, A.J.
Mathematics ofComputation (104), 745–762. Delplace, P., Marston, J.B. & Venaille, A.
Science , 1075–1077.
Dziuk, G. & Elliott, C.M.
Acta Numer. ,289–396. Ellis, P.W., Pearce, D.J.G., Chang, Y.-W., Goldsztein, G., Giomi, L. & Fernandez-Nieves, A.
Nature Physics , 85–90. Fries, T.-P.
Int. J. Num. Meth. Fluids , 55–78. Giomi, L.
Phys. Rev. X ,031003. Giomi, L. & Bowick, M.J.
Europ. Phys. J.E , 275–296. Golovaty, D., M.J., Alberto & Sternberg, P.
J. Nonlin. Sci. , 1905–1932. Gross, B. & Atzberger, P.J.
J. Comput. Phys. , 663–689.
Hansbo, P., Larson, M.G. & Larsson, K.
IMA J. Num. Anal. , 1652–1701. Heidenreich, S., Dunkel, J., Klapp, S.H.L. & B¨ar, M.
Phys. Rev. E , 020601. Heisenberg, C.-P. & Bellaiche, Y.
Cell , 948–962.
James, M., Bos, W.J.T. & Wilczek, M.
Phys. Rev. Fluids , 061101. Jankuhn, T., Olshanskii, M.A. & Reusken, A. M. Rank and A. Voigt embedded surfaces: Modeling and variational formulations.
Interf. Free Bound. , 353–377. Jesenek, D., Kralj, S., Rosso, R. & Virga, E.G.
Soft Matter , 2434–2444. Jost, J.
Riemannian Geometry and Geometric Analysis . Springer.
Keber, F.C., Loiseau, E., Sanchez, T., DeCamp, S.J., Giomi, L., Bowick, M.J.,Marchetti, M.C., Dogic, Z. & Bausch, A.R.
Science , 1135–1139.
Kralj, S., Rosso, R. & Virga, E.G.
SoftMatter , 670–683. Lederer, P.L., Lehrenfeld, C. & Sch¨oberl, J.
Int. J. Num. Meth. Eng. ,2503–2533.
Linkmann, M., Boffetta, G., Marchetti, M.C. & Eckhardt, B.
Phys.Rev. Lett. , 214503.
Linkmann, M., Hohmann, M. & Eckhardt, B. a Non-universal transitions to two-dimensional turbulence.
J. Fluid Mech. , A18.
Linkmann, M., Marchetti, M.C., Boffetta, G. & Eckhardt, B. b Condensateformation and multiscale dynamics in two-dimensional active suspensions.
Phys. Rev.E , 022609.
Mayer, M., Depken, M., Bois, J.S., J¨ulicher, F. & Grill, S.W.
Nature , 617–621.
Mickelin, Oscar, S(cid:32)lomka, Jonasz, Burns, Keaton J., Lecoanet, Daniel, Vasil,Geoffrey M., Faria, Luiz M. & Dunkel, J¨orn
Phys. Rev. Lett. , 164503.
Miura, T.-H.
Quart. Appl. Math. , 215–251. Mueller, R., Yeomans, J. M. & Doostmohammadi, A.
Phys. Rev. Lett. , 048004.
Nestler, M., Nitschke, I., L¨owen, H. & Voigt, A.
Soft Matter , 4032–4042. Nestler, M., Nitschke, I., Praetorius, S. & Voigt, A.
J. Nonlin. Sci. , 147–191. Nestler, M., Nitschke, I. & Voigt, A.
J. Comput. Phys. , 48–61.
Nitschke, I., Nestler, M., Praetorius, S., L¨owen, H. & Voigt, A.
Proc. Roy. Soc. A , 20170686.
Nitschke, I., Reuther, S. & Voigt, A.
Transport Processes at Fluidic Interfaces (ed. Bothe, D andReusken, A), pp. 177–197.
Nitschke, I., Reuther, S. & Voigt, A.
Phys. Rev. Fluids , 044002. Nitschke, I., Reuther, S. & Voigt, A. a Liquid crystals on deformable surfaces.
Proc.Roy. Soc. A , 20200313.
Nitschke, I., Reuther, S. & Voigt, A. b Vorticity-stream function approaches areinappropriate to solve the surface Navier-Stokes equation on a torus.
Proc. Appl. Math.Mech. , e202000006. Nitschke, I., Voigt, A. & Wensch, J.
J. Fluid Mech. , 418–438.
Pearce, D.J.G.
New J. Phys. ,063051. Pearce, D.J.G., Ellis, P.W., Fernandez-Nieves, A. & Giomi, L.
Phys. Rev. Lett. , 168002.
Rannacher, R.
The Navier-Stokes Equations II — Theory and Numerical Methods .Springer, Berlin, Heidelberg. ctive flows on curved surfaces Reusken, A.
IMA J. Num. Anal. , 109–139. Reuther, S., Nitschke, I. & Voigt, A.
J. Fluid Mech. , R8.
Reuther, S. & Voigt, A.
Multiscale Model. Sim. , 632–643. Reuther, S. & Voigt, A.
Phys. Fluids , 012107. Sahu, A., Omar, Y.A.D., Sauer, R.A. & Mandadapu, K.K.
J. Comput. Phys. , 109253.
Scriven, L.E.
Chem. Eng. Science , 98–108. Segatti, A., Snarski, M. & Veneroni, M.
Math. Mod. Meth. Appl. Sci. , 1865. Sipos, O., Nagy, K., Di Leonardo, R. & Galajda, P.
Phys. Rev. Lett. , 258104.
Slomka, J. & Dunkel, J.
Europ. Phys. J.- Spec. Top. , 1349–1358.
Slomka, J. & Dunkel, J.
Phys. Rev. Fluids , 043102. Supekar, R., Heinonen, V., Burns, K.J. & Dunkel, J.
J. Fluid Mech. , A30.
Tan, T.H., Liu, J., Miller, P.W., Tekant, M., Dunkel, J. & Fakhri, N.
Nature Physics , 657–662. Toner, J & Tu, Y.H.
Phys.Rev. E , 4828–4858. Torres-Sanchez, A., Millan, D. & Arroyo, M.
J. Fluid Mech. , 218–271.
Torres-Sanchez, A., Santos-Olivan, D. & Arroyo, M.
J. Comput. Phys. , 109168.
Tribelsky, M.I. & Tsuboi, K.
Phys. Rev.Lett. , 1631–1634. Turner, A.M., Vitelli, V. & Nelson, D.R.
Rev. Mod.Phys. , 1301–1348. Vey, S. & Voigt, A.
Comput. Vis. Sci. , 57–67. Van der Vorst, H.
SIAM Journal on Scientific and StatisticalComputing , 631. Wensink, H.H., Dunkel, J., Heidenreich, S., Drescher, K., Goldstein, R.E., L¨owen,H. & Yeomans, J.M.
Proc. Nat. Acad. Sci.(USA) , 14308–14313.
Witkowski, T., Ling, S., Praetorius, S. & Voigt, A.
Adv. Comput. Math.41