Quantifying the checkerboard problem to reduce numerical dissipation

Jou nal o Compu a ional Physics 521 (2025) 113537
A ailable online 29 Oc obe 2024
0021-9991/© 2025 The Au ho s. Published by Else ie Inc. This is an open access a icle unde he CC BY license
(h p://c ea i ecommons.o g/licenses/by/4.0/).
Con en s lis s a ailable a ScienceDi ec
Jou nal o Compu a ional Physics
jou nal homepage: www.else ie .com/loca e/jcp
Quan i ying he checke boa d p oblem o educe nume ical
dissipa ion
J.A. Hopmana,∗, D. San osa, À. Alsal i-Baldelloua,b, J. Rigolaa, F.X. T iasa
aHea and Mass T ans e Technological Cen e , Technical Uni e si y o Ca alonia, ESEIAAT, c/Colom 11, Te assa, 08222, Ba celona, Spain
bTe mo Fluids SL, c/Magí Cole 8, Sabadell, 08204, Ba celona, Spain1
A R T I C L E I N F O A B S T R A C T
Keywo ds:
Checke boa ding
Colloca ed g ids
Conse a i e disc e isa ion
This wo k p o ides a comp ehensi e explo a ion o a ious me hods in sol ing incomp essible
flows using a p ojec ion me hod, and hei ela ion o he occu ence and managemen o
checke boa d oscilla ions. I employs an algeb aic symme y-p ese ing amewo k, cla i ying
he de i a ion and implemen a ion o disc e e ope a o s while also add essing he associa ed
nume ical e o s. The lack o a p ope defini ion o he checke boa d p oblem is add essed
by p oposing a physics-based coefficien . This coefficien , oo ed in he dispa i y be ween he
compac - and wide-s encil Laplacian ope a o s, is able o quan i y oscilla o y solu ion fields wi h
a physics-based, global, no malised, non-dimensional alue. The influence o mesh and ime-s ep
efinemen on he occu ence o checke boa ding is highligh ed. The e o e, single measu emen s
using his coefficien should be conside ed wi h cau ion, as he alue p esen s li le use wi hou
any con ex and can ei he sugges mesh efinemen o use o a diffe en sol e . In addi ion, an
example is gi en on how o employ his coefficien , by es ablishing a nega i e eedback be ween
he le el o checke boa ding and he inclusion o a p essu e p edic o , o dynamically balance he
checke boa ding and nume ical dissipa ion. This me hod is es ed o lamina and u bulen flows,
demons a ing i s capabili ies in ob aining his dynamical balance, wi hou equi ing use inpu .
The me hod is able o achie e low nume ical dissipa ion in absence o oscilla ions o diminish
oscilla ion on skew meshes, while i shows minimal loss in accu acy o a u bulen es case.
Despi e i s ad an ages, he me hod exhibi s a sligh dec ease in he second-o de ela ion be ween
ime-s ep size and p essu e e o , sugges ing ha o he eedback mechanisms could be o in e es .
1. In oduc ion
The checke boa d p oblem a ises due o he non- i ial p essu e- eloci y coupling o incomp essible New onian fluid flows. The
go e ning equa ions o such flows a e gi en by he momen um and con inui y equa ions:
𝜕𝑡𝐮+(𝐮⋅∇)𝐮=𝜈∇2𝐮−1
𝜌∇𝑝, (1)
∇⋅𝐮=0,(2)
* Co esponding au ho .
E-mail add ess: [email p o ec ed] (J.A. Hopman).
URL: h ps://gi hub.com/janneshopman (J.A. Hopman).
1www . e mofluids .com.
h ps://doi.o g/10.1016/j.jcp.2024.113537
Recei ed 13 June 2024; Recei ed in e ised o m 24 Sep embe 2024; Accep ed 24 Oc obe 2024
Jou nal o Compu a ional Physics 521 (2025) 113537
2
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Fig. 1. Schema ic d awing o how he Cen al Diffe ence Scheme leads o decoupling o nodes.
which, o h ee spa ial dimensions, p o ide only h ee independen equa ions and ou unknowns. Fo incomp essible flows, he
p essu e canno be ob ained om he equa ion o s a e, and o mula ing an equa ion o p essu e becomes non- i ial, a opic ha
has been widely s udied in he field o compu a ional fluid dynamics (CFD) [1–4]. Many me hods sol e his p oblem i e a i ely by
p edic ing a eloci y field and hen finding a p essu e field o which he g adien p ojec s i on o a di e gence- ee space, using he
Helmhol z-Hodge heo em [5–9]. I he disc e e g adien a node 𝑖is de i ed h ough cen al diffe encing, i s alue will only depend
on he p essu e a neighbou ing nodes 𝑗and no on he p essu e a node 𝑖i sel . Mo eo e , i done consis en ly, he disc e e Laplacian
ope a o de i ed om his g adien will be based on a wide s encil. In his s encil, node 𝑖is coupled o nodes 𝑘, which a e neighbou s
o nodes 𝑗, esul ing in a decoupling be ween node 𝑖and di ec ly neighbou ing nodes 𝑗, see Fig. 1. The use o his so-called wide-
s encil Laplacian in he Poisson equa ion and he a o emen ioned g adien in he eloci y co ec ion esul s in a decoupling o he
p essu e field be ween neighbou ing cells. This in u n can lead o non-physical oscilla ions and checke boa d-like pa e ns in he
p essu e field, also known as he checke boa d p oblem [1].
On s uc u ed Ca esian meshes, his p oblem can be ci cum en ed by using a s agge ed g id a angemen in which he eloci ies a
he cell aces a e coupled o a compac -s encil g adien o p essu e be ween di ec ly neighbou ing cell-cen e ed nodes [10]. Howe e ,
o many indus ial applica ions he use o CFD in ol es complex geome ies ha equi e uns uc u ed meshes. Ex ension o he
s agge ed g id me hod o such cases is no s aigh - o wa d and leads o complex co ec ion schemes and inc eased compu a ional
cos s [11,12].
One commonly used s a egy o add ess p essu e field oscilla ions, is he use o he weigh ed in e pola ion me hod (WIM) o
e alua e he eloci y a he cell aces, in con as o he di ec in e pola ion me hod (DIM). This me hod was o iginally in oduced
as he p essu e-weigh ed in e pola ion me hod and has allowed he wide-sp ead usage o he colloca ed g id a angemen [13,14].
Since i s in oduc ion his me hod has been ex ended in many ways, e.g. o accoun o decoupling caused by he use o small ime-
s eps ound in ansien flows [15–18]and o accoun o unde - elaxa ion ac o s [19–21]. Gene alisa ion and unifica ion o he
a o emen ioned p oblems and accompanying solu ions can be ound in mo e ecen wo ks [22–24], which offe a comp ehensible
o e iew o hese me hods. One una oidable consequence ha hese solu ions ha e in common is he in oduc ion o a nume ical
e o in he o m o a non-ze o con ibu ion o he kine ic ene gy balance o ei he he p essu e e m o he con ec i e e m [25].
Ne e heless, many comme cially a ailable and open-sou ce codes a ou s abili y a he cos o accu acy [26–29], by applying a
dissipa i e o m o he WIM. This a oids p essu e field oscilla ions while a he same ime inc easing s abili y and accep ing he
consequen ial nume ical e o . Usually his is done implici ly h ough he use o a compac -s encil Laplacian, which di ec ly couples
node 𝑖 o i s neighbou ing nodes 𝑗ins ead o i s second-neighbou s 𝑘. In his me hod, di e gence- ee alues o he eloci ies a he
aces ollow di ec ly om he Poisson equa ion, wi hou in oducing any nume ical e o . Howe e , he coupling effec o he WIM
is in oduced implici ly by he Laplace ope a o and he nume ical e o is in oduced o he colloca ed eloci ies in his case. This
me hod in oduces a dissipa i e p essu e e o , wi h he benefi o g ea ly educing compu a ional complexi y and cos when using
uns uc u ed g ids.
Howe e , wi h he inc ease in compu a ional esou ces a ailable and he ensuing ise o high-fideli y simula ions, highe accu acy
is desi ed o desc ibe mo ion o fluids. Nume ical dissipa ion is limi ing his accu acy by dis up ing fluid mo ion, especially a he
smalle scales, which a e essen ial o accu a ely depic u bulen flows, making he use o conse a i e schemes mo e impo an . The
symme y-p ese ing me hod elimina es nume ical dissipa ion and conse es he physical p ope ies o he flow, while a he same
ime wa an ing uncondi ional s abili y, by mimicking he p ope ies o he con inuous ope a o s in hei disc e e coun e pa s [30].
This me hod has been ex ended o colloca ed g id a angemen s [31]and implemen ed in he open-sou ce code OpenFOAM [32].
The p essu e e o ha is in oduced by he applica ion o his me hod on colloca ed g ids emains as he la ges sou ce o nume ical
dissipa ion.
O he me hods o fil e p essu e oscilla ions wi hou he in oduc ion o nume ical dissipa ion ha e been a emp ed, such as
fil e ing p essu e modes ha lie on he ke nel o he wide-s encil Laplacian ope a o [33]. Howe e , h ough examina ion o he
connec ion be ween he mesh and his ke nel, i was shown ha his me hod only wo ks on Ca esian meshes [34], as he oscilla o y
pa o he ke nel anishes o mos uns uc u ed meshes and complex geome ies. A me hod ha can elimina e, o a leas balance,
bo h he checke boa d p oblem and nume ical dissipa ion is he e o e s ill sough a e . Mo eo e , he inadequacy o he ke nel
fil e ing me hod illus a ed ha a b oade and clea defini ion o checke boa ding is lacking, and mos published wo ks in li e a u e
use a quali a i e desc ip ion o he phenomenon.
The p esen wo k con ibu es o his opic by in oducing a clea defini ion o he checke boa d p oblem. This defini ion can be
used o quan i y he le el o checke boa ding on any geome y o mesh. Fu he mo e, a sol e algo i hm is de eloped which uses his
quan ifica ion me hod o dynamically balance he occu ence o p essu e oscilla ions and nume ical dissipa ion. The s uc u e o his
Jou nal o Compu a ional Physics 521 (2025) 113537
3
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Fig. 2. Geome ic pa ame e s o he ope a o s in Table 1, needed o es ablish he symme y-p ese ing scheme.
Table 1
Ma ices con aining he geome ic pa ame e s shown in Fig. 2, which o m he building blocks o all symme y-p ese ing ope a o s. 𝑚
and 𝑛deno e he numbe o aces and con ol olumes espec i ely.
Ope a o Dimensions Desc ip ion
𝑇𝑓𝑜,𝑇
𝑓𝑛 𝑛×𝑚 ace-owne and ace-neighbou connec i i y ma ices, con aining en y (𝑖, 𝑓) =1i cell 𝑖is connec ed o
ace 𝑓as an owne o as a neighbou espec i ely.
𝑁𝑠𝑚×3𝑚 ace-no mal ma ix wi h diagonal blocks
(𝑁𝑠𝑥,𝑁
𝑠𝑦,𝑁
𝑠𝑧)con aining he
(𝑥, 𝑦, 𝑧)-componen s o he
ace-no mals, 𝐧𝑓, espec i ely.
𝐴𝑠𝑚×𝑚diagonal ma ix con aining ace a eas, 𝐴𝑓.
𝛿𝑜
𝑛𝑠,𝛿𝑛
𝑛𝑠 𝑚×𝑚 ace-owne and ace-neighbou no mal dis ance diagonal ma ices con aining 𝛿𝑜
𝑛𝑓 and 𝛿𝑛
𝑛𝑓 espec i ely,
deno ing he absolu e alues o he 𝐧𝑓-p ojec ed ec o s om ace-cen oid o owne o neighbou
cen oid espec i ely.
Ω𝑐𝑛×𝑛cell- olume diagonal ma ix.
pape is as ollows: Sec ion 2gi es an o e iew o he equa ions and diffe en algo i hms ha a e used and discusses he occu ence
o checke boa ding and nume ical e o s. Sec ion 3discusses me hods o quan i y checke boa ding and in oduces a sol e ha
dynamically balances nume ical dissipa ion and checke boa d oscilla ions. Sec ion 4shows esul s o he new sol e compa ed o
exis ing me hods. Finally, sec ion 5discusses he conclusions o he wo k and he ou look o u u e wo k.
2. Nume ical amewo k
2.1. Symme y-p ese ing me hod
The semi-disc e ised o mula ion o equa ions (1)and (2) o a bi a y colloca ed g ids, using he ma ix- ec o no a ion o [31],
is gi en by:
Ω𝜕𝑡𝐮𝑐+𝐶(𝐮𝑠)𝐮𝑐=−𝐷𝐮𝑐−Ω𝐺𝑐𝐩𝑐,(3)
𝑀𝐮𝑠=𝟎𝑐.(4)
In h ee dimensions, he disc e e colloca ed kinema ic p essu e and eloci y fields a e gi en by 𝐩𝑐∈ℝ𝑛and 𝐮𝑐=(𝐮𝑇
𝑐,𝑥,𝐮𝑇
𝑐,𝑦,𝐮𝑇
𝑐,𝑧)𝑇
∈
ℝ3𝑛, whe eas 𝐮𝑠∈ℝ𝑚gi es he s agge ed eloci ies, in which 𝑛and 𝑚gi e he numbe o con ol olumes and aces espec i ely. A
small se o ma ices con aining geome ic in o ma ion abou he mesh, gi en in Table 1, is enough o de i e all he ma ix ope a o s
needed o om an algeb aic symme y-p ese ing amewo k, gi en in Table 2.
Midpoin in e pola ion in he con ec i e e m is necessa y o main ain he skew-symme y o he con inuous ope a o [30], whe eas
i was shown in [35–38] ha employing olume ic in e pola ion in he o he ope a o s leads o uncondi ional s abili y, e en on
highly-dis o ed meshes. In he ollowing ex he e o e, i supe sc ip 𝛾is d opped, olume ic in e pola ion is employed, e.g. 𝐿𝑐=
𝑀𝑐𝐺𝑐=−𝑀Γ𝑐𝑠Ω−1Γ𝑇
𝑐𝑠𝑀𝑇=−𝑀Γ𝑉
𝑐𝑠Ω−1Γ𝑉𝑇
𝑐𝑠 𝑀𝑇.
2.2. Time-s epping algo i hm
The empo al in eg a ion, which was le undisc e ised in equa ion (3), is aken ca e o by he ac ional s ep me hod [39], in which
a p edic o eloci y is calcula ed and co ec ed h ough sol ing a Poisson equa ion o p essu e. Usage o a wide-s encil Laplacian
o he p essu e Poisson equa ion in combina ion wi h he DIM o calcula e he s agge ed eloci ies leads o checke boa ding. By
applying a compac s encil and/o he WIM, his me hod can be adjus ed o deal wi h p essu e field oscilla ions. Table 3gi es an
o e iew o he ou possible ac ional s ep me hods ound by combining he compac - and wide-s encil me hods wi h he DIM and
Jou nal o Compu a ional Physics 521 (2025) 113537
4
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Table 2
Full se o ma ix ope a o s o o m an algeb aic symme y-p ese ing amewo k. 𝑚
and 𝑛deno e he numbe o aces and con ol olumes espec i ely. 𝛾indica es an un-
specified in e pola ion me hod, o which 𝐿, 𝑀and 𝑉a e gi en as op ions, indica ing
linea , midpoin and olume ic in e pola ion espec i ely.
Ope a o defini ion Dimensions Desc ip ion
Ω=𝐼3⊗Ω𝑐3𝑛×3𝑛colloca ed olumes
𝛿𝑛𝑠 =𝛿𝑜
𝑛𝑠 +𝛿𝑛
𝑛𝑠 𝑚×𝑚 ace-no mal dis ances
Ω𝑠=𝛿𝑛𝑠𝐴𝑠𝑚×𝑚s agge ed olumes
𝑆𝑠=𝐴𝑠𝑁𝑠𝑚×3𝑚su ace ec o s
𝑀=(𝑇𝑓𝑜 −𝑇𝑓𝑛)𝐴𝑠𝑛×𝑚di e gence
𝐺=𝛿−1
𝑛𝑠 (𝑇𝑓𝑛 −𝑇𝑓𝑜)𝑇=−Ω
−1
𝑠𝑀𝑇𝑚×𝑛g adien
𝐿=𝑀𝐺 =−𝑀Ω−1
𝑠𝑀𝑇𝑛×𝑛compac -s encil Laplacian
𝑊𝛾
𝑜=⎧
⎪
⎨
⎪
⎩
𝛿−1
𝑛𝑠 𝛿𝑛
𝑛𝑠,𝛾=𝐿
1
2𝐼𝑚,𝛾=𝑀
𝛿−1
𝑛𝑠 𝛿𝑜
𝑛𝑠,𝛾=𝑉
𝑚×𝑚in e pola ion weigh s
𝑊𝛾
𝑛=𝐼𝑚−𝑊𝛾
𝑜
Π𝛾
𝑐𝑠 =𝑊𝛾
𝑜𝑇𝑇
𝑓𝑜 +𝑊𝛾
𝑛𝑇𝑇
𝑓𝑛 𝑚×𝑛cell- o- ace in e pola o
Γ𝛾
𝑐𝑠 =𝑁𝑠(𝐼3⊗Π𝛾
𝑐𝑠)𝑚×3𝑛cell- o- ace do -in e pola o
Γ𝛾
𝑠𝑐 =Ω
−1Γ𝛾𝑇
𝑐𝑠 Ω𝑠3𝑛×𝑚 ace- o-cell in e pola o
𝑀𝛾
𝑐=𝑀Γ𝛾
𝑐𝑠 𝑛×3𝑛colloca ed di e gence
𝐺𝛾
𝑐=Γ
𝛾
𝑠𝑐 𝐺=−Ω
−1𝑀𝛾𝑇
𝑐3𝑛×𝑛colloca ed g adien
𝐿𝛾
𝑐=𝑀𝛾
𝑐𝐺𝛾
𝑐=−𝑀Γ𝛾
𝑐𝑠Ω−1Γ𝛾𝑇
𝑐𝑠 𝑀𝑇𝑛×𝑛wide-s encil Laplacian
𝐶𝑐(𝐮𝑠)=𝑀diag(𝐮𝑠)Π𝑀
𝑐𝑠 𝑛×𝑛con ec i e block
𝐶(𝐮𝑠)=𝐼3⊗𝐶
𝑐(𝐮𝑠)3𝑛×3𝑛con ec i e ope a o
𝐷𝑐=−𝜈𝐿 𝑛 ×𝑛diffusi e block
𝐷=𝐼3⊗𝐷
𝑐3𝑛×3𝑛diffusi e ope a o
Table 3
O e iew o he ou possible ac ional s ep me hods ound by combining he compac - and wide-
s encil me hods wi h he DIM and he WIM.
Wide s encil Compac s encil
DIM WIM DIM WIM
𝐮𝑝
𝑐=(𝐮𝑐,𝐮𝑠)𝐮𝑝∗
𝑐=𝐮𝑝
𝑐−𝐺𝑐
𝐩𝑝
𝑐
𝐿𝑐
𝐩𝑛+1
𝑐=𝑀𝑐𝐮𝑝
𝑐𝐿
𝐩′
𝑐=𝑀𝑐𝐮𝑝∗
𝑐

𝐩𝑛+1
𝑐=
𝐩𝑝
𝑐+
𝐩′
𝑐
𝐮𝑛+1
𝑐=𝐮𝑝
𝑐−𝐺𝑐
𝐩𝑛+1
𝑐𝐮𝑛+1
𝑠=Γ
𝑐𝑠𝐮𝑝
𝑐−𝐺
𝐩𝑛+1
𝑐
𝐮𝑛+1
𝑠=Γ
𝑐𝑠𝐮𝑛+1
𝑐𝐮𝑛+1
𝑠=Γ
𝑐𝑠𝐮𝑝
𝑐−𝐺
𝐩𝑛+1
𝑐𝐮𝑛+1
𝑐=Γ
𝑠𝑐 𝐮𝑛+1
𝑠𝐮𝑛+1
𝑐=𝐮𝑝
𝑐−𝐺𝑐
𝐩𝑛+1
𝑐
he WIM. The empo al disc e isa ion is no he main in e es in his wo k and is simply deno ed by a unc ion, (𝐮𝑐,𝐮𝑠), which
calcula es he eloci y p edic o , 𝐮𝑝
𝑐. As an example, o Fo wa d Eule ime-s epping, his e m is gi en by:
𝐮𝑝
𝑐=(𝐮𝑐,𝐮𝑠)
=𝐮𝑛
𝑐+Δ𝑡𝑅 (𝐮𝑛
𝑐,𝐮𝑛
𝑠)
=𝐮𝑛
𝑐−Δ𝑡Ω−1 (𝐶(𝐮𝑛
𝑠)+𝐷)𝐮𝑛
𝑐.
(5)
The alue o he p essu e p edic o , 
𝐩𝑝
𝑐is usually chosen o be 𝟎𝑐o 
𝐩𝑛
𝑐, co esponding o he classical Cho in p ojec ion me hod
[40]and he second-o de Van Kan p ojec ion me hod [9], espec i ely. No e ha his p essu e p edic o only has an effec when a
compac -s encil Laplacian is used. Fo he wide-s encil me hods, 𝐿𝑐
𝐩𝑝
𝑐could simply be added o bo h sides o he Poisson equa ion,
esul ing in 𝐿𝑐
𝐩𝑛+1
𝑐on he le hand side (LHS) and 𝑀Γ𝑐𝑠𝐮𝑝
𝑐on he igh hand side (RHS).
Upon close inspec ion, he compac -s encil DIM can immedia ely be disca ded o being less accu a e han he compac -s encil
WIM. This is ue because he only diffe ence be ween hese me hods is he back-and- o h in e pola ion o he p edic o eloci y.
This ope a ion can be iewed as he applica ion o a Laplacian fil e , as: 𝐮𝑛+1
𝑐=Γ
𝑠𝑐Γ𝑐𝑠𝐮𝑝
𝑐−𝐺𝑐
𝐩′
𝑐=𝑓(𝐮𝑝
𝑐)−𝐺𝑐
𝐩′
𝑐. This fil e can
be unde s ood mos easily by conside ing i s effec on a uni o m Ca esian mesh, whe e he classical [1,−2,1] coefficien s o he
Laplacian a e e ie ed:
𝑓(𝐮𝑝
𝑐)=𝐿𝑓𝐮𝑝
𝑐(uni o m Ca esian),(6)
𝐿𝑓=𝐼+1
4diag(𝐿𝑓,𝑥,𝐿
𝑓,𝑦,𝐿
𝑓,𝑧),(7)
[𝐿𝑓,𝑥𝐮𝑝
𝑐,𝑥]𝑖,𝑗,𝑘 =[𝐮𝑝
𝑐,𝑥]𝑖−1,𝑗,𝑘 −2[𝐮𝑝
𝑐,𝑥]𝑖+[𝐮𝑝
𝑐,𝑥]𝑖+1,𝑗,𝑘
,(8)
Jou nal o Compu a ional Physics 521 (2025) 113537
5
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
whe e subsc ip s 𝑖, 𝑗, 𝑘indica e cell numbe ing in 𝑥, 𝑦, 𝑧-di ec ions espec i ely. This unnecessa y ex a fil e ing is smoo hing he
eloci y field and he eby causing undesi ed nume ical dissipa ion. This me hod is he e o e no conside ed any u he .
The emaining me hods each ha e hei own p oblem, which makes choosing he igh me hod a ade-off be ween hese ac o s.
As discussed in he in oduc ion, he wide-s encil DIM is e y p one o he occu ence o checke boa ding. The emainde o his
sec ion discusses occu ence o checke boa ding o he WIMs, and he nume ical e o s ha hey in oduce when ying o diminish
he p oblem. As he e is no nume ical e o o he wide-s encil DIM o his ca ego y, his me hod is discussed no mo e he ea e . In
he wide-s encil WIM, a co ec ion is applied o he s agge ed eloci ies, which makes hem non-di e gence- ee as a esul . Al hough
ca ied ou in exac ly he same way, he calcula ion o he s agge ed eloci ies o he compac -s encil WIM is wi hou co ec ion,
and ollows di ec ly om he Poisson equa ion. In his me hod, howe e , he co ec ion is applied o he colloca ed eloci ies, which
in u n makes hem non-di e gence- ee. The di e gence o ei he he s agge ed o colloca ed eloci ies in oduces a nume ical e o
o he simula ion, which is u he discussed in sec ion 2.4.
2.3. Occu ence o checke boa ding
Al hough he WIM, o bo h he compac - and wide-s encil me hods, in oduces a coupling o he p essu e field be ween neigh-
bou ing nodes a he cos o a nume ical e o , hese me hods can s ill show oscilla ing p essu e fields, mos commonly caused by
he usage o a small ime-s ep in uns eady simula ions [25], o by he inclusion o a p essu e p edic o as 
𝐩𝑝
𝑐=
𝐩𝑛
𝑐in case o he
compac -s encil WIM [32]. To illus a e his, no e ha i Δ𝑡 →0+, hen 𝐮𝑝
𝑐→𝐮𝑛
𝑐in equa ion (5), since he effec o he con ec i e and
diffusi e e ms a e p opo ional o Δ𝑡. This emo es he coupling es ablished in he wide-s encil WIM, since i is dependen on 𝐮𝑠
and he con ec i e e m. To show he decoupling o he compac -s encil WIM, fi s ega d he case o 
𝐩𝑝
𝑐=𝟎𝑐, wi h Fo wa d Eule
empo al disc e isa ion as an example, which leads o:
𝐮𝑝∗
𝑐=𝐮𝑛
𝑐,(9)
𝐮𝑛
𝑐=𝐮𝑛−1
𝑐−𝐺𝑐
𝐩𝑛
𝑐=𝐮0
𝑐−𝐺𝑐
𝑛
∑
𝑖

𝐩𝑖
𝑐,(10)
𝐿
𝐩𝑛+1
𝑐=𝑀𝑐𝐮𝑝∗
𝑐=𝑀𝑐𝐮0
𝑐−𝐿𝑐
𝑛
∑
𝑖

𝐩𝑖
𝑐,(11)
𝐿ℙ𝑛+1
𝑐=𝑀𝑐𝐮0
𝑐+(𝐿−𝐿𝑐)ℙ𝑛
𝑐,(12)
whe e ℙ𝑛
𝑐=∑𝑛
𝑖
𝐩𝑖
𝑐and 𝐿ℙ𝑛
𝑐is added on bo h sides o equa ion (11) o each equa ion (12). Equa ion (12)gi es a s a iona y i e a i e
me hod o sol e he wide-s encil Poisson equa ion, o which he solu ion is a decoupled p essu e field. The coupling ha he compac -
s encil Laplacian ga e is he e o e los i he ime-s ep becomes oo small. Using he Van Kan me hod, whe e 
𝐩𝑝
𝑐=
𝐩𝑛
𝑐, leads o mo e
se e e oscilla ions in he p essu e field. In his case he wide-s encil Laplacian ope a es on a la ge pa o he p essu e field:
𝐿
𝐩′
𝑐=𝑀𝑐𝐮𝑝∗
𝑐=𝑀𝑐𝐮𝑝
𝑐−𝐿𝑐
𝐩𝑝
𝑐,(13)
which weakens he coupling ha he compac -s encil Laplacian p o ided.
Since he p essu e oscilla ions a e “in isible” o he wide-s encil g adien , he colloca ed eloci ies emain unaffec ed and he
algo i hm ad ances as i hey we e no he e. The oscilla ions a e e ained in his case and could g ow un il hey each a poin in
which hey cause nume ical issues and uns able solu ions. Al hough hei e en ion is e iden , hei ac ual o igins and me hod o
g ow h a e no en i ely unde s ood. In [31]i was a gued ha he con ec i e e m causes spu ious modes in he eloci y field ha
lead o checke boa ding, bu an analysis o his p ocess was no pe o med. A ela ion o he me hod wi h which he Poisson equa ion
is sol ed, including p econdi ione s, was discussed in [41], bu he exac mechanisms emain unclea .
2.4. Nume ical e o s o he WIM
To analyse he nume ical e o s o he WIM, he global disc e e kine ic ene gy is conside ed, gi en by 𝐸𝐾=1
2𝐮𝑇
𝑐Ω𝐮𝑐. The empo al
e olu ion o his e m is de i ed using he p oduc ule and equa ion (3):
𝑑
𝑑𝑡𝐸𝐾=−1
2⎛⎜⎜⎝
𝐮𝑇
𝑐(𝐶(𝐮𝑠)+𝐶𝑇(𝐮𝑠))𝐮𝑐
+𝐮𝑇
𝑐(𝐷+𝐷𝑇)𝐮𝑐
+𝐮𝑇
𝑐Ω𝐺𝑐𝐩𝑐+𝐩𝑇
𝑐𝐺𝑇
𝑐Ω𝑇𝐮𝑐⎞⎟⎟⎠
.(14)
I 𝐶(𝐮𝑠)is symme y-p ese ing, i.e. 𝐶(𝐮𝑠) =−𝐶𝑇(𝐮𝑠), mimicking he con inuous ope a o , hen he con ibu ion o he con ec i e
e m equals ze o. Simila ly, he con ibu ion o he p essu e e m equals ze o i Ω𝐺𝑐=−𝑀𝑇
𝑐and 𝑀𝑐𝐮𝑐=𝟎𝑐. The global kine ic
ene gy dissipa ion is he e o e educed o he iscous e m, equaling:
𝑑
𝑑𝑡𝐸𝐾=−1
2𝐮𝑇
𝑐(𝐷+𝐷𝑇)𝐮𝑐≤0,(15)
which is s ic ly dissipa i e i 𝐷is cons uc ed using he symme y-p ese ing disc e isa ion, e aining he posi i e-defini eness o he
con inuous ope a o .

Jou nal o Compu a ional Physics 521 (2025) 113537
6
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
2.4.1. Wide-s encil WIM -con ec i e e o
Fo he wide-s encil WIM, 𝑀𝑐𝐮𝑛+1
𝑐=𝟎𝑐and he e o e he p essu e e o equals ze o. Con e sely, he con ec i e e o is non-ze o.
This is caused by he co ec ion o he s agge ed eloci ies and he esul ing ac ha 𝑀𝐮𝑛+1
𝑠≠0, since i is gi en by:
𝑀𝐮𝑛+1
𝑠=𝑀(Γ𝑐𝑠𝐮𝑝
𝑐−𝐺
𝐩𝑛+1
𝑐)
=

: 𝟎𝑐
𝑀Γ𝑐𝑠 (𝐮𝑝
𝑐−𝐺𝑐
𝐩)+𝑀(Γ𝑐𝑠Γ𝑠𝑐 −𝐼)𝐺
𝐩𝑛+1
𝑐
=(𝐿𝑐−𝐿)
𝐩𝑛+1
𝑐.
(16)
To see he effec on he e olu ion o kine ic ene gy, ecall ha 𝐶𝑐(𝐮𝑠) =𝑀diag(𝐮𝑠)Π𝑚
𝑐𝑠. The appea ance o 𝑀in his defini ion
and he incomp essibili y cons ain o equa ion (4), is ela ed o he simplified ad ec i e o m in which equa ion (1)is w i en. The
off-diagonals o 𝐶𝑐(𝐮𝑠)con ain he fluxes be ween cells, and a e skew-symme ic by cons uc ion. In con as he diagonal en ies a e
gi en by:
𝐶𝑖,𝑖(𝐮𝑠)= 1
2[𝑀𝐮𝑠]𝑖,(17)
which a e non-ze o in his case. Equa ion (17)is consis en wi h he si ua ion o comp essible flows, see o example [42]. Howe e ,
in he incomp essible amewo k, his e m is no balanced by he mass conse a ion equa ion. I 𝐶𝑐(𝐮𝑠)is spli in o i s non-ze o
diagonal and skew-symme ic off-diagonals as 𝐶𝑐(𝐮𝑠) =1
2diag(𝑀𝐮𝑠)+𝐶𝑂𝐷
𝑐(𝐮𝑠), hen he nume ical e o o he e olu ion o he
global kine ic ene gy a ime-s ep 𝑛 +1due o he con ec i e e o is gi en by:
𝑑
𝑑𝑡𝐾=− 1
2𝐮𝑛+1𝑇
𝑐⎛⎜⎜⎜⎝
: 𝟎3𝑛×3𝑛
𝐶𝑂𝐷(𝐮𝑛+1
𝑠)+𝐶𝑂𝐷𝑇(𝐮𝑛+1
𝑠)⎞⎟⎟⎟⎠
𝐮𝑛+1
𝑐
−1
2𝐮𝑛+1𝑇
𝑐[𝐼3⊗diag(𝑀𝐮𝑛+1
𝑠)]𝐮𝑛+1
𝑐.
(18)
The con ec i e e o in he wide-s encil WIM is he e o e caused by he di e gence o he s agge ed eloci ies and p opo ional o
Δ𝑡
(𝐿𝑐−𝐿)𝐩𝑛+1
𝑐.
2.4.2. Compac -s encil WIM - p essu e e o
On he o he hand, o he compac -s encil WIM, 𝑀𝐮𝑛+1
𝑠=𝟎𝑐and he e o e he con ec i e e o equals ze o. Howe e , in his
case he co ec ion o he colloca ed eloci ies leads o 𝑀𝑐𝐮𝑛+1
𝑐≠𝟎𝑐, which is gi en by:
𝑀𝑐𝐮𝑛+1
𝑐=𝑀Γ𝑐𝑠 (𝐮𝑝
𝑐−𝐺𝑐
𝐩′
𝑐)
=

: 𝟎𝑐
𝑀(Γ𝑐𝑠𝐮𝑝
𝑐)−𝐺
𝐩′
𝑐+𝑀(𝐼−Γ
𝑐𝑠Γ𝑠𝑐 )𝐺
𝐩′
𝑐
=(𝐿−𝐿𝑐)
𝐩′
𝑐.
(19)
Then he ensuing nume ical e o in he e olu ion o he global kine ic ene gy a ime-s ep 𝑛 +1due o he p essu e e o is gi en by:
𝑑
𝑑𝑡𝐾=−𝐮𝑛+1𝑇
𝑐Ω𝐺𝑐𝐩′
𝑐
=𝐩′𝑇
𝑐𝑀𝑐𝐮𝑛+1
𝑐
=𝐩′𝑇
𝑐(𝐿−𝐿𝑐)
𝐩′
𝑐.
(20)
In he compac -s encil WIM, he p essu e e o is caused by he di e gence o he colloca ed eloci ies and is, simila ly o he con ec i e
e o , p opo ional o Δ𝑡
(𝐿−𝐿𝑐)𝐩′
𝑐.
In p ac ice, mos nume ical codes e e o he compac -s encil WIM. The main easons o his a e h ee old: (i) The inclu-
sion o a p essu e p edic o as 
𝐩𝑝
𝑐=
𝐩𝑛
𝑐 educes he o de o he p essu e e o om 
(Δ𝑡2) o 
(Δ𝑡4)since 
𝐩𝑛+1
𝑐∼
(Δ𝑡)and
(
𝐩𝑛+1
𝑐−
𝐩𝑛
𝑐)∼
(Δ𝑡2), making he non-physical con ibu ion qui e small. This educ ion is no possible wi h he wide-s encil Lapla-
cian. (ii) The compac -s encil Laplacian educes he compu a ional complexi y and cos when using uns uc u ed g ids. And finally,
(iii) he con ec i e e o e m, as shown in equa ion (18), is no s ic ly dissipa i e. Since he global di e gence o he s agge ed
eloci ies will always be equal o ze o, he e will always be cells wi h nega i e di e gence o compensa e cells ha ha e posi i e
di e gence. This means ha kine ic ene gy can be added in some pa s o he solu ion, leading o ins abili ies. Con e sely, he e-
sul ing e m o equa ion (20)will be dependen on he e m
(𝐼−Γ
𝑐𝑠Γ𝑠𝑐 ), which in u n depends on he choice in in e pola o s and
meshing. These ac o s a e much easie o con ol, and mo eo e , i was shown in [35–38] ha choosing olume ic in e pola ion o
his ope a o causes his e m o be s ic ly dissipa i e on ci cumcen e meshes. In mos p ac ical meshes, s abili y can be wa an ed
in his way.
Jou nal o Compu a ional Physics 521 (2025) 113537
7
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
3. Defining he checke boa d p oblem
Pas wo ks ha ex ensi ely discuss he opic use diffe en names o he checke boa d p oblem, such as checke boa ding [43,17,
44], spu ious p essu e modes [22], odd-e en decoupling [45], o zigzagness [46], bu none o hem gi e a quan ifiable defini ion o
he p oblem. In hose wo ks a cons an and uni o m solu ion is applied ha does no need quan ifica ion o he p oblem, whe eas in
his wo k a solu ion is sough a e ha is applied p opo ionally o he se e i y o he p oblem. Fo his eason a clea quan ifica ion
me hod has o be in oduced fi s .
3.1. Using he ke nel o he disc e e Laplacian o define checke boa ding
One way o define checke boa d modes is by pe o ming an eigen ec o decomposi ion o 𝐿𝛾
𝑐and defining he checke boa d
modes as 𝐩−
𝑐∈𝐾𝑒𝑟(𝐿𝛾
𝑐). Since 𝐿𝛾
𝑐=𝑀𝛾
𝑐𝐺𝛾
𝑐=𝑀Γ𝛾
𝑐𝑠Γ𝛾
𝑠𝑐𝐺, he 𝐩−
𝑐modes will depend on he choice o in e pola o and he mesh. The
cons an mode ec o will always lie on he ke nel, in addi ion o any spu ious mode ec o s. This me hod was o ins ance used by
[33], whe e he me hod was adequa e since midpoin in e pola ion and Ca esian meshes we e used. By looking a he defini ions o
𝑀𝛾
𝑐and 𝐺𝛾
𝑐i becomes clea why his null-space exis s o hese specific condi ions, since:
[Γ𝛾
𝑠𝑐𝐺𝜙
𝜙
𝜙𝑐]𝑖=1
Ω𝑖∑
𝑓∈𝐹𝑓(𝑖)
𝑤𝛾
𝑓𝑖Ω𝑓
𝜙
𝜙
𝜙𝑗−𝜙
𝜙
𝜙𝑖
𝛿𝑛𝑓
𝐧𝑓(𝑖)
=1
Ω𝑖∑
𝑓∈𝐹𝑓(𝑖)
𝑤𝛾
𝑓𝑖𝜙
𝜙
𝜙𝑗𝐬𝑓(𝑖)−1
Ω𝑖∑
𝑓∈𝐹𝑓(𝑖)
𝑤𝛾
𝑓𝑖𝜙
𝜙
𝜙𝑖𝐬𝑓(𝑖),
(21)
gi es he wide-s encil g adien a cell 𝑖. Whe e 𝜙
𝜙
𝜙𝑐∈ℝ𝑛×1 and 𝜙
𝜙
𝜙𝑖=[𝜙
𝜙
𝜙𝑐]𝑖. Cell 𝑗is neighbou ing cell 𝑖 h ough ace 𝑓. 𝐹𝑓(𝑖)deno es
he se o aces ha define cell 𝑖. Ω𝑖=[Ω𝑐]𝑖,𝑖, Ω𝑓=[Ω𝑠]𝑓,𝑓 and 𝛿𝑛𝑓 =[𝛿𝑛𝑠]𝑓. 𝑤𝑖𝑓 equals
[𝑊𝛾
𝑜]𝑓o
[𝑊𝛾
𝑛]𝑓depending on whe he 𝑖is
he owne o he neighbou o ace 𝑓, espec i ely. Finally, 𝐧𝑓(𝑖)and 𝐬𝑓(𝑖) espec i ely gi e he ou wa d-poin ing ace-no mal ec o
and ou wa d-poin ing su ace ec o a ace 𝑓wi h espec o cell 𝑖, no e ha 𝐬𝑓(𝑖)=𝐴𝑓𝐧𝑓(𝑖). The second e m o he RHS o equa ion
(21) anishes in case 𝑤𝛾
𝑓𝑖 is cons an o e all aces o 𝑖, since
∑𝑓∈𝐹𝑓(𝑖)𝐬𝑓(𝑖)=𝟎. This is he case o midpoin in e pola ion and o
uni o m Ca esian meshes. Simila ly o he di e gence ope a o :
[𝑀Γ𝛾
𝑐𝑠𝜓
𝜓
𝜓𝑐]𝑖=∑
𝑓∈𝐹𝑓(𝑖)(𝑤𝛾
𝑖𝑓 𝜓
𝜓
𝜓𝑖+𝑤𝛾
𝑗𝑓𝜓
𝜓
𝜓𝑗)⋅𝐬𝑓(𝑖)
=∑
𝑓∈𝐹𝑓(𝑖)
𝑤𝛾
𝑖𝑓 𝜓
𝜓
𝜓𝑖⋅𝐬𝑓(𝑖)+∑
𝑓∈𝐹𝑓(𝑖)
𝑤𝛾
𝑗𝑓𝜓
𝜓
𝜓𝑗⋅𝐬𝑓(𝑖),
(22)
in which, again, he fi s e m o he RHS anishes o midpoin in e pola ion and o uni o m Ca esian meshes. When combining
equa ions (21)and (22)i becomes e iden ha 𝐿𝛾
𝑐does no connec cell 𝑖 o i s di ec neighbou s 𝑗, bu only o i s second-neighbou s
𝑘, when midpoin in e pola ion is used:
[𝐿𝑀
𝑐]𝑖,𝑗 =0,[𝐿𝑀
𝑐]𝑖,𝑘 =𝐴𝑓𝐴𝑔
4Ω𝑗
𝐬𝑓(𝑖)⋅𝐬𝑔(𝑗),(23)
in which ace 𝑓lies be ween cells 𝑖and 𝑗, and ace 𝑔lies be ween cells 𝑗and 𝑘. I his odd-e en pa i y is sus ained h oughou he
en i e mesh, wo disconnec ed g oups o cells will exis . In addi ion o he cons an ke nel ec o , his gi es ise o a spu ious ke nel
ec o :
[𝐩−
𝑐]𝑖=(−1)
𝑝(𝑖),(24)
whe e 𝑝(𝑖)deno es he pa i y o cell 𝑖, equal o 0 o 1. In he special case o Ca esian meshes, he do p oduc in equa ion (23)will be
ze o o any diagonal second-neighbou pai ing, gi ing ise o 2𝑁𝑑𝑖𝑚 ke nel ec o s, wi h numbe o dimensions 𝑁𝑑𝑖𝑚. Fo example,
in h ee-dimensional Ca esian meshes he esul ing se o eigh ec o s is gi en by [33]:
[𝐩−
𝑐(𝐼𝐽𝐾)]𝑖,𝑗,𝑘 =(−1)
𝑖𝐼+𝑗𝐽+𝑘𝐾 ,(25)
in which 𝐼, 𝐽, 𝐾∈{0, 1} and indices 𝑖, 𝑗, 𝑘indica e he cell numbe ing in each o he Ca esian di ec ions. In his no a ion, 𝐩−
𝑐(000) gi es
he cons an ec o .
[𝐿𝛾
𝑐]𝑖,𝑗 is gene ally only ze o o midpoin in e pola ion. Whe eas, in gene al,
[𝐿𝛾
𝑐]𝑖,𝑗 ≠0 o linea o olume ic
in e pola ion, excep on uni o m meshes whe e all in e pola o s equal midpoin in e pola ion. Howe e , a se o ke nel ec o s was
de i ed by [34] o any Ca esian mesh wi h midpoin , linea o olume ic in e pola ion. To do so, he ac ha 𝐾𝑒𝑟(𝐺𝛾
𝑐) ∈𝐾𝑒𝑟(𝐿𝛾
𝑐)
was used, in addi ion o ew i ing 𝐺𝛾
𝑐 o 𝐺𝐺Π𝛾
𝑐𝑠. 𝐺𝐺deno es a Gauss-g adien and he o e ba deno es swapping he in e pola ion
weigh s be ween owne s and neighbou s o he aces, such ha Π𝑀
𝑐𝑠 = Π𝑀
𝑐𝑠 , Π𝑉
𝑐𝑠 = Π𝐿
𝑐𝑠 and Π𝐿
𝑐𝑠 = Π𝑉
𝑐𝑠. Fo de ails on his equali y, he
eade is e e ed o Appendix A. By aking he in e pola ion be o e he Gauss-g adien , cell-cen e ed alues can be chosen such ha
se s o opposing in e pola ed ace alues a e always equal, o which he cell-cen e ed Gaus-g adien equals ze o. The esul ing se o
eigh ec o s o h ee-dimensional Ca esian meshes is hen gi en by:
Jou nal o Compu a ional Physics 521 (2025) 113537
8
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Fig. 3. P essu e oscilla ion in a one-dimensional pe iodic domain, o hogonal o 𝐾𝑒𝑟(𝐿𝑐),bu de ec ableby𝐶𝑐𝑏.
[𝐩−
𝑐(𝐼𝐽𝐾,𝛾)]𝑖,𝑗,𝑘 =(−1)
𝑖𝐼+𝑗𝐽+𝑘𝐾 ([Δ𝑥]𝐼
𝑖[Δ𝑦]𝐽
𝑗[Δ𝑧]𝐾
𝑘)𝛼
,(26)
in which 𝛼={−1, 0, 1} o linea , midpoin and olume ic in e pola ions espec i ely. This se is no necessa ily mu ually o hogonal,
especially o linea and olume ic in e pola ions, howe e , in all cases he se is linea ly independen and he e o e spans he null-
space o 𝐿𝛾
𝑐. The de i a ion o his se o ec o s wi h an example is gi en in Appendix B.
Despi e his ex ension, calcula ing he ke nel o 𝐿𝛾
𝑐 o non-Ca esian meshes in ol es pe o ming a singula alue decomposi ion,
o which he compu a ional cos g ows exponen ially wi h he numbe o g id poin s as (𝑁3
𝑔𝑟𝑖𝑑), which quickly becomes unaffo dable
[47]. Mo eo e , o mos meshes he ank o he ke nel educes o one i he pa i y and he o hogonal aces disappea , which is nea ly
always he case o any uns uc u ed mesh, lea ing only he cons an mode ec o . The e o e using he ke nel o 𝐿𝛾
𝑐 o define and
quan i y he checke boa d p oblem is o en insufficien , leading o he sea ch o a mo e gene ally applicable defini ion.
3.2. Applying a gene al defini ion o checke boa ding
In he p e ious sec ion a e y es ic i e defini ion o checke boa ding was used, o which some easy examples we e gi en in
which case he defini ion would no be ui ul. A mo e use ul defini ion can be ound i he p oblem is ega ded a con ol olume
le el. The essen ial p oblem o he decoupled con ol olume is ha he p essu e g adien o e a gi en ace, [𝐺𝐩𝑐]𝑓, migh gi e a
significan non-ze o alue, while he alue a he adjacen cell [𝐺𝑐𝐩𝑐]𝑖can be (close o) ze o. This p oblem migh occu only in a ew
cells and he e o e lie mos ly ou side o he ke nel o 𝐿𝑐, ha is, i he ke nel e en con ains spu ious ec o s.
I u ns ou ha he a io be ween he 𝐿2no ms o ec o s 𝐺𝑐𝐩𝑐and 𝐺𝐩𝑐gi es a good global indica ion o his decoupling. The
exp ession o he 𝐿2no ms o cell-cen e ed and ace-cen e ed fields is gi en by ‖𝑎𝑐‖=𝑎𝑇
𝑐Ω𝑐𝑎𝑐and ‖𝑎𝑠‖=𝑎𝑇
𝑠Ω𝑠𝑎𝑠, espec i ely.
Fu he mo e, in he lowe limi he field lies ully inside he ke nel, esul ing in a a io o ze o, whe eas a pe ec ly smoo h field
gi es he uppe limi , esul ing in a a io o one. Since a highe coefficien should indica e mo e p e alen checke boa ding, he
checke boa d coefficien is finally ound by sub ac ing his a io om one, esul ing in:
𝐶𝑐𝑏 (𝐩𝑐)=1−‖𝐺𝑐𝐩𝑐‖
‖𝐺𝐩𝑐‖=1−𝐩𝑇
𝑐𝐺𝑇
𝑐Ω𝐺𝑐𝐩𝑐
𝐩𝑇
𝑐𝐺𝑇Ω𝑠𝐺𝐩𝑐
=𝐩𝑇
𝑐(𝐿−𝐿𝑐)𝐩𝑐
𝐩𝑇
𝑐𝐿𝐩𝑐
.(27)
𝐶𝑐𝑏 (𝐩𝑐)is defined as ze o o he cons an p essu e field in which case ‖𝐺𝐩𝑐‖ =0. As he 𝐿2no ms o bo h p essu e g adien fields a e
non-nega i e, he coefficien ’s uppe bound is 1. The lowe bound is 0 as long as 𝐩𝑇
𝑐(𝐿−𝐿𝑐)𝐩𝑐is non-posi i e, i.e. 𝐿 −𝐿𝑐is nega i e
semidefini e. This holds ue o olume ic in e pola ion and ci cumcen e meshes [38]. Fo i ually all simula ions in p ac ice, he
eigen alues o 𝐿 −𝐿𝑐we e non-posi i e o a leas e y close o ze o and non-p oblema ic, as long as olume ic in e pola ion was
used, leading o a alue o 𝐶𝑐𝑏 in he in e al [0,1]. The checke boa d coefficien can be calcula ed o any colloca ed scala field.
As p essu e field oscilla ions a e he ocus o his wo k, 𝐶𝑐𝑏 (𝐩𝑐)will hence o h simply be deno ed by 𝐶𝑐𝑏. Aside om being a
non-dimensional and no malised coefficien , ew i ing he coefficien as done in equa ion (27) e eals some o he p ope ies which
suppo his defini ion.
Fi s ly, i e eals ha he magni ude o he p essu e fields does no ma e , since he cons an p essu e ec o lies inside he ke nel
o 𝐿and 𝐿𝑐and will he e o e be fil e ed in he ope a ion o equa ion (27). This should always be he case o any defini ion ha is
employed, since he alue o p essu e is no wha ma e s in incomp essible flows, a he i s g adien .
Secondly, he ma ix
(𝐿−𝐿𝑐)=𝑀(𝐼−Γ
𝑐𝑠Γ𝑠𝑐 )𝐺is inna ely linked o he undamen al p oblem o he colloca ed g id a ange-
men , which is ha he g adien o p essu e and he eloci ies a e no defined a he same loca ion. The ensuing in e pola o s will
always apply some smoo hing o a field such ha only 𝐼≈Γ
𝑐𝑠Γ𝑠𝑐 .
Thi dly, his defini ion has a physical meaning, since he nume ical e o s in he con ec i e and p essu e e ms, as discussed in
sec ions 2.4.1 and 2.4.2 espec i ely, a e p opo ional o he e m ±
(𝐿−𝐿𝑐)𝐩𝑐. Fo he wide-s encil WIM,
(𝐿𝑐−𝐿)𝐩 =𝑀𝐮𝑠, which
is non-ze o and leads o he e o exp essed in equa ion (18). Simila ly, o he compac -s encil WIM,
(𝐿𝑐−𝐿)𝐩 =𝑀𝑐𝐮𝑐, which is
non-ze o and leads o he p essu e e o as seen in equa ion (20). One no able diffe ence is ha 𝐶𝑐𝑏 is independen o he ime-s ep,
which is also a necessa y p ope y, since i should be possible o calcula e he coefficien wi hou any knowledge o he empo al
disc e isa ion.
Finally, 𝐶𝑐𝑏 is able o de ec local oscilla ions which lie ou side o he ke nel o 𝐿𝑐, which i should be able o do, as men ioned
be o e. A simple example o illus a e his poin can be p o ided wi h an oscilla ion on a one-dimensional pe iodic domain, as seen
in Fig. 3. The ollowing calcula ions apply o his configu a ion:
Jou nal o Compu a ional Physics 521 (2025) 113537
9
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
𝐩𝑐=[0010
10
],(28)
𝐩−
𝑐(1,𝛾)=[111111],(29)
𝐺𝐩𝑐=[01
111 0
],(30)
𝐺𝑐𝐩𝑐=[01
2010 1
2],(31)
𝐩𝑇
𝑐𝐩−
𝑐(1,𝑣)=0,(32)
𝐶𝑐𝑏 =1 − ‖𝐺𝑐𝐩𝑐‖
‖𝐺𝐩𝑐‖=5
8,(33)
whe e column ec o s a e ep esen ed ho izon ally o eadabili y. In ui i ely, his wiggle should no ha e a ze o alue when quan-
i ying he checke boa d p oblem, he e o e 𝐶𝑐𝑏 gi es a mo e desi able ou come han he ke nel ec o me hod.
3.3. One possible applica ion o he checke boa d coefficien
Since 𝐶𝑐𝑏 gi es a global, non-dimensional, no malised coefficien o checke boa ding, i can di ec ly be applied in he sol e
algo i hm o diminish he occu ence o he p oblem. I could, o example, be used o shi be ween he possible ac ional s ep me hod
algo i hms displayed in Table 3, whe e choices a e made conce ning: (i) he wid h o he Laplacian s encil, (ii) he in e pola ion
me hod o he eloci y co ec ion and (iii) he calcula ion o he p essu e p edic o . Since mos colloca ed fini e olume ac ional s ep
me hod codes apply he compac -s encil WIM, and since he inclusion o he p essu e p edic o is a known cause o checke boa ding,
one way o use 𝐶𝑐𝑏 is o deli e a nega i e eedback h ough his p edic o alue. To his end, he ollowing exp ession o 𝐩𝑝
𝑐is used:
𝐩𝑝
𝑐=𝜃𝑝𝐩𝑛
𝑐,(34)
𝜃𝑝=1−𝐶𝑐𝑏 =𝐩𝑛𝑇
𝑐𝐿𝑐𝐩𝑛
𝑐
𝐩𝑛𝑇
𝑐𝐿𝐩𝑛
𝑐
.(35)
By doing so, he sol e will con e ge o 𝜃𝑝=1in absence o checke boa ding in he p essu e field, benefi ing om he lowe nume ical
dissipa ion his offe s. Whe eas i he case o mesh is p one o checke boa ding, he algo i hm will end o 𝜃𝑝=0, in which case he
inc ease in nume ical dissipa ion can damp he oscilla ions. This eedback on o he p essu e p edic o offe s a dynamical balance
be ween wo p oblems, so ha he use o he code does no ha e o decide which p oblem will be mo e p e alen , c ea ing a unified
sol e ha should pe o m well in bo h si ua ions. Fo his ini ial a emp a es ablishing such a nega i e eedback, a simple linea
ela ion is used o explo e i s effec s, o he ela ions migh also wo k and could gi e be e esul s.
4. Resul s
In his sec ion, he new sol e wi h he p edic o p essu e gi en by equa ions (34)and (35)is es ed and compa ed o he Cho in
and Van Kan me hods, which use 𝜃𝑝=0and 𝜃𝑝=1, espec i ely. The alues o 𝜃𝑝 o he esul ing h ee me hods, labelled 𝜃0, 𝜃1
and 𝜃𝑝, a e gi en in Table 4. These h ee me hods a e es ed and compa ed in his sec ion. To do so, fi s ly, nume ical dissipa ion
is measu ed using a wo-dimensional Taylo -G een o ex. Secondly, a empo al and spa ial con e gence s udy is pe o med using a
wo-dimensional lid-d i en ca i y. And finally, a u bulen channel flow case is used o measu e o e all accu acy o he sol e and
o s udy he beha iou o he checke boa d coefficien in uns eady flows. All cases we e un wi h OpenFOAM, using he symme y-
p ese ing spa ial disc e isa ion and Runge-Ku a 3 empo al in eg a ion, which a e bo h implemen ed in he RKSymFoam sol e
de eloped o [32,48].
Table 4
Se ings o he es ed sol e s.
𝜃0𝜃1𝜃𝑑𝑦
𝜃𝑝011−𝐶𝑐𝑏
In each case he newly implemen ed sol e s a e compa ed o a eadily implemen ed sol e a ailable in OpenFOAM, which is
ei he icoFoam o lamina flows o pisoFoam o u bulen flows [28]. Apa om some choices in spa ial disc e isa ion, he bigges
diffe ence o he OpenFOAM sol e s is how hey emo e ime-s ep dependency o he compac -s encil WIM, which is simila o and
de i ed om he me hods used by [16,21]. Fo explici ime-s epping i is gi en by:
𝐮𝑝
𝑠=Γ
𝛾
𝑐𝑠𝐮𝑝
𝑐−𝐶𝑢𝑠𝐮𝑠,𝑐𝑜𝑟𝑟,(36)
𝐮𝑠,𝑐𝑜𝑟𝑟 =𝐮𝑛
𝑠−Γ
𝛾
𝑐𝑠𝐮𝑛
𝑐,(37)
[𝐶𝑢𝑠]𝑓,𝑓 =1−𝑚𝑖𝑛 (||[𝐮𝑠,𝑐𝑜𝑟𝑟]𝑓||
||[𝐮𝑛
𝑠]𝑓||+𝜖,1),(38)
whe e 𝜖is a e y small numbe o a oid di ision by ze o. This e m s abilises he esul s a he cos o in oducing a sizable amoun
o nume ical dissipa ion, which was shown by [49].
Jou nal o Compu a ional Physics 521 (2025) 113537
16
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
𝐷𝜖
𝑘=−𝜈⟨𝜕𝑗𝑢𝑖𝜕𝑗𝑢𝑖⟩+𝜈⟨𝜕𝑗𝑢𝑖⟩⟨𝜕𝑗𝑢𝑖⟩,
𝐶𝑇
𝑘=𝐶𝑘−𝐶𝑃
𝑘+⟨𝑢𝑗⟩𝜕𝑗𝑘,
𝐷𝑣
𝑘=𝐷𝑘−𝐷𝜖
𝑘,(48)
such ha he only choice in disc e isa ion ha emains o be made, is how o ake a cell-cen e ed g adien o he eloci y, 𝜕𝑗𝑢𝑖, as
⟨𝑢𝑗⟩𝜕𝑗𝑘 =0 o he channel configu a ion. In his wo k, he cell-cen e ed g adien is aken o each componen o 𝐮𝑐sepa a ely using
𝐺𝑐, c ea ing a enso field wi h 9 componen s a he cell-cen e s. The de i a ion o he equa ions and he implemen a ion o he
me hod in OpenFOAM can be ound in he Gi Hub eposi o y o [55].
The implemen ed me hods ely on a Runge-Ku a based sol e , whe eas pisoFoam elies on he empo al in eg a ion me hods
implemen ed in OpenFOAM. To make a ai compa ison be ween hese me hods he Backwa d Eule empo al in eg a ion was chosen,
since i is one o ew schemes a ailable o bo h sol e s. The e o e, he esul s on he coa se mesh, shown in Figs. 12a h ough 13e,
con ain a ai ly la ge e o , which is mainly caused by he Backwa d Eule ime scheme. The emphasis o hese figu es, howe e , is o
show he diffe ence be ween he implemen ed me hods and pisoFoam, no he o e all quali y o he solu ion. To show he influence
o he empo al e o , he esul s o he 𝜃0me hod wi h a Runge-Ku a 3 scheme a e shown in Figs. 12a, 12band 12c. F om hese
figu es i can also be seen ha e en on a coa se mesh, he implemen ed sol e s a e able o a ain ai ly decen esul s. F om Fig. 12a
i can be seen ha he e he mean s eam-wise eloci y o all sol e s is gene ally much oo high. pisoFoam shows he highes mean
s eam-wise eloci y, indica ing ha i is less able o de elop u bulence on his mesh, which may be due o i s dissipa i e na u e.
This finding is confi med by Fig. 12b, in which he iscous shea s ess, 𝜌𝜈 𝑑⟨𝑢𝑥⟩
𝑑𝑦 , and he Reynolds s ess, 𝜌⟨𝑢′
𝑥𝑢′
𝑦⟩a e plo ed. The
highes alue o he iscous s ess and he lowes alue o he Reynolds s ess o he pisoFoam sol e confi ms ha his sol e is leas
able o de elop u bulence on his specific mesh. The oo -mean-squa e (RMS) eloci ies a e depic ed in Fig. 12c, which shows ha
he alues o 𝑅𝑀𝑆(𝑢𝑥)a e oo high and he RMS eloci ies in o he di ec ions a e oo low, he leas accu a e solu ion again gi en
by pisoFoam. These inaccu acies can be caused by coa se mesh esolu ion, especially in he s eam-wise di ec ion.
The u bulen kine ic ene gy budge s a e gi en in Figs. 13a h ough 13e. In gene al, he implemen ed sol e s show beha iou
esembling he e e ence solu ion, howe e , because o he coa se mesh hey ha e no a ained an accu a e solu ion. The lowe
u bulence o pisoFoam can be seen in hese figu es by he lowe absolu e alues o he peaks in p oduc ion, anspo , iscous
diffusion and p essu e diffusion. The dissipa ion budge also shows lowe absolu e alues, and a lesse exp ession o he cha ac e is ic
bump in he line a ound 𝑦+=15. Fo all sol e s, he loca ion o he peak in p oduc ion seems o occu mo e o less a he loca ion
whe e he Reynolds shea s ess and he iscous s ess a e equal, as is expec ed [56]. The alue o he peak is sligh ly oo high o oo
low o he implemen ed me hods, wi h he 𝜃0me hod ha ing he lowes p oduc ion o he h ee. pisoFoam has he lowes p oduc ion,
in line wi h wha was seen be o e. The inaccu acies in p oduc ion seem o be mainly compensa ed by he iscous diffusion a ound
he peak and he dissipa ion close o he wall. The unde -es ima ion o he anspo and iscous diffusion e ms u he away om
he wall also seem o be compensa ed by he dissipa ion, o which he alues a e oo high in his egion. The cha ac e is ic bump
in he dissipa ion e m is no well esol ed on his mesh, as he p ofiles a e oo smoo h. Finally, he p essu e diffusion e m shows
ela i ely la ge de ia ions om he e e ence alue, since p essu e diffusion is small in gene al, he de ia ions appea mo e isibly.
In gene al, he implemen ed sol e s show simila esul s, and he de ia ions om he e e ence alues a e expec ed o be caused by
he coa se mesh esolu ion.
Finally, he esul s o he 𝑋180 mesh using he Runge-Ku a 3 empo al in eg a ion me hod a e shown in Figs. 14a h ough 14c.
Fig. 14a shows ha he mean s eam-wise eloci y is e y accu a ely p edic ed, indica ing ha he u bulence is p ope ly de eloped.
The RMS eloci ies also depic his, al hough a sligh de ia ion can s ill be seen in he alues o 𝑅𝑀𝑆(𝑢𝑥). This can be explained by
some unde - esolu ion in he s eam-wise di ec ion, since he mesh is no o ull di ec nume ical simula ion (DNS) quali y. Finally,
Fig. 14c shows ha he implemen ed sol e s con e ge o one ano he and lie e y closely o he e e ence alues. Con e gence o he
e e ence is no ye a ained o all alues, bu he igh end is isible and sugges s ha addi ional mesh efinemen would u he
inc ease con e gence o he esul s. Especially mesh efinemen in he s eam-wise di ec ion is needed o imp o e he 𝑅𝑀𝑆(𝑢𝑥)
alues.
O e all, wi h mesh efinemen , he implemen ed sol e s seem o con e ge o he same solu ion, and checke boa ding seems o
diminish. The ange in which he 𝜃𝑑𝑦 me hod ope a es is he e o e also smalle . On coa se meshes and in lamina flows howe e ,
he me hod s ill has significan benefi s and i s abili y o balance dynamically be ween adding nume ical dissipa ion o allowing
checke boa ding. In flows ha show ansi ion om lamina flow o u bulence o e ime, he me hod should pe o m well wi h
espec o he o he me hods. Finally, i local quan ifica ion and local balancing is de eloped, flows wi h u bulence ansi ion a eas
could benefi om his me hod as well.
5. Concluding ema ks
Nume ical analysis
This wo k has p esen ed an elabo a e o e iew o he diffe en me hods in ob aining and a oiding checke boa d oscilla ions
when using he p ojec ion me hod o sol e incomp essible flows on colloca ed g ids, by using he symme y-p ese ing amewo k
o which he de i a ions and implemen a ions o he disc e e ope a o s we e also explained in de ail. Dis inc ion be ween usage o a
wide-s encil o a compac -s encil Laplacian combined wi h a di ec (DIM) o weigh ed (WIM) in e pola ion me hod lead o a o al o
ou me hods. The wide-s encil DIM is he me hod o which one would a i e when using a cen al diffe encing scheme, which is e y
p one o checke boa ding and he e o e usually a oided. The compac -s encil DIM in oduces an unnecessa y dissipa i e smoo hing

Jou nal o Compu a ional Physics 521 (2025) 113537
17
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Fig. 12. (a) Mean s eam-wise eloci y on he 𝑋40 mesh. (b) P ofiles o he iscous and Reynolds shea s ess on he 𝑋40 mesh, showing ha pisoFoam is no able o
ully de elop u bulence. (c) RMS eloci ies on he 𝑋40 mesh. Influence o he empo al in eg a ion scheme is shown by plo ing he 𝜃0me hod using he RK3 scheme.
o he eloci y fields, and is he e o e also undesi able. Be ween he WIMs, he wide-s encil WIM leads o an e o in he con ec i e
e m, which is no s ic ly-dissipa i e and can lead o uns able solu ions. The p essu e e o o he compac -s encil WIM is a o able
because i s o de can be educed by employing he Van Kan me hod. Mo eo e , he compac -s encil Laplacian is compu a ionally
a o able.
Quan i ying checke boa ding
The p oblem o he lack o defini ion o checke boa ding in li e a u e was add essed and a solu ion o his p oblem was sugges ed
by quan i ying he phenomenon as 𝐶𝑐𝑏 =1 −𝐩𝑇
𝑐(𝐿𝑐)𝐩𝑐
𝐩𝑇
𝑐𝐿𝐩𝑐
, which is a physics-based, global, non-dimensional, no malised coefficien . The
basis o his coefficien is ied o he diffe ence be ween he wide-s encil and compac -s encil Laplacian ope a o s. This diffe ence
is linked o he undamen al p oblem o he colloca ed g id a angemen , which is una oidably accompanied by in e pola ions, o
which he e canno be an in e se such ha he o iginal field is ees ablished, i.e. Γ𝑐𝑠Γ𝑠𝑐 ≈𝐼𝑠. This diffe ence also o ms he basis o
he nume ical e o s ha a e in oduced by he diffe en ac ional s ep me hods o coun e -ac checke boa d oscilla ions.
Applica ion o he checke boa d coefficien and pe o mance
This coefficien was es ed using lamina and u bulen flows, and compa ed o quali a i e and o he quan i a i e assessmen s,
in which i showed o gi e an in ui i e and seemingly co ec es ima ion o he le els o checke boa ding. An example o a possible
usage o he coefficien was also in oduced, by employing i as a pa ame e ha g adually includes a p essu e p edic o in o he
momen um p edic o equa ion. Since he inclusion o a p essu e p edic o i sel can cause checke boa ding, his effec i ely es ablishes
a dynamic balancing be ween checke boa ding and nume ical dissipa ion h ough nega i e eedback. In he pe o med es cases, his
sol e showed i s possibili y in achie ing his balance wi hou equi ing any use inpu . By doing so i was able o achie e low le els
o nume ical dissipa ion in cases wi hou any oscilla ions, whe eas i was also able o educe he amoun o checke boa ding in mo e
Jou nal o Compu a ional Physics 521 (2025) 113537
18
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Fig. 13. (a) Tu bulen kine ic ene gy p oduc ion on he 𝑋40 mesh. (b) Tu bulen kine ic ene gy anspo on he 𝑋40 mesh. (c) Viscous diffusion on he 𝑋40 mesh. (d)
Dissipa ion on he 𝑋40 mesh. (e) P essu e diffusion on he 𝑋40 mesh.
Jou nal o Compu a ional Physics 521 (2025) 113537
19
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Fig. 14. (a) Mean s eam-wise eloci y on he 𝑋180 mesh. (b) RMS eloci ies on he 𝑋180 mesh. (c) Tu bulen kine ic ene gy budge s on he 𝑋180 mesh.
challenging cases. In he ansien es case, he sol e showed no significan diffe ence in accu acy, while also showing con e gence
o he same solu ion as he Cho in and Van Kan me hods.
Scope o applicabili y
Cu en ly, he me hod has i s g ea es use o lamina flows, flows on coa se meshes, cases wi h ansi ion be ween lamina
and u bulen flows o mo e gene ally, simula ions ha a e un wi hou knowing he esul ing flow a p io i. Fo efined meshes and
u bulen flows, he checke boa ding seems o diminish in gene al, and he ope a ing window o he dynamic sol e becomes smalle ,
while he accu acy o he esul s a e simila o he Van Kan and Cho in me hods, making he me hod less use ul.
C i ical conside a ions
One d awback o his me hod o conside is he loss o i s second-o de ela ion be ween ime-s ep size and p essu e e o , o
which i could be in e es ing o conside a diffe en eedback mechanism ins ead o he p oposed linea eedback, such ha he
me hod s ops u ilising a p essu e p edic o only a ce ain le els o checke boa ding. Ano he conside a ion o he quan ifica ion
me hod is ha he le el o checke boa ding dec eases g ea ly when he mesh is efined. The e o e, a single measu emen o he le el
o checke boa ding wi hou any con ex can be misleading o p esen li le use. Howe e , i s ill emains as a alid quan ifica ion
me hod whe e no o he me hods a e widely-known o -used.
Fu u e wo k
An in e es ing sugges ion o u u e wo ks would be he measu emen o checke boa ding on a local o e en cell le el. This could
p o ide mo e insigh in o p oblema ic a eas, sugges ing mesh efinemen o e en eedback mechanisms ha ac on a local le el, o
example by local inclusion o a p essu e p edic o o hyb id usage o a compac - and wide-s encil scheme depending on he le el
o checke boa ding. This could hen also make his me hod use ul in cases ha ha e lamina o u bulen flow ansi ional egions.
Ano he po en ial a ea o in e es is using 𝐶𝑐𝑏 in o he eedback loops, such as o he de e mina ion o he ime-s ep size, since i
Jou nal o Compu a ional Physics 521 (2025) 113537
20
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
affec s he le el o checke boa ding. Howe e , since a dec easing ime-s ep size inc eases he le els o checke boa ding, he nega i e
eedback would end o inc ease he ime-s ep. Gene ally speaking, a la ge ime-s ep is al eady sough a e , bu is limi ed by o he
ac o s, such as he CFL condi ion.
CRediT au ho ship con ibu ion s a emen
J.A. Hopman: W i ing – e iew & edi ing, W i ing – o iginal d a , Visualiza ion, Valida ion, So wa e, Me hodology, In es iga-
ion, Fo mal analysis, Da a cu a ion, Concep ualiza ion. D. San os: W i ing – e iew & edi ing, Me hodology, Fo mal analysis. À.
Alsal i-Baldellou: W i ing – e iew & edi ing, Me hodology, Fo mal analysis. J. Rigola: Supe ision, P ojec adminis a ion, Fund-
ing acquisi ion. F.X. T ias: W i ing – e iew & edi ing, Supe ision, P ojec adminis a ion, Me hodology, In es iga ion, Funding
acquisi ion, Fo mal analysis, Concep ualiza ion.
Decla a ion o compe ing in e es
The au ho s decla e ha hey ha e no known compe ing financial in e es s o pe sonal ela ionships ha could ha e appea ed o
influence he wo k epo ed in his pape .
Appendix A. A no e on colloca ed g adien s
Using he defini ion 𝐺𝛾
𝑐=Γ
𝛾
𝑠𝑐𝐺=Ω
−1Γ𝛾𝑇
𝑐𝑠 Ω𝑠𝐺can be qui e cumbe some in nume ical analyses and in e ms o implemen a ion
in o codes. In his sec ion his e m is ew i en o a mo e in ui i e no a ion which can be easily implemen ed in o mos codes. The
s a ing poin o his de i a ion is a gene al colloca ed g adien , 𝐺𝛾
𝑐, o which he in e pola ion me hod is no ye specified. Using
he defini ions om 2, he ollowing de i a ion is made:
𝐺𝛾
𝑐=Γ
𝛾
𝑠𝑐𝐺
=Ω
−1Γ𝛾𝑇
𝑐𝑠 Ω𝑠𝐺
=−Ω
−1 [𝐼3⊗Π𝛾𝑇
𝑐𝑠 ]𝑁𝑇
𝑓𝑀𝑇
=Ω
−1 [𝐼3⊗(𝑇𝑓𝑜𝑊𝛾
𝑜+𝑇𝑓𝑛𝑊𝛾
𝑛)]𝑆𝑇
𝑓(𝑇𝑇
𝑓𝑛 −𝑇𝑇
𝑓𝑜)
=Ω
−1 ⎛⎜⎜⎜⎜⎝
−[𝐼3⊗(𝑇𝑓𝑜 (𝐼𝑚−𝑊𝛾
𝑛))]𝑆𝑇
𝑓𝑇𝑇
𝑓𝑜
[𝐼3⊗(𝑇𝑓𝑜𝑊𝛾
𝑜)] 𝑆𝑇
𝑓𝑇𝑇
𝑓𝑛
−[𝐼3⊗(𝑇𝑓𝑛𝑊𝛾
𝑛)] 𝑆𝑇
𝑓𝑇𝑇
𝑓𝑜
[𝐼3⊗(𝑇𝑓𝑛 (𝐼𝑚−𝑊𝛾
𝑜))]𝑆𝑇
𝑓𝑇𝑇
𝑓𝑛
⎞⎟⎟⎟⎟⎠
=Ω
−1 [𝐼3⊗(𝑇𝑓𝑜 −𝑇𝑓𝑛)]𝑆𝑇
𝑓(𝑊𝛾
𝑛𝑇𝑇
𝑓𝑜 +𝑊𝛾
𝑜𝑇𝑛𝑓 )
=Ω
−1 [𝐼3⊗𝑀]𝑁𝑇
𝑓(𝑊𝛾
𝑛𝑇𝑇
𝑓𝑜 +𝑊𝛾
𝑜𝑇𝑛𝑓 )
=𝐺𝐺Π𝛾
𝑐𝑠,
(A.1)
whe e a each cell 𝑖 he ollowing simplifica ion is made:
[−[𝐼3⊗𝑇
𝑓𝑜𝐼𝑚]𝑆𝑇
𝑓𝑇𝑇
𝑓𝑜 +[𝐼3⊗𝑇
𝑓𝑛𝐼𝑚]𝑆𝑇
𝑓𝑇𝑇
𝑓𝑛]𝑖=− ∑
𝑓∈𝐹(𝑖)
𝐬𝑓(𝑖)=𝟎,(A.2)
which ollows om he ac ha he sum o all ou wa d-poin ing su ace ec o s o a closed geome y equals ze o, which holds o any
numbe o dimensions. In he final line, wo new ope a o s a e in oduced, he Gauss g adien ope a o , 𝐺𝐺, and he in e se-weigh ed
cell- o- ace in e pola o , Π𝛾
𝑐𝑠. The Gauss g adien can be explained as assigning he no mal di ec ion o each ace o he scala alue
a he ace mul iplied by he ace a ea, summing hese ace- ec o s o he cell-cen e , hen di iding by he cell olume, which is a
common me hod o ake g adien s in he fini e olume me hod. The in e se-weigh ed cell- o- ace in e pola o ope a es exac ly as he
no mal cell- o- ace in e pola o , bu he weigh s o neighbou and owne a e swapped. This has he no ewo hy p ope y ha :
Π𝑀
𝑐𝑠 = Π𝑀
𝑐𝑠 ,Π𝑉
𝑐𝑠 = Π𝐿
𝑐𝑠,Π𝐿
𝑐𝑠 = Π𝑉
𝑐𝑠,(A.3)
such ha finally:
𝐺𝑐=Γ
𝑉
𝑠𝑐𝐺=𝐺𝐺Π𝐿
𝑐𝑠.(A.4)
In conclusion, using he de i a ion p o ided by equa ion (A.1), in combina ion wi h he defini ions o equa ion (A.3), i becomes
easy o ew i e any wide-s encil g adien o a combina ion o a scala in e pola ion and a Gauss g adien , which can be help ul in
nume ical analysis. Mo eo e , as discussed below, i g ea ly simplifies implemen a ion in o nume ical algo i hms.
Jou nal o Compu a ional Physics 521 (2025) 113537
21
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
A.1. A no e on implemen a ion
Aside om simpli ying nume ical analyses, equa ion (A.4)g ea ly simplifies he implemen a ion o he colloca ed g adien in o
code. Fo example, in he popula open-sou ce fini e olume code OpenFOAM, implemen ing 𝐺𝑐as Γ𝑠𝑐𝐺is no s aigh - o wa d,
especially in symme y-p ese ing codes. Fo hese codes, he choice o in e pola ion has o be ha d-coded in o he sol e , o emo e
any deg ee o eedom in choice o in e pola o . Ins ead o w i ing a unc ion o he ace- o-cell in e pola o , Γ𝑠𝑐 , equa ion (A.4) can be
used, esul ing in mo e con en ional unc ions ha a e usually eadily a ailable. Fo example, he symme y-p ese ing OpenFOAM
sol e RKSymFoam uses he ollowing syn ax o ha d-code 𝐺𝑉
𝑐𝐩𝑛
𝑐as 𝐺𝐺Π𝐿
𝑐𝑠𝐩𝑛
𝑐[48]:
cons olVec o Field g adpn
(
:: gaussG ad<scala >::g ad
(
linea <scala >(mesh). in e pola e(pn),
"g adpn"
)
);
which uses he g ad () and in e pola e() unc ions ha a e p o ided in he sou ce code o OpenFOAM.
Appendix B. Ke nel ec o s o a bi a y Ca esian meshes
In his sec ion, a se o ke nel ec o s o a bi a y Ca esian meshes is de i ed o he linea , midpoin and olume ic wide-s encil
Laplacian ope a o s, 𝐿𝐿
𝑐, 𝐿𝑀
𝑐, 𝐿𝑉
𝑐. The de i a ion is shown o he wo-dimensional case. In his case, ou linea ly independen ke nel
ec o s a e equi ed o span he null-space; he cons an ke nel ec o , 𝐩−
𝑐(00,𝛾), he e ical ke nel ec o , 𝐩−
𝑐(10,𝛾), he ho izon al ke nel
ec o , 𝐩−
𝑐(01,𝛾)and finally, he checke ed ke nel ec o , 𝐩−
𝑐(11,𝛾). Since he case o 𝛾=𝑀was gi en in [33], only he cases o 𝛾=𝐿
and 𝛾=𝑉will be de i ed he e. Fo he de i a ions, he ac ha 𝐾𝑒𝑟(𝐺𝛾
𝑐) ∈𝐾𝑒𝑟(𝐿𝛾
𝑐)was used in addi ion o he equali y 𝐺𝛾
𝑐=𝐺𝐺Π𝛾
𝑐𝑠,
which was de i ed in equa ion (A.1). A e applying hese equali ies, a se o ec o s, 𝜙
𝜙
𝜙𝑐needs o be ound such ha :
𝐺𝛾
𝑐𝜙
𝜙
𝜙𝑐=𝐺𝐺Π𝛾
𝑐𝑠𝜙
𝜙
𝜙𝑐=𝟎𝑐.(B.1)
Since he cons an ke nel ec o is a i ial solu ion o his equa ion, he ke nel ec o ha is conside ed fi s is he e ical linea
ke nel ec o , 𝐩−
𝑐(10,𝐿). I he alues o his ec o a e chosen as:
[𝐩−
𝑐(10,𝐿)]𝑖,𝑗 =(−1)
𝑖(Δ𝑥𝑖)−1 ,(B.2)
a e ically s iped pa e n will be ob ained, as seen in Fig. B.15. Since he e ically neighbou ing cells ha e equal alues, he alues
a aces 𝑠and 𝑛a e i ially gi en by:
[Π𝐿
𝑐𝑠𝐩−
𝑐(10,𝐿)]𝑠=[Π𝐿
𝑐𝑠𝐩−
𝑐(10,𝐿)]𝑛=−1
Δ𝑥1
.(B.3)
Gi en he ac ha Π𝐿
𝑐𝑠 =Π
𝑉
𝑐𝑠, he alue a ace 𝑤is calcula ed as:
[Π𝐿
𝑐𝑠𝐩−
𝑐(10,𝐿)]𝑤=[Π𝑉
𝑐𝑠𝐩−
𝑐(10,𝐿)]𝑤=Δ𝑥0𝜙𝑊+Δ𝑥1𝜙𝑃
Δ𝑥0+Δ𝑥1
=Δ𝑥0∕Δ𝑥0−Δ𝑥1∕Δ𝑥1
Δ𝑥0+Δ𝑥1
=0.
(B.4)
Simila ly, [Π𝐿
𝑐𝑠𝐩−
𝑐(10,𝐿)]𝑒=0. Since his leads o opposing ace pai s wi h equal alues, i becomes immedia ely e iden ha
[𝐺𝐺Π𝐿
𝑐𝑠𝐩−
𝑐(10,𝐿)]𝑃=0. This equali y holds h oughou he whole mesh and he e o e he ec o gi en by equa ion (B.2)lies on he
ke nel o 𝐿𝐿
𝑐. The ho izon al linea ke nel ec o is de i ed simila ly by swapping he axes. This will be shown by de i ing he alues
o he ho izon al olume ic ke nel ec o , 𝐩−
𝑐(01,𝑉 ), in which, addi ionally, he alues a each cell a e eplaced by hei in e se. This
ke nel ec o has i s alues gi en by:
[𝐩−
𝑐(01,𝑉 )]𝑖,𝑗 =(−1)
𝑗Δ𝑦𝑗,(B.5)
so ha :

Jou nal o Compu a ional Physics 521 (2025) 113537
22
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Fig. B.15. Ve ical linea ke nel ec o , 𝐩−
𝑐(10,𝑉 ), wi h 𝐼=0,𝐽=0and 𝛾=𝑉.
[Π𝑉
𝑐𝑠𝐩−
𝑐(01,𝑉 )]𝑤=[Π𝑉
𝑐𝑠𝐩−
𝑐(01,𝑉 )]𝑒=−Δ𝑦1,(B.6)
[Π𝑉
𝑐𝑠𝐩−
𝑐(01,𝑉 )𝑐]𝑠=[Π𝐿
𝑐𝑠𝐩−
𝑐(01,𝑉 )]𝑠=Δ𝑦1𝜙𝑆+Δ𝑦0𝜙𝑃
Δ𝑦0+Δ𝑦1
(B.7)
=Δ𝑦1Δ𝑦0−Δ𝑦0Δ𝑦1
Δ𝑦0+Δ𝑦1
=0,
[Π𝑉
𝑐𝑠𝐩−
𝑐(01,𝑉 )]𝑛=0,(B.8)
which also leads o equal opposi e ace pai alues. The checke ed ke nel ec o s a e sligh ly mo e difficul o see di ec ly. Fo 𝐿𝑉
𝑐,
he alues o ke nel ec o 𝐩−
𝑐(11,𝑉 )a e gi en by:
[𝐩−
𝑐(11,𝑉 )]𝑖,𝑗 =(−1)
𝑖+𝑗Δ𝑥𝑖Δ𝑦𝑗,(B.9)
leading o:
[Π𝑉
𝑐𝑠𝜙
𝜙
𝜙𝑐]𝑤=[Π𝐿
𝑐𝑠𝜙
𝜙
𝜙𝑐]𝑤=Δ𝑥1𝜙𝑊+Δ𝑥0𝜙𝑃
Δ𝑥0+Δ𝑥1
(B.10)
=Δ𝑥1(Δ𝑥0Δ𝑦1)−Δ𝑥0(Δ𝑥1Δ𝑦1)
Δ𝑥0+Δ𝑥1
=0,
[Π𝑉
𝑐𝑠𝜙
𝜙
𝜙𝑐]𝑒=[Π𝑉
𝑐𝑠𝜙
𝜙
𝜙𝑐]𝑠=[Π𝑉
𝑐𝑠𝜙
𝜙
𝜙𝑐]𝑛=0.(B.11)
Which gi es simila esul s o he checke ed olume ic ke nel ec o . Ex ending hese a gumen s o all combina ions o pa e ns
and in e pola o s and o he h ee-dimensional case, a ull se o linea ly independen ke nel ec o s o a bi a y Ca esian meshes
can be de i ed. Fo he wide-s encil Laplacian ope a o s 𝐿𝛾
𝑐 his se is gi en by:
[𝐩−
𝑐(𝐼𝐽𝐾,𝛾)]𝑖,𝑗,𝑘 =(−1)
𝑖𝐼+𝑗𝐽+𝑘𝐾 ([Δ𝑥]𝐼
𝑖[Δ𝑦]𝐽
𝑗[Δ𝑧]𝐾
𝑘)𝛼
,∀𝐼,𝐽,𝐾 ∈{0,1},
𝛼=⎧
⎪
⎨
⎪
⎩
−1 i 𝛾=𝐿
0i 𝛾=𝑀
1i 𝛾=𝑉
.
(B.12)
In conclusion, a se o ec o s ha spans he ke nel o 𝐿𝛾
𝑐can always be de i ed o any Ca esian mesh and o any o he discussed
in e pola o s, using equa ion (B.12). Since he ke nel is some imes used in desc ibing he checke boa d phenomenon, using equa ion
(B.12)sa es a lo o compu a ional effo compa ed o pe o ming a singula alue decomposi ion.
Da a a ailabili y
Da a will be made a ailable on eques .
Jou nal o Compu a ional Physics 521 (2025) 113537
23
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
Re e ences
[1] S.V. Pa anka , Nume ical Hea T ans e and Fluid Flow, Hemisphe e Publishing Company, 1980.
[2] P. Wesseling, P inciples o Compu a ional Fluid Dynamics, ol. 29, Sp inge Science & Business Media, 2009.
[3] J.H. Fe zige , M. Pe ić, Compu a ional Me hods o Fluid Dynamics, Sp inge -Ve lag, Be lin, Heidelbe g, New Yo k, 2002.
[4] H.K. Ve s eeg, W. Malalaseke a, An In oduc ion o Compu a ional Fluid Dynamics: The Fini e Volume Me hod, Pea son Educa ion, 2007.
[5] A.J. Cho in, Nume ical solu ion o he Na ie –S okes equa ions, Ma h. Compu . 22 (104) (1968) 745–762.
[6] R. Temam, Su l’app oxima ion de la solu ion des equa ions de Na ie —S okes pa la me hode des pas ac ionnai es ii, A ch. Ra ion. Mech. Anal. 33 (1969)
377–385.
[7] S.V. Pa anka , D.B. Spalding, A calcula ion p ocedu e o hea , mass and momen um ans e in h ee-dimensional pa abolic flows, in: Nume ical P edic ion o
Flow, Hea T ans e , Tu bulence and Combus ion, Else ie , 1983, pp. 54–73.
[8] J. Kim, P. Moin, Applica ion o a ac ional-s ep me hod o incomp essible Na ie –S okes equa ions, J. Compu . Phys. 59 (2) (1985) 308–323.
[9] J. Van Kan, A second-o de accu a e p essu e-co ec ion scheme o iscous incomp essible flow, SIAM J. Sci. S a . Compu . 7(3) (1986) 870–891.
[10] F.H. Ha low, J.E. Welch, Nume ical calcula ion o ime-dependen iscous incomp essible flow o fluid wi h ee su ace, Phys. Fluids 8 (12) (1965) 2182–2189.
[11] B. Pe o , Conse a ion p ope ies o uns uc u ed s agge ed mesh schemes, J. Compu . Phys. 159 (1) (2000) 58–89.
[12] A. Pascau, Cell ace eloci y al e na i es in a s uc u ed coloca ed g id o he uns eady Na ie –S okes equa ions, In . J. Nume . Me hods Fluids 65 (7) (2011)
812–833.
[13] C.M. Rhie, W.L. Chow, Nume ical s udy o he u bulen flow pas an ai oil wi h ailing edge sepa a ion, AIAA J. 21 (11) (1983) 1525–1532.
[14] M. Pe ić, R. Kessle , G. Scheue e , Compa ison o fini e- olume nume ical me hods wi h s agge ed and coloca ed g ids, Compu . Fluids 16 (4) (1988) 389–403.
[15] W.Z. Shen, J.A. Michelsen, J.N. Sø ensen, Imp o ed Rhie-Chow in e pola ion o uns eady flow compu a ions, AIAA J. 39 (12) (2001) 2406–2409.
[16] S.K. Choi, No e on he use o momen um in e pola ion me hod o uns eady flows, Nume . Hea T ans ., Pa A, Appl. 36 (5) (1999) 545–550.
[17] B. Yu, Y. Kawaguchi, W.Q. Tao, H. Ozoe, Checke boa d p essu e p edic ions due o he unde elaxa ion ac o and ime s ep size o a nons agge ed g id wi h
momen um in e pola ion me hod, Nume . Hea T ans ., Pa B, Fundam. 41 (1) (2002) 85–94.
[18] B. Yu, W.-Q. Tao, J.-J. Wei, Y. Kawaguchi, T. Tagawa, H. Ozoe, Discussion on momen um in e pola ion me hod o colloca ed g ids o incomp essible flow,
Nume . Hea T ans ., Pa B, Fundam. 42 (2) (2002) 141–166.
[19] S. Majumda , Role o unde elaxa ion in momen um in e pola ion o calcula ion o flow wi h nons agge ed g ids, Nume . Hea T ans . 13 (1) (1988) 125–132.
[20] T.F. Mille , F.W. Schmid , Use o a p essu e-weigh ed in e pola ion me hod o he solu ion o he incomp essible Na ie –S okes equa ions on a nons agge ed
g id sys em, Nume . Hea T ans . 14 (2) (1988) 213–233.
[21] B. Yu, W.Q. Tao, J.J. Wei, Y. Kawaguchi, T. Tagawa, H. Ozoe, Discussion on momen um in e pola ion me hod o colloca ed g ids o incomp essible flow, Nume .
Hea T ans ., Pa B, Fundam. 42 (2) (2002) 141–166.
[22] J.L. Gue mond, P. Mine , J. Shen, An o e iew o p ojec ion me hods o incomp essible flows, Compu . Me hods Appl. Mech. Eng. 195 (44–47) (2006)
6011–6045.
[23] S. Zhang, X. Zhao, S. Bayyuk, Gene alized o mula ions o he Rhie–Chow in e pola ion, J. Compu . Phys. 258 (2014) 880–914.
[24] P. Ba holomew, F. Denne , M.H. Abdol-azis, A. Ma quis, B.G.M. Van Wachem, Unified o mula ion o he momen um-weigh ed in e pola ion o colloca ed
a iable a angemen s, J. Compu . Phys. 375 (2018) 177–208.
[25] F.N. Fel en, T.S. Lund, Kine ic ene gy conse a ion issues associa ed wi h he colloca ed mesh scheme o incomp essible flow, J. Compu . Phys. 215 (2) (2006)
465–484.
[26] J.E. Ma son, An In oduc ion o Ansys Fluen 2023, 1s edi ion, SDC Publica ions, 2023.
[27] S.D.I. So wa e, Simcen e STAR-CCM+ Use Guide . 2306, 2023.
[28] C. G eenshields, OpenFOAM 11 Use Guide, London, UK, h ps://doc .c d .di ec /open oam /use -guide - 11, 2023.
[29] F. A chambeau, N. Méchi oua, M. Sakiz, Code sa u ne: a fini e olume code o he compu a ion o u bulen incomp essible flows-indus ial applica ions, In . J.
Fini e Vol. 1(1) (2004).
[30] R.W. Ve s appen, A.E. Veldman, Symme y-p ese ing disc e iza ion o u bulen flow, J. Compu . Phys. 187 (1) (2003) 343–368.
[31] F.X. T ias, O. Lehmkuhl, A. Oli a, C.D. Pé ez-Sega a, R.W. Ve s appen, Symme y-p ese ing disc e iza ion o Na ie –S okes equa ions on colloca ed uns uc u ed
g ids, J. Compu . Phys. 258 (2014) 246–267.
[32] E.M. Komen, J.A. Hopman, E.M. F ede ix, F.X. T ias, R.W. Ve s appen, A symme y-p ese ing second-o de ime-accu a e piso-based me hod, Compu . Fluids
225 (2021) 104979.
[33] J. La sson, G. Iacca ino, e al., A co-loca ed incomp essible Na ie –S okes sol e wi h exac mass, momen um and kine ic ene gy conse a ion in he in iscid
limi , J. Compu . Phys. 229 (12) (2010) 4425–4430.
[34] J.A. Hopman, F. T ias, J. Rigola, On a conse a i e solu ion o checke boa ding: examining he disc e e laplacian ke nel using mesh connec i i y, in: ERCOFTAC
Wo kshop Di ec and La ge Eddy Simula ion, Sp inge , 2023, pp. 306–311.
[35] J.A. Hopman, F.X. T ias Miquel, J. Rigola Se ano, Symme y-p ese ing disc e isa ion me hods o magne ohyd odynamics, in: Wo ld Cong ess in Compu a ional
Mechanics and ECCOMAS Cong ess, Oslo, No way, 2022.
[36] D. San os, F.X. T ias, G. Colome , C.D. Pé ez-Sega a, An ene gy-p ese ing uncondi ionally s able ac ional s ep me hod on colloca ed g ids, in: Wo ld Cong ess
in Compu a ional Mechanics and ECCOMAS Cong ess, Oslo, No way, 2022.
[37] D. San os, F.X. T ias, J.A. Hopman, C.D. Pé ez-Sega a, P essu e- eloci y coupling on uns uc u ed colloca ed g ids: econciling s abili y and ene gy-conse a ion,
in: P oceedings o he Ten h In e na ional Symposium on Tu bulence, Hea and Mass T ans e , Rome, I aly, 2023, pp. 259–262.
[38] D. San os, J.A. Hopman, C.D. Pé ez-Sega a, F.X. T ias, On an ene gy-p ese ing uncondi ionally s able p ojec ion me hod on colloca ed uns uc u ed g ids,
A ailable a SSRN: h ps://doi .o g /10 .2139 /ss n .4788673, 2024.
[39] N.N. Yanenko, The Me hod o F ac ional S eps, Sp inge -Ve lag, New Yo k, 1971.
[40] A.J. Cho in, A nume ical me hod o sol ing incomp essible iscous flow p oblems, J. Compu . Phys. 2(1) (1967) 12–26.
[41] J.A. Hopman, A. Alsal i-Baldellou, F.X. T ias, J. Rigola, On a conse a i e solu ion o checke boa ding: examining he causes o non-physical p essu e modes, in:
P oceedings o he 14 h In e na ional ERCOFTAC Symposium on Enginee ing Tu bulence Modelling and Measu emen s, Ba celona, Spain, 2023, pp. 599–603.
[42] A.E.P. Veldman, A gene al condi ion o kine ic-ene gy p ese ing disc e iza ion o flow anspo equa ions, J. Compu . Phys. 398 (2019) 108894.
[43] M.M. La maei, J. Behzadi, T.F. Mahdi, T ea men o checke boa d p essu e in he colloca ed uns uc u ed fini e- olume scheme, Nume . Hea T ans ., Pa B,
Fundam. 58 (2) (2010) 121–144.
[44] C.M. Klaij, On he s abiliza ion o fini e olume me hods wi h co-loca ed a iables o incomp essible flow, J. Compu . Phys. 297 (2015) 84–89.
[45] P. Rauwoens, J. Vie endeels, B. Me ci, A solu ion o he odd-e en decoupling p oblem in p essu e-co ec ion algo i hms o a iable densi y flows, J. Compu .
Phys. 227 (1) (2007) 79–99.
[46] A.W. Da e, Fluid dynamical iew o p essu e checke boa ding p oblem and smoo hing p essu e co ec ion on meshes wi h coloca ed a iables, In . J. Hea Mass
T ans . 46 (25) (2003) 4885–4898.
[47] G.H. Golub, C.F. Van Loan, Ma ix Compu a ions, Johns Hopkins Uni e si y P ess, 1996.
[48] J.A. Hopman, E.M.A. F ede ix, RKSymFoam Gi Hub page, h ps://gi hub .com /janneshopman /RKSymFoam, 2023.
[49] V. Vuo inen, J.-P. Keskinen, C. Duwig, B.J. Boe sma, On he implemen a ion o low-dissipa i e Runge–Ku a p ojec ion me hods o ime dependen flows using
OpenFOAM®, Compu . Fluids 93 (2014) 153–163.
Jou nal o Compu a ional Physics 521 (2025) 113537
24
J.A. Hopman, D. San os, À. Alsal i-Baldellou e al.
[50] G.I. Taylo , A.E. G een, Mechanism o he p oduc ion o small eddies om la ge ones, P oc. R. Soc. Lond. Se . A, Ma h. Phys. Sci. 158 (895) (1937) 499–521.
[51] A.W. V eman, J.G. Kue en, Compa ison o di ec nume ical simula ion da abases o u bulen channel flow a 𝑅𝑒𝜏= 180, Phys. Fluids 26 (1) (2014).
[52] E.M.J. Komen, E.M.A. F ede ix, T.H.J. Coppen, V. D’alessand o, J.G.M. Kue en, Analysis o he nume ical dissipa ion a e o diffe en Runge-Ku a and eloci y
in e pola ion me hods in an uns uc u ed colloca ed fini e olume me hod in open oam, Compu . Phys. Comun. (2020).
[53] H. Zhang, F.X. T ias, A. Go obe s, Y. Tan, A. Oli a, Di ec nume ical simula ion o a ully de eloped u bulen squa e duc flow up o 𝑅𝑒𝜏= 1200, In . J. Hea
Fluid Flow 54 (2015) 258–267.
[54] P. Du bin, B. Rei , S a is ical Theo y and Modeling o Tu bulen Flows, Wiley & Sons, L d, 2011.
[55] J.A. Hopman, unTimeChannelBudge s Gi Hub page, 2023, h ps://gi hub .com /janneshopman / unTimeChannelBudge s.
[56] S.B. Pope, Tu bulen Flows, Camb idge Uni e si y P ess, 2000.