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
A high-o de imme sed bounda y me hod o app oxima e flow
p oblems in domains wi h cu ed bounda ies
S. Colomboa, ,∗, G. Rubioa,b, J. Kou c, E. Vale oa,b, R. Codinad,e, E. Fe e a,b
aETSIAE-UPM-School o Ae onau ics, Uni e sidad Poli écnica de Mad id, Mad id, Spain
bCen e o Compu a ional Simula ion, Uni e sidad Poli écnica de Mad id, Mad id, Spain
cSchool o Ae onau ics, No hwes e n Poly echnical Uni e si y, Xi’an, China
dUni e si a Poli ècnica de Ca alunya, Ba celona, Spain
eCen e In e nacional de Mè odes Numè ics en Enginye ia (CIMNE), Ba celona, Spain
A R T I C L E I N F O A B S T R A C T
Keywo ds:
Imme sed bounda y me hod
Cu ed bounda y condi ions
High-o de h/p sol e s
Discon inuous Gale kin
Ho ses3D
High-o de ℎ∕𝑝sol e s in compu a ional fluid dynamics offe scalabili y, efficiency, and supe io
e o educ ion compa ed o adi ional low-o de me hods. Imme sed bounda y me hods
elimina e he need o body-fi ed meshes bu o en deg ade he o de o he solu ion nea
bounda ies, which can damage he o e all accu acy o he high-o de sol e . This pape p esen s
a new app oach o impose bounda y condi ions in high-o de fini e elemen o fini e olume flow
sol e s ha e ain high-o de 𝑃+1con e gence, whe e 𝑃is he polynomial o de . Fu he mo e,
he me hodology akes in o accoun cu ed bounda y condi ions wi hou loss in accu acy. I
in oduces a su oga e bounda y ha elimina es ins abili ies due o badly cu elemen s.
We es he me hodology using a high-o de discon inuous Gale kin amewo k o sol e pu ely
ellip ic p oblems and he comp essible Na ie -S okes equa ions (2D and 3D), o show ha we
e ain he o mal o de o con e gence 𝑃+1. Finally, we compa e he esul s wi h a olume
penaliza ion app oach and show ha spu ious p essu e oscilla ions on he imme sed bounda y
a e elimina ed when he p oposed me hodology is used.
Con en s
1. In oduc ion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2. High-o de discon inuous Gale kin sol e s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3. A high-o de imme sed bounda y me hod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.1. Iden ifica ion o he shi ed aces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.2. Compu a ion o he shi ed aces s a e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.3. Di ichle bounda y condi ion o he Poisson equa ion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.4. Neumann bounda y condi ion o he Poisson equa ions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.5. Su ace ep esen a ion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
4. Tes cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4.1. 1D linea con ec ion-diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4.2. 2D Poisson equa ion wi h Di ichle bounda y condi ions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
*Co esponding au ho .
E-mail add ess: [email p o ec ed] (S. Colombo).
h ps://doi.o g/10.1016/j.jcp.2025.113807
Recei ed 17 July 2024; Recei ed in e ised o m 28 Janua y 2025; Accep ed 30 Janua y 2025
Jou nal o Compu a ional Physics 528 (2025) 113807
A ailable online 3 Feb ua y 2025
0021-9991/© 2025 The Au ho s. Published by Else ie Inc. This is an open access a icle unde he CC BY-NC-ND license
( h p://c ea i ecommons.o g/licenses/by-nc-nd/4.0/ ).
S. Colombo, G. Rubio, J. Kou e al.
4.3. 2D Poisson equa ion wi h Neumann bounda y condi ions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.4. S eady hea equa ion and cu ed bounda ies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.5. Na ie -S okes equa ions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.5.1. Two dimensional ci cula cylinde . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.5.2. Th ee dimensional sphe e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
5. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
CRediT au ho ship con ibu ion s a emen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Decla a ion o compe ing in e es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Acknowledgemen s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Appendix A. Volume penaliza ion imme sed bounda y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Appendix B. Bassi-Rebay 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Appendix C. Some conside a ion abou he iden ifica ion o he shi ed aces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Appendix D. Effec o he a ia ion o he pa ame e 𝛼. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Da a a ailabili y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
Re e ences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1. In oduc ion
In he field o compu a ional fluid dynamics, high-o de h/p sol e s (e.g., based on discon inuous Gale kin) ha e eme ged as a
p omising app oach due o hei scalabili y in mode n a chi ec u es [1] and supe io efficiency o a gi en le el o accu acy [2,3].
These me hods equi e an accu a e geome ic desc ip ion o he bounda ies o main ain high-o de accu acy nea hem. This is
accomplished by gene a ing a body-fi ed mesh a ound cu ed walls, whe e he o de o cu a u e used o ep esen he wall is he
same as he o de o he solu ion inside he elemen (i.e., an isopa ame ic ep esen a ion is assumed). Al hough he gene a ion o
linea meshes has eached a high le el o ma u i y, e en o complex geome ies, he same canno be said o cu ilinea meshes
in complex geome ies. In ac , he gene a ion o body-fi ed cu ed meshes is a bo leneck in he simula ion p ocess and emains a
ib an a ea o esea ch; see, o example, [4,5].
Imme sed bounda y me hods (IBMs) (o embedded me hods o unfi ed bounda y me hods) do no equi e body-fi ed meshes.
O iginally in oduced by Peskin [6], IBM simula es he p esence o he geome ies h ough a ificial mechanisms, hus a oiding he
need o body-fi ed cu ed meshes and allowing complex flows using simple Ca esian g ids. IBMs ha e seen a su ge in popula i y
in ecen yea s [7–9], and a ious echniques ha e been de eloped. In pa icula , cu en IBM app oaches include cu -cell [10–12],
di ec o cing [13–16], ghos -cell [17,18], olume penaliza ion [19–22] and i s ex ension o high-o de sol e s and e o analysis
[22–26]. While IBMs ci cum en he need o body-fi ed meshes, hey o en do so a he expense o accu acy nea su aces; see,
o example, he olume penaliza ion app oach in [22,20,27]. O he me hods can main ain accu acy a he cos o he inc eased
complexi y and deg aded condi ioning in oduced by cu cells. In ac , cu cells a e p one o ins abili ies when small cu elemen s
(a he in e sec ion be ween he backg ound mesh and he imme sed bounda y) eme ge and a emp s ha e been made o bypass he
p oblem by in oducing s abiliza ion [28–30].
To add ess hese difficul ies, he fini e- olume communi y has de eloped in e pola ion schemes o eco e second-o de accu acy
nea bounda ies [31–33], while he fini e-elemen communi y has p oposed, among o he s, he shi ed bounda y me hod (SBM)
[34,35] o e ie e highe accu acy nea imme sed bounda ies. He e, he idea is o shi he bounda y, whe e he bounda y condi ions
a e applied, om he ac ual su ace o a su oga e su ace. The condi ions imposed on he su oga e su ace a e compu ed h ough a
Taylo expansion.
In his pape , we p opose a new echnique, based on he ecen wo k by Funada e al. [36], in which ideas om shi ed bounda ies
and in e pola ion om fini e olumes a e combined. The key ideas o he me hod we p opose a e he ollowing. Suppose ha he
flow domain is Ω, wi h bounda y Γand we ha e a mesh o a backg ound domain Ω𝑠 ha co e s Ω. The s eps o ollow a e:
•C ea e a su oga e bounda y Γ𝑠made o elemen bounda ies o he backg ound mesh,’close’ o he physical bounda y Γ. I can be
bo h in e io o ex e io o he compu a ional domain Ω. Wo king wi h his bounda y will elimina e possible nume ical ins abili ies
due o badly cu elemen s.
•Recons uc he unknowns o he p oblem along a s aigh line no mal o he ue bounda y Γ om bo h in e io deg ees o
eedom o a chosen dono cell (elemen ) and he bounda y alues o be p esc ibed, ei he o he unknowns o o hei de i a i es.
Thus, bo h Di ichle and Neumann bounda y condi ions can be conside ed.
•E alua e he econs uc ed unknowns on he su oga e bounda y Γ𝑠and p esc ibe hese alues weakly, o example by modi ying
app op ia ely he flux ec o o he nume ical scheme.
Using high-o de Lag ange in e pola o s, we enable consis en high-o de solu ions. The p oposed econs uc ion can be seen as an
adap a ion o he me hod de eloped o fini e olume in [37,38] o high-o de sol e s o as a high-o de modifica ion o he shi ed
bounda y me hod (as de eloped by he fini e elemen communi y). Fu he mo e, he me hod is applicable o cu ed bounda ies
wi hou loss o accu acy. To he bes o ou knowledge, his is he fi s ime such an app oach has been de eloped o high-o de ℎ∕𝑝
me hods (he e a discon inuous Gale kin me hod). In his wo k we de elop he me hod in he con ex o he comp essible Na ie -
S okes equa ions, desc ibing also i s applica ion o he simple Poisson’s p oblem; ne e heless, he idea is gene al and can be applied
Jou nal o Compu a ional Physics 528 (2025) 113807
2
S. Colombo, G. Rubio, J. Kou e al.
Fig. 1. Discon inuous Gale kin bounda ies include wo s a es u+and u−, which a e o e alua e he nume ical fluxes a he elemen s in e aces.
o o he p oblems as well. This me hod is obus and allows o he eco e y o high-o de accu acy in s aigh and cu ed geome ic
bounda ies. Finally, in his wo k, we ocus on discon inuous Gale kin disc e iza ions. Howe e , he p oposed me hod can easily
accommoda e o he nume ical s a egies (e.g., fini e olumes, o high o de con inuous fini e elemen s) when using weak en o cemen
o bounda y condi ions. Le us also ema k ha he me hod leads o non-symme ic algeb aic p oblems e en o symme ic bounda y
alue p oblems, as he SBM o he pu ely algeb aic app oach p oposed in [39]. Ne e heless, we shall conside explici ime in eg a ion
schemes e en o app oxima e s a iona y p oblems, and his non-symme y will no be an issue.
The emainde o he pape is o ganized as ollows. Sec ion 2p o ides a b ie o e iew o he nume ical amewo k used h ough-
ou he pape and how he bounda y condi ions a e imposed. In Sec ion 3, he inno a i e me hodology is p esen ed and explained.
Sec ion 3.5 gi es some conside a ion abou he disc e iza ion o he body su ace, while Sec ion 4is dedica ed o a ious es cases
o inc easing complexi y. Finally, Sec ion 5p esen s he conclusions.
2. High-o de discon inuous Gale kin sol e s
All wo k desc ibed h oughou he a icle has been implemen ed in he open sou ce Ho ses3D sol e [40], which is a high-o de
spec al discon inuous Gale kin sol e capable o sol ing a a ie y o flow applica ions, including comp essible flows (wi h o wi h-
ou shocks), incomp essible flows, a ious RANS and LES u bulence models, pa icle dynamics, mul iphase flows, and ae oacous ics.
Ho ses3D handles body-fi ed cu ed bounda ies and also includes a olume penaliza ion IBM. Bo h o hese unc ionali ies will be
compa ed wi h he p oposed me hodology. He e, we p o ide only a b ie desc ip ion o he discon inuous Gale kin me hod imple-
men ed in Ho ses3D in he con ex o flow p oblems.
We conside he gene al case o a second-o de ime-dependen pa ial diffe en ial equa ion:
u𝑡+∇
x⋅
F𝑒=∇
x⋅
F𝑣+
S, x∈Ω, 𝑡>0(1)
whe e
uis he s a e ec o o conse a i e a iables,
F𝑒and
F𝑣a e he in iscid and iscous fluxes espec i ely, and
Sis a sou ce e m.
The physical domain is subdi ided in o nono e lapping cu ilinea hexahed al elemen s, 𝑒, which a e geome ically ans o med o
a e e ence elemen , 𝑒𝑙. T ans o ma ion is pe o med using a polynomial ansfini e mapping ha ela es he physical coo dina es
xand he local e e ence coo dina es 𝝃∈[−1,1]𝑑, whe e 𝑑is he spa ial dimension o he p oblem. Addi ionally, he s a e ec o ,
fluxes, and sou ce e m a e app oxima ed by polynomials using a enso p oduc o he one-dimensional Lag ange basis, esul ing in:
𝐽u𝑡+∇
𝝃⋅F𝑒=∇
𝝃⋅F𝑣+𝐽S,(2)
whe e 𝐽is he Jacobian o he ansfini e mapping, ∇𝝃is he diffe en ial ope a o in he e e ence space, Sis a sou ce e m and F𝑒,
F𝑣a e he con a a ian fluxes [41].
The weak o mula ion o he discon inuous Gale kin is ob ained by mul iplying Eq. (2)by a se o es unc ions 𝜙𝑖 o 0≤𝑖≤𝑃
and in eg a ing o e 𝑒𝑙:
∫
𝑒𝑙
𝜙𝑖𝐽u𝑡=−∮
𝜕𝑒𝑙
𝜙𝑖(F𝑒−F𝑣)⋅
n+∫
𝑒𝑙
∇𝝃𝜙𝑖⋅(F𝑒−F𝑣)+∫
𝑒𝑙
𝜙𝑖𝐽S(3)
whe e
nis he uni no mal ec o o he in e ace 𝜕𝑒𝑙; as cus oma y o he discon inuous Gale kin me hods, he es unc ions a e
aken equal o he basis unc ions. Finally, a summa ion is pe o med o e all he elemen s in he mesh, and Eq. (3)becomes:
∑
𝑒𝑙 ∫
𝑒𝑙 [𝜙𝑖𝐽u𝑡−∇
𝝃𝜙𝑖⋅(F𝑒−F𝑣)−𝜙𝑖𝐽S]+∑
𝜕𝑒𝑙 ∮
𝜕𝑒𝑙
𝜙𝑖(F𝑒−F𝑣)=0,(4)
whe e * ep esen s he jump o he fluxes (including no mals) ac oss he in e ace 𝜕𝑒𝑙, in oducing a coupling be ween 𝑒𝑙 and he
neighbo ing elemen s. Specifically, he jump ac oss he in e ace is compu ed h ough he nume ical fluxes.
The in iscid nume ical flux F𝑒=F⋆
𝑒can be compu ed using a upwind app oach F⋆
𝑒→F⋆
𝑒(u+,u−,
n)and is a unc ion o he
s a e o he in e nal in e ace u+and he s a e o he in e ace o he elemen sha ing he ace wi h 𝑒𝑙, u−, see Fig. 1. Simila ly,
he iscous nume ical flux F𝑣=F⋆
𝑣is a unc ion o bo h u+, u−and he g adien o conse a i e a iables (∇𝝃u)+, (∇𝝃u)−: F⋆
𝑣→
F⋆
𝑣(u+,u−,(∇𝝃u)+,(∇𝝃u)−,
n). In his s udy, we employ he Lax-F ied ichs flux o in iscid e ms and use he Bassi-Rebay 1 (BR1)
[42], he In e io -Penal y (IP) [43], and local discon inuous Gale kin (LDG) [44] fluxes o disc e ize iscous e ms. The pa icula
selec ion o fluxes is de ailed o each es case.
In discon inuous Gale kin me hods, bounda y condi ions on walls ( o body-fi ed meshes) a e imposed weakly h ough he nu-
me ical fluxes. The e o e, no explici imposi ion o bounda y condi ions appea s in he o mula ion. In he case o con inuous fini e
Jou nal o Compu a ional Physics 528 (2025) 113807
3
S. Colombo, G. Rubio, J. Kou e al.
Fig. 2. Ci cle embedded in a simple Ca esian mesh; Fig. 2b shows a close up o he su ace. The compu a ional domain is he ex e io o he ci cle.
elemen me hods, a echnique o weakly impose bounda y condi ion needs o be chosen, such as Ni sche’s me hod. In ou case, we
assume a ghos s a e o u−, see Fig. 1. When he elemen 𝑒𝑙 lies on he physical su ace, 𝜕Ω, he in e ace s a e u−and he in e ace
g adien (∇𝝃u)−a e ob ained om he bounda y condi ions o he p oblem conside ed: u−=u𝑏and (∇𝝃u)−=(∇
𝝃u)𝑏. To impose
bounda y condi ions, an accu a e ep esen a ion o he body su ace is essen ial o ob ain he high-o de accu acy o he DG. In ac ,
he geome y is desc ibed h ough he no mal
n ha appea s in he in e ace in eg als o Eq. (3)and ailing o accoun o he bound-
a y cu a u e can comp omise he accu acy. As a consequence, a cu ilinea mesh is equi ed o p ope ly desc ibe he cu ed su ace
o he body. Conce ning he mesh gene a ion p ocess, while c ea ing linea meshes has become p oficien e en o complex shapes,
he same canno be said o cu ilinea meshes. The e o e, he e is conside able in e es in employing high-o de me hods ha do no
equi e a high-o de cu ed mesh. In he ollowing sec ion, we p esen a echnique capable o achie ing high-o de p ecision wi h
linea Ca esian meshes, bypassing he explici cu a u e o he mesh. This is achie ed by using an IBM as he one p oposed nex .
3. A high-o de imme sed bounda y me hod
We aim o de elop an IBM (see Fig. 2) capable o p ese ing he high-o de accu acy o he unde lying scheme, namely he
discon inuous Gale kin me hod. The s anda d app oach o de eloping IBMs in ol es modi ying he equa ions; o ins ance, in he
case o olume penaliza ion, a sou ce e m is added o each deg ee o eedom loca ed inside he body, which is ea ed as a po ous
medium wi h infini e pe meabili y (see Appendix Aand p o ided e e ences). Howe e , his me hod ypically loses accu acy nea
he su ace [22,23]. As an al e na i e, we p opose o iden i y a se o aces (shi ed aces) o ming a closed piecewise line (2D) o a
su ace (3D). This will be he su oga e bounda y, Γ𝑠, on which he s a e is p ope ly modified. This bounda y can be inside he flow
domain, as o he SBM, o ou side i . We a e in e es ed in flow p oblems and in bounda ies defined by solid bodies, and we will
conside Γ𝑠ou side he flow domain and he e o e inside he body ha con o ms he physical bounda y; howe e he al e na i e o
in e io su oga e bounda ies is also possible wi h he p esen ed o mula ion. As we shall see, ou app oach implies an ex apola ion
o he unknowns, whe eas he choice o in e io Γ𝑠would amoun o an in e pola ion.
Specifically, he changed s a e is such ha he bounda y condi ions a e sa isfied on he body su ace and he nume ical fluxes
F⋆
𝑒and F⋆
𝑣(Eq. (4)) on he shi ed aces a e compu ed. In his way, he equi ed s a e on hese aces is no s ongly en o ced bu is
indi ec ly imposed h ough an in e ace flux, and hus weakly sa isfied.
The iden ifica ion o he shi ed aces and he me hod o adjus ing he s a e a e explained in de ail in he ollowing sec ion.
3.1. Iden ifica ion o he shi ed aces
The se o shi ed aces mus sa is y wo cons ain s. Fo he op ion o ex e io su oga e bounda y ha we shall ollow, he fi s
is ha he ace has o be loca ed en i ely inside he body ha con o ms he bounda y; he second is ha he opposi e ace has o be
pa ially ou side he body. Fig. 3shows examples o possible scena ios o 2D p oblems.
To iden i y he shi ed aces, we use a ay- acing echnique [45], which is also used o compu e he no mal ec o s o he poin s
o he physical bounda y; in pa icula , each deg ee o eedom belonging o a gene ic ace is agged when inside he body. Once
he deg ees o eedom a e selec ed, each ace ha has all he deg ees o eedom inside he body (ou side he flow domain) is sa ed
as a shi ed ace (i he second cons ain men ioned abo e is also ulfilled). Two addi ional checks a e equi ed: i , ollowing he
p ocedu e p e iously explained a ace is agged wice as a shi ed ace, i is disca ded; he second check is done o all he elemen s
ha ing mo e han one shi ed ace and one (o a leas one in 3D) ace ha ing one o mo e deg ees o eedom in he fluid egion. Fo
wha conce ns he la e poin , all he shi ed aces belonging o he elemen a e disca ded, whe eas he only one ha is s o ed is he
one in on o he ace whose deg ees o eedom a e all o in pa inside he fluid domain (see Appendix C).
Jou nal o Compu a ional Physics 528 (2025) 113807
4
S. Colombo, G. Rubio, J. Kou e al.
Fig. 3. Diffe en cases o shi ed aces (in blue). The solid body (shaded g ay a ea) is assumed o be on he igh o he body bounda y ( ed). The shi ed aces a e
inside he solid, and he e o e in he ex e io o he flow domain whe e he p oblem is sol ed. The s a e o he shi ed ace is co ec ed. (Fo in e p e a ion o he
colo s in he figu e(s), he eade is e e ed o he web e sion o his a icle.)
Fig. 4. Shi ed aces and deg ees o eedom lying on hese aces. The flow domain is he ex e io o he ci cle.
A Ca esian mesh wi h an imme sed bounda y ep esen ing a ci cula cylinde is depic ed in Fig. 4. The figu e de ails he deg ees
o eedom belonging o he shi ed aces shown in black. The s a e on each o hese deg ees o eedom is modified o ake in o
accoun he p esence o he embedded body.
3.2. Compu a ion o he shi ed aces s a e
Once he se o shi ed aces has been selec ed, he nex c i ical s ep is o compu e he new s a e o be imposed o weakly sa is y
he bounda y condi ions on he body su ace. Each o he shi ed aces can be seen as an in e ace sepa a ing an elemen ully inside
he body (and he e o e ou side he flow domain) and an elemen pa ially inside he body, which is now called a shi ed elemen .
Mo eo e , hey a e cha ac e ized by he usual in e ace s a es u+and u−(see Fig. 1): u+is he s a e coming om he shi ed elemen
(Fig. 5) while u−is he s a e o he neighbo elemen , ully inside he body. Specifically, he la e s a e u−is he one add essed
by he p oposed me hodology and is p ope ly modified so ha he solu ion on he su ace is he equi ed one, i.e., u𝑏. The newly
compu ed s a e u−=u𝑠𝑏 ( he subsc ip (.)𝑠𝑏 s ands o shi ed bounda y) is hen eplaced in he nume ical flux defini ion o ob ain
F⋆
𝑒(u+,u−,
n)→F⋆
𝑒(u+,u𝑠𝑏,
n). Fo con inuous fini e elemen in e pola ions, he alue o u𝑠𝑏 is wha would be weakly p esc ibed on
he su oga e bounda y. In gene al, i is impo an o unde line he ac ha u𝑠𝑏 ≠u𝑏, see Fig. 5.
We now ocus on he compu a ion o he s a e u𝑠𝑏 o be imposed on each o he selec ed aces, which is he c i ical poin o he
o mula ion. While in he SBM u𝑠𝑏 would be ob ained om u𝑏using a Taylo expansion, o compu e he alue u𝑠𝑏, we econs uc he
solu ion on he shi ed ace by using a one-dimensional Lag angian in e pola ion o o de 𝑃along he body no mal
n𝑠, see Fig. 7. No e ha ,
in gene al,
n𝑠is diffe en om he no mal o he e e ence elemen
n. We compu e
n𝑠as ex e io o he body, i.e., poin ing inwa ds
he compu a ional domain.
The 𝑃𝑡ℎ-o de Lag angian polynomial needs he iden ifica ion o 𝑃+1poin s (dono s encil) on which he solu ion is known.
Fig. 6shows how hese nodes a e selec ed along a ay om a poin on he physical bounda y wi h a known no mal. The geome ical
coo dina es o he gene ic 𝑗𝑡ℎ-node a e:
Jou nal o Compu a ional Physics 528 (2025) 113807
5
S. Colombo, G. Rubio, J. Kou e al.
Fig. 5. Mesh egion su ounding he imme sed bounda y. The ed line ep esen s a gene ic body (imme sed bounda y). The g ay pa is he in e io o he body and
shi ed aces a e highligh ed in black. The zoomed image shows he s a es u+and u𝑠𝑏 on one o he shi ed aces. Fo comple eness, Gauss nodes on he ace ( o 𝑃=
2) a e also shown.
Fig. 6. The s encil poin s a e shown along wi h he coo dina es in he local e e ence ame. The local e e ence ame is such ha o 𝜉=−1, 𝑥=𝑥𝑏and o 𝜉=1,
𝑥=𝑥𝑎=+𝛼√𝑑⋅ℎ+𝑑𝐿. I is impo an o no e ha nodes 𝜉1, 𝜉2and 𝜉3a e chosen so ha , on he e e ence leng h 𝑑𝐿, hey co espond o he Gauss nodes 𝜉𝐺
1, 𝜉𝐺
2
and 𝜉𝐺
3. The poin 𝜉0a ising om he in e sec ion be ween he body bounda y and he no mal 𝑛 is pic u ed in g een.
Fig. 7. The in e pola ion elemen 𝑒𝑖𝑛𝑡 is shaded. The s encil is he same as he one o Fig. 6.
Jou nal o Compu a ional Physics 528 (2025) 113807
6
S. Colombo, G. Rubio, J. Kou e al.
⎧
⎪
⎨
⎪
⎩
x𝑗=x𝑠𝑏 +𝜒𝑗
n𝑠
𝜒𝑗=𝐿+(1+𝜉𝐺
𝑗)
2
𝑑𝐿, wi h 𝑗=1,…,𝑃
𝜒0=,
(5)
whe e 𝜉𝐺
𝑗a e he Gauss nodes dis ibu ed on 𝑑𝐿. We defined he segmen 𝑑𝐿 as |𝜉𝐺
0−𝜉𝐺
1|⋅ℎ, whe e ℎis he cha ac e is ic size o
he mesh and 𝐿is he e e ence leng h. The bounda y coo dina es a e ob ained o 𝜒0=, wi h he dis ance be ween he body
su ace and he poin on he shi ed bounda y x𝑠𝑏.
The dono s encil mus be a enough om he su ace o a oid in e dependencies. A e some es s, we figu ed ou expe imen ally
ha a good ade off o he defini ion o he e e ence leng h 𝐿seems o be:
𝐿=+𝛼√𝑑⋅ℎ, (6)
whe e we ecall ha 𝑑is he dimension o he p oblem and 𝛼is a use -defined pa ame e . Finally, when selec ing he alue o he
pa ame e 𝛼, diffe en choices a e possible; howe e , we ound ha a solid choice is o se 𝛼=1+
ℎ
( o mo e de ails, see Appendix D).
Ha ing fixed one poin on he body su ace, we a e le wi h 𝑃addi ional poin s o pe o m he in e pola ion a he shi ed
bounda y gi en by Eq. (5). The (𝑃+1)
𝑡ℎ-poin co esponds o he bounda y poin x𝑏=x𝑠𝑏 +
n𝑠whe e he solu ion is known.
Once he s encil has been compu ed, we pe o m a mapping om he physical space along he no mal
n𝑠(defined by he scala 𝜒
in Eq. (5)) o a local e e ence ame: Ψ(𝜒)→𝜉, wi h 𝜉∈[−1,1]. A simple linea map sa is ying Ψ(0) = −1, Ψ(𝐿+𝑑𝐿)=1is used:
Ψ(𝜒)=2 𝜒
𝐿+𝑑𝐿 −1.(7)
F om he map, one can ge he local e e ence ame coo dina es o each poin :
𝜉𝑗=2 𝜒𝑗
(𝐿+𝑑𝐿)−1
wi h 𝑗=0,…,𝑃. (8)
Fig. 6shows he mapping including he o iginal s encil in he physical space and he one ob ained a e he applica ion o he map.
The Gauss nodes 𝜉𝐺
𝑗 o an app oxima ion o 𝑃=3a e shown oge he wi h all he quan i ies defined in Eq. (5). Once all geome ic
quan i ies a e compu ed, he s a e is e alua ed a each in e pola ion poin by simply e alua ing he polynomial ep esen a ion o he
solu ion a hose loca ions. Recall ha he s a e o he fi s poin 𝜉0, coming om he in e sec ion be ween he body no mal
n𝑠and
he body su ace, is known and comes om he bounda y condi ions. A single elemen , 𝑒𝑖𝑛𝑡, can be used o compu e he s a e a each
in e pola ion poin e en i he s encil does no belong o ha specific elemen ; Fig. 7shows he s encil and he elemen 𝑒𝑖𝑛𝑡 used o
econs uc he da a.
In he local coo dina e 𝜉, he in e pola ion cons uc ed can be w i en as
u(𝜉)=
𝑃
∑
𝑖=0
𝑙𝑖(𝜉)u𝑖,(9)
𝑙𝑖(𝜉)being he in e pola ion unc ions. The known alues o his in e pola ion a e:
u(𝜉0)=u𝑏, u(𝜉𝑖)=u𝑖, 𝑖=1,…,𝑃,
whe eas he bounda y condi ions o be p esc ibed on he su oga e bounda y a e:
u𝑠𝑏 =u(−1) =
𝑃
∑
𝑖=0
𝑙𝑖(−1)u𝑖.
Exp ession (9)is he key o he p oposed me hod. I allows one o ob ain u𝑠𝑏 in e ms o u𝑏(o o he de i a i es o u, see below) and
o he alues o he dono s encil u𝑖, 𝑖=1,…,𝑃. Since hese a e unknown, hey will con ibu e o he s iffness ma ix o he p oblem,
hus making i no symme ic; indeed o each DoF belonging o he shi ed aces, Eq. (9) will con ibu e o he Jacobian, and once
he sys em o equa ions is sol ed, he alues on he shi ed bounda ies a e upda ed. Howe e , his is i ele an i he equa ions a e
sol ed using an explici ime s epping, as men ioned ea lie . The effec is simila o ha o he SBM, in which he Taylo expansion
used o ob ain u𝑠𝑏 in ol es he de i a i es o u, which need o be exp essed in e ms o he nodal alues o ui sel and con ibu e o
he s iffness ma ix, when i is assembled.
An impo an ema k is ha ins ead o knowing u𝑏we may know i s no mal de i a i e, i.e., he de i a i e in he di ec ion o
𝜉. This also allows us o compu e u𝑏using Eq. (9)and, om his, ob ain again u𝑠𝑏. This also opens he doo o p esc ibe Neumann
bounda y condi ions. In he ollowing subsec ion, we explain how o apply his p ocedu e o Poisson’s p oblem.
We nex de ail he implemen a ion o bounda y condi ions o he Poisson equa ion and o he Na ie -S okes equa ions. Ou
app oach enables us o assign a p ede e mined alue o he s a e ec o , u, on he body su ace, 𝜕Ω, exp essed as u|𝜕Ω=u𝑏.
3.3. Di ichle bounda y condi ion o he Poisson equa ion
Le us conside a gene ic Poisson equa ion:
−∇2
x𝑢=𝑆, x∈Ω,(10)
Jou nal o Compu a ional Physics 528 (2025) 113807
7
S. Colombo, G. Rubio, J. Kou e al.
wi h Di ichle bounda y condi ions:
𝑢=𝑢𝑏on 𝜕Ω,(11)
whe e 𝑢is now a scala unc ion and 𝑆is a sou ce e m.
We can econs uc 𝑢along he no mal di ec ion using he one-dimensional Lag angian polynomial ou lined in he p eceding
sec ion, gi en in Eq. (9).
No e ha he in e pola ion poin s a e he poin s inside he fluid, i.e., 𝑖=1,…,𝑃, and he poin on he body co esponds o 𝑖=0.
Since he fi s in e pola ion poin (𝑖=0) lies on he body, i s s a e is known:
𝑢0=𝑢𝑏,(12)
while 𝑢𝑖, wi h 𝑖>0, a e he nodal alues a he poin s o he dono cell desc ibed in he p e ious sec ion. The alue on he shi ed
bounda y can now be ound as
𝑢𝑠𝑏 =𝑢|𝜉=−1 =𝑙0(−1)𝑢𝑏
⏟ ⏞⏟ ⏞⏟
Di ichle B.C.
+
𝑃
∑
𝑖=1
𝑙𝑖(−1)𝑢𝑖
⏟ ⏞⏞⏞⏟ ⏞⏞⏞⏟
Fluid e alua ion
.(13)
3.4. Neumann bounda y condi ion o he Poisson equa ions
The Neumann bounda y condi ion can be w i en as:
𝜕𝑢
𝜕
n𝑠
=𝐺, (14)
wi h 𝐺gi en. Using he chain ules o he de i a i e, he le -hand side o Eq. (14)we ge :
𝜕𝑢
𝜕
n𝑠
=𝜕𝑢
𝜕𝜉
𝜕𝜉
𝜕
n𝑠
;(15)
in pa icula , keeping in mind ha he econs uc ion akes place along
n𝑠 h ough a one-dimensional s encil and conside ing Eq. (8)
one ob ains:
𝜕𝑢
𝜕𝜉
𝜕𝜉
𝜕
n𝑠
=𝜕𝑢
𝜕𝜉
d𝜉
d𝜒=𝜕𝑢
𝜕𝜉
2
(𝐿+𝑑𝐿).(16)
Finally, combining Eq. (14)and Eq. (16)and disc e izing he solu ion using he one-dimensional Lag angian polynomial as i was
done o he Di ichle bounda y condi ions, we can w i e:
𝑃
∑
𝑖=0
𝑙′
𝑖(𝜉)𝑢𝑖=(𝐿+𝑑𝐿)
2
𝐺, (17)
whe e 𝑙′
𝑖(𝜉)is he de i a i e o he Lag angian polynomial. I is now possible o compu e he alue 𝑢0so ha Eq. (17)is sa isfied,
ob aining:
𝑢0=
(𝐿+𝑑𝐿)
2
𝐺−∑𝑃
𝑖=1 𝑙′
𝑖(𝜉0)𝑢𝑖
𝑙′
0(𝜉0)
.(18)
Once he s a e a 𝑢0is known, we p oceed as compu ing he alue on he shi ed bounda y:
𝑢𝑠𝑏 =𝑢|𝜉=−1 =𝑙0(−1)𝑢0
⏟ ⏞⏟ ⏞⏟
Neumann B.C.
+
𝑃
∑
𝑖=1
𝑙𝑖(−1)𝑢𝑖
⏟ ⏞⏞⏞⏟ ⏞⏞⏞⏟
Fluid e alua ion
.(19)
No e ha he econs uc ion o he unknown can be based on ei he knowing i o i s no mal de i a i e on he physical bounda y,
bu in bo h cases 𝑢𝑠𝑏 is weakly p esc ibed on he su oga e bounda y.
3.5. Su ace ep esen a ion
The disc e e ep esen a ion o he body su ace has a c ucial ole in he amewo k o high-o de me hods. This is also ue o ou
me hodology, whe e he posi ion o he body su ace and i s no mal
n𝑠, a e equi ed. When an analy ical desc ip ion o he shape
is a ailable, he exac no mal body is
n𝑠. Mo e gene ally, in eal applica ion he analy ical o mula is no accessible and hus he
geome y and no mals a e p o ided by a CAD so wa e. In his wo k, we use he STL o ma o he geome y, and we calcula e
no mals based on his geome y inside ou sol e . An STL is a file made o uns uc u ed linea iangula elemen s (Fig. 8) and can
be p o ided by mos o he CAD so wa e a ailable oday. This choice is due o he ac ha STL a e widely suppo ed by mos o he
Jou nal o Compu a ional Physics 528 (2025) 113807
8
S. Colombo, G. Rubio, J. Kou e al.
Fig. 8. Examples o STL files.
so wa e, a e ligh in e ms o size, and a e open sou ce. The d awback o STLs is ha he iangula ion is made o linea elemen s,
and hence he cu a u e migh no be accu a e enough, bu in hese cases one can inc ease he numbe o iangles o he STL.
4. Tes cases
The aim o his sec ion is o e i y he accu acy o he me hod. We will check he accu acy o he me hod o a ange o p oblems:
a 1D ad ec ion-diffusion equa ion, a Poisson p oblem including cu elemen s and cu ed elemen s, and finally a ci cula cylinde
d i en by he comp essible Na ie -S okes equa ions.
The fi s es is a 1D linea con ec ion-diffusion equa ion, which shows he capabili y o e ie ing high-o de accu acy e en when
aking in o accoun he bounda ies. Addi ionally, i is compa ed wi h he olume penaliza ion IBM o show he supe io i y o he new
me hod. The second case is a ully 2D Poisson equa ion sol ed on a squa ed bounda y whe e no cu a u e is p esen . In pa icula , we
show wha happens when he shi ed aces coincide wi h he physical bounda y and when hey do no leading o cu -cells. To check
bounda y condi ions on cu ed bounda ies, a hi d es case is p oposed in which he s eady-s a e hea equa ion is sol ed on wo
concen ic ci cles. The las case compa es ou me hod, he body-fi ed and olume penaliza ion o he comp essible Na ie -S okes.
In all es cases, he Runge-Ku a (RK3) ma ching scheme is used o ad ance in ime, e en when dealing wi h s eady p oblems. In he
fi s es cases whe e he con e gence is s udied, he geome y is analy ical and he body no mals a e analy ically compu ed; on he
o he hand, in he Na ie -S okes case, a STL file is used.
4.1. 1D linea con ec ion-diffusion
We fi s es ou me hodology using a 1D linea ad ec ion-diffusion equa ion, inspi ed om [22]:
𝜕𝑢
𝜕𝑡
+𝑐𝜕𝑢
𝜕𝑥 −𝜈𝜕2𝑢
𝜕𝑥2=0 𝑥∈[0,1],(20)
whe e 𝑐=0.1is he ad ec ion speed and 𝜈=0.01 is he diffusi i y. Homogeneous bounda y condi ions a e applied: 𝑢(0,𝑡)=0and
𝑢(1,𝑡)=0. The ime s ep used o he simula ions is se o Δ𝑡=1.0⋅10−7 o ensu e ime accu acy. The analy ical solu ion o Eq. (20)
is
𝑢(𝑥,𝑡)=e
5𝑥−𝑡(𝜈𝜋2+25𝜈)sin(𝜋𝑥).(21)
F om Eq. (21)one can easily ob ain he ini ial condi ion by se ing 𝑡=0.
To make his es sui able o he IBM, wo addi ional elemen s a e added a he ex eme o he compu a ional domain. In he
classical olume penaliza ion, a sou ce e m is added o Eq. (20) o simula e he p esence o he bounda ies. The sou ce e m is
applied o all he deg ees o eedom lying ou side he compu a ional domain, i.e., all he poin s such ha 𝑥<0and 𝑥>1, see Fig. 9.
Eq. (20)becomes:
𝜕𝑢
𝜕𝑡
+𝑐𝜕𝑢
𝜕𝑥 −𝜈𝜕2𝑢
𝜕𝑥2=−1
𝜂(𝑢−𝑢𝑠),(22)
whe e 𝜂is he penaliza ion e m applied ou side (0,1) ha , in his simula ion, is se o 𝜂=Δ𝑡( he ime s ep size) and 𝑢𝑠is he
solu ion we wan o impose on he bounda ies, in his case 𝑢𝑠=0. The final simula ion ime is 𝑡𝑚𝑎𝑥 =0.01. Fig. 10 shows he esul s
ob ained o adi ional olume penaliza ion (Fig. 10a) and hose coming om he high-o de app oach (Fig. 10b), when disc e izing
he domain wi h 𝑁𝑒𝑙 elemen s. As expec ed [22], he olume penaliza ion me hod loses he high o de accu acy nea bounda ies,
while o ou new p oposed me hod he high o de a e o con e gence (𝑁−(𝑃+1)
𝑒𝑙 )is main ained o all polynomial o de s.
Jou nal o Compu a ional Physics 528 (2025) 113807
9
S. Colombo, G. Rubio, J. Kou e al.
Fig. 21. Sphe e a Re = 100: P essu e coefficien compa ison be ween high-o de imme sed bounda y and olume penaliza ion.
Table 2
Sphe e a Re = 100: Compa ison be ween
he d ag coefficien 𝑐𝐷.
Sou ce 𝑐𝐷
Ho ses3D Body-fi ed 1.087
Ho ses3D IBM VP 1.045
Ho ses3D High-o de IBM 1.081
Fadlun e al. [13] 1.0794
Fo nbe g [49] 1.0852
Finally, Table 2summa izes he d ag coefficien s o he olume penaliza ion and he new me hodology oge he wi h published
da a. The d ag p edic ed wi h he new me hodology is close o he e e ence alues, han he one ob ained wi h he olume penal-
iza ion app oach.
5. Conclusions
We ha e p oposed a IBM capable o p ese ing high-o de accu acy, pa icula ly in he ea men o bounda y condi ions. Ou
app oach main ains he o mal o de o he scheme (𝑃+1) when using polynomials o o de 𝑃, a oiding he deg ada ion nea
geome ies obse ed in o he unfi ed bounda y me hods (e.g., olume penaliza ion). By using Lag angian polynomials along he
body no mal o econs uc he s a e on specific aces, he me hod efficien ly exploi s Ca esian con o ming meshes o simula e
cu ed geome ies. Implemen ed wi hin a discon inuous Gale kin amewo k, using he open sou ce Ho ses3D sol e , his s a egy
has been es ed on he con ec ion-diffusion equa ion, he Poisson equa ion, and he comp essible Na ie -S okes equa ions. These es s
confi med he me hod’s abili y o achie e high-o de accu acy nea he body, e en in he p esence o cu elemen s including cu ed
su aces. Fu he mo e, ou s udy demons a ed ha he p oposed me hod ou pe o ms olume penaliza ion o he comp essible
Na ie -S okes equa ions, as e idenced by he imp o ed p essu e dis ibu ion on he body and mo e accu a e d ag coefficien s. In
conclusion, his wo k p esen s a uly high-o de IBM. Fu u e esea ch will ocus on applying his me hodology o mo e complex flow
egimes, including u bulence and mo ing geome ies, o u he alida e and ex end i s applicabili y.
CRediT au ho ship con ibu ion s a emen
S. Colombo: W i ing – e iew & edi ing, W i ing – o iginal d a , Visualiza ion, Valida ion, So wa e, Me hodology, Fo mal
analysis, Da a cu a ion, Concep ualiza ion. G. Rubio: W i ing – e iew & edi ing, W i ing – o iginal d a , Supe ision, P ojec
adminis a ion, Me hodology, Funding acquisi ion, Fo mal analysis. J. Kou: W i ing – e iew & edi ing, W i ing – o iginal d a ,
Supe ision. E. Vale o: P ojec adminis a ion, Funding acquisi ion. R. Codina: W i ing – e iew & edi ing, W i ing – o iginal d a ,
Supe ision, Me hodology, In es iga ion, Fo mal analysis. E. Fe e : W i ing – e iew & edi ing, W i ing – o iginal d a , Supe ision,
Resou ces, P ojec adminis a ion, Me hodology, In es iga ion, Funding acquisi ion, Fo mal analysis.
Decla a ion o compe ing in e es
The au ho s decla e he ollowing financial in e es s/pe sonal ela ionships which may be conside ed as po en ial compe ing
in e es s: S e ano Colombo epo s financial suppo was p o ided by Eu opean Union. I he e a e o he au ho s, hey 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 .
Jou nal o Compu a ional Physics 528 (2025) 113807
16
S. Colombo, G. Rubio, J. Kou e al.
Acknowledgemen s
SC acknowledges he unding om he Eu opean Union’s Ho izon 2020 esea ch and inno a ion p og amme unde he Ma ie
Skłodowska Cu ie g an ag eemen No 955923-SSECOID. RC g a e ully acknowledges he suppo ecei ed om he ICREA Acadèmia
P og am, om he Ca alan Go e nmen . EF would like o hank he suppo o Agencia Es a al de In es igación ( o he g an “Eu opa
Excelencia 2022” P oyec o EUR2022-134041/AEI/10.13039/501100011033) y del Mecanismo de Recupe ación y Resiliencia de la
Unión Eu opea. EF and GR acknowledge he unding ecei ed by he G an DeepCFD (P ojec No. PID2022-137899OB-I00) unded by
MCIN/ AEI/10.13039/501100011033 and by ERDF A way o making Eu ope. This esea ch has ecei ed unding om he Eu opean
Union (ROSAS, p ojec numbe 101138319). This esea ch has ecei ed unding om he Eu opean Union (ERC, Off-cous ics, p ojec
numbe 101086075). Views and opinions exp essed a e, howe e , hose o he au ho s only and do no necessa ily eflec hose
o he Eu opean Union o he Eu opean Resea ch Council. Nei he he Eu opean Union no he g an ing au ho i y can be held
esponsible o hem. Finally, all au ho s g a e ully acknowledge he Uni e sidad Poli écnica de Mad id (www.upm.es) o p o iding
compu ing esou ces on Mage i Supe compu e and he compu e esou ces a Ma eNos um and he echnical suppo po p o ided
by Ba celona Supe compu ing Cen e (RES-IM-2024-1-0003).
Appendix A. Volume penaliza ion imme sed bounda y
A b ie explana ion o he olume penaliza ion imme sed bounda y is he e p oposed o cla i y.
The olume penaliza ion IBM consis s in imposing a sou ce e m o he Na ie -S okes equa ions o simula e he p esence o a
body. As is common in his me hod, a mask 𝜁(x,𝑡)is in oduced in o de o dis inguish he egion o he domain Ωinside he body
Ω𝑏 om he one ou side Ω𝑓, wi h Ω=Ω
𝑏∪Ω
𝑓:
𝜁(x,𝑡)={1, i x∈Ω
𝑏
0 i x∈Ω
𝑓.(A.1)
The mask is compu ed using a ay- acing app oach ha , acco ding o he numbe o in e sec ions, allows one o check i a deg ee o
eedom is inside he body (odd numbe o in e sec ions) o ou side (e en numbe o in e sec ions).
The Na ie -S okes equa ions, wi h he addi ion o he sou ce e m, become:
u𝑡+∇
x⋅(
F𝑒−
F𝑣)=
S, x∈Ω, 𝑡>0.(A.2)
Fo he olume penaliza ion, he sou ce e m o Di ichle bounda y condi ions is:
S=𝜁
𝜂⎛⎜⎜⎜⎜⎜⎝
0
𝜌𝑢𝑠−𝜌𝑢
𝜌𝑣𝑠−𝜌𝑣
𝜌𝑤𝑠−𝜌𝑤
𝜌
2
(𝑢2
𝑠+𝑣2
𝑠+𝑤2
𝑠)− 𝜌
2
(𝑢2+𝑣2+𝑤2)
⎞⎟⎟⎟⎟⎟⎠
,(A.3)
whe e 𝑠=(𝑢𝑠,𝑣𝑠,𝑤𝑠)𝑇is he eloci y o be imposed on he body. Fo no-slip bounda y condi ions, 𝑠=(0,0,0)𝑇; 𝜂is he penaliza ion
pa ame e . The penaliza ion pa ame e should be high enough o p ope ly simula e he po osi y o he body. Howe e , his inc eases
he s iffness o he p oblem and educes he ime s ep o he simula ion. The common p ac ice in he imme sed bounda y communi y
is o se he penaliza ion pa ame e equal o he explici ime s ep o he scheme, i.e., 𝜂=Δ𝑡. Mo e de ails can be ound in [23], [22],
[50], [51], among o he s.
Appendix B. Bassi-Rebay 1
As i is common p ac ice o he Bassi-Rebay 1 o mula ion (BR1) [52], an addi ional a iable is added o he Na ie -S okes
equa ions, g=∇
xu, leading o a coupled sys em o equa ions:
𝐽g−∇
𝝃u=0(B.1)
𝐽u𝑡+∇
𝝃⋅(F𝑒−F𝑣)−𝐽S=0.(B.2)
Conside ing he weak o mula ion o Eq. (B.1), a jump in he s a e uac oss he in e ace be ween wo elemen s appea s:
∑
𝑒𝑙 ∫
𝑒𝑙
[𝜙𝑖𝐽g+∇
𝝃𝜙⋅u]−∑
𝜕𝑒𝑙 ∮
𝜕𝑒𝑙
𝜙𝑖u⋅
n=0.(B.3)
No e ha , when ans o med back o he physical domain, a ac o 1∕ℎappea s in he las e m. The jump appea ing in Eq. (B.3)is
defined, o he BR1 case, as:
u=1
2(u+−u−).(B.4)
Jou nal o Compu a ional Physics 528 (2025) 113807
17
S. Colombo, G. Rubio, J. Kou e al.
Fig. 22. Fi s case o in e es in he e alua ion o he shi ed aces; he shi ed aces a e depic ed in blue.
Fig. 23. Second case o in e es in he e alua ion o he shi ed aces; he shi ed aces a e depic ed in blue.
When a shi ed ace is conside ed, he jump appea ing in Eq. (B.4)is modified and eplaced by he compu ed shi ed s a e u𝑠𝑏:
u𝑠𝑏 =1
2(u𝑠𝑏 −u−).(B.5)
Consequen ly, he g adien g=1
𝐽∇𝝃uis e alua ed aking in o accoun he p esence o he body and is used o compu e he iscous
nume ical fluxes F⋆
𝑣(u+,u−,g+,g−,
n); in pa icula , when conside ing a shi ed ace, we se g𝑠𝑏 =g−, ob aining F⋆
𝑣(u𝑠𝑏,u−,g𝑠𝑏,g−,
n).
No e ha simila app oaches can be ollowed o o he iscous fluxes, such as in e io penal y o Bassi-Rebay 2.
Appendix C. Some conside a ion abou he iden ifica ion o he shi ed aces
In o de o cla i y some aspec s abou he iden ifica ion o he shi ed aces, some in e es ing cases a e he e epo ed and analyzed.
In pa icula , we p o ide some examples o he wo possible condi ions ha can occu du ing he p ocess o selec ion o he shi ed
aces, i.e. 1 -when a ace is agged wice as a shi ed ace and 2 -when one elemen has one ace wi h a leas one deg ee o eedom
inside he fluid domain (see sec ion 3.1).
An example o he fi s si ua ion is epo ed in Fig. 22. We ocus on he le ace o elemen 𝐵highligh ed in blue in Fig. 22a. Following
he p ocedu e epo ed in sec ion 3.1, his ace is agged wice as a shi ed ace and, he e o e i is disca ded, leading o he final
shi ed bounda y Γ𝑠is shown in Fig. 22b.
The o he possible condi ion is epo ed in Fig. 23. He e, he igh and le aces o elemen 𝐵a e agged as shi ed aces (Fig. 23a)
leading o a shi ed bounda y ha is no closed. I his occu s, hese wo aces a e dele ed om Γ𝑠as hey a e no in on o he ace
ha has one o mo e deg ees o eedom in he fluid domain. The final bounda y Γ𝑠is epo ed in Fig. 23b.
I is impo an o no ice ha hese condi ions a ise when complex geome ies a e disc e ized using coa se meshes. Refining he mesh
locally helps o alle ia e he occu ence o hese p oblems and also allows o be e cap u e he shape o he body.
Appendix D. Effec o he a ia ion o he pa ame e 𝜶
This sec ion is dedica ed o a desc ip ion o he beha io o he p oposed me hod as a unc ion o he pa ame e 𝛼. In pa icula ,
we use as a e e ence he s eady hea equa ion, which is desc ibed in de ail in Sec ion 4.4, and pe o m simula ions o a ious
pa ame e s 𝛼, on a mesh wi h 25 × 25 elemen s. The esul s a e shown in Fig. 24 whe e 𝛼=𝜖(1+
ℎ
)√𝑑⋅ℎ.
Fig. 24 shows ha he alue o 𝛼sligh ly affec s he e o . When 𝛼dec eases, he e o educes o he same polynomial o de . No e
ha he e is a limi a ion on he minimum achie able alue o he pa ame e . In Fig. 24 he black diamonds ep esen he minimum
alue o 𝛼(𝛼𝑚) below which he scheme does no con e ge: no ably, he highe he polynomial o de , he g ea e he minimum
allowable alue o 𝛼. Fu he mo e, as we educe he alue o 𝛼, he posi ion o he s encil poin s ge s close o he su ace. The use
Jou nal o Compu a ional Physics 528 (2025) 113807
18
S. Colombo, G. Rubio, J. Kou e al.
Fig. 24. Beha io o he e o as 𝛼changes.
o a 𝛼𝑚is jus ified by he ac ha a ecu si e dependency is expec ed as we mo e close o he body, and his ac o helps o eco e
s abili y and con e gence. Rega ding 𝑃=1, no alue o 𝛼𝑚is ound, which sugges s ha he e m 𝑑𝐿 (see Sec ion 3.2) is always la ge
enough o a oid a ecu si e dependency. In ac , we can check ha educing he alue o 𝑑𝐿 leads o he same beha io as epo ed
o highe polynomial o de s. In his wo k, he alue o 𝛼co esponds o 𝜖=1, which is a good comp omise in e ms o e o and
sa e y ma gin om 𝛼𝑚, since he same can be used o all polynomials.
Da a a ailabili y
Da a will be made a ailable on eques .
Re e ences
[1] N. Chalme s, G. Agbaglah, M. Ch us , C. Ma iplis, A pa allel ℎ𝑝-adap i e high o de discon inuous Gale kin me hod o he incomp essible Na ie -S okes
equa ions, J. Compu . Phys. X 2 (2019) 100023, h ps://doi.o g/10.1016/j.jcpx.2019.100023.
[2] Z. Wang, K. Fidkowski, R. Abg all, F. Bassi, D. Ca aeni, A. Ca y, H. Deconinck, R. Ha mann, K. Hillewae , H. Huynh, N. K oll, G. May, P.-O. Pe sson, B. an Lee ,
M. Visbal, High-o de CFD me hods: cu en s a us and pe spec i e, In . J. Nume . Me hods Fluids 72 (8) (2013) 811–845, h ps://doi.o g/10.1002/fld.3767,
h ps://onlinelib a y.wiley.com/doi/abs/10.1002/fld.3767.
[3] G. Ka niadakis, S. She win, Spec al/hp Elemen Me hods o Compu a ional Fluid Dynamics, Ox o d Uni e si y P ess, 2005, h ps://doi.o g/10.1093/acp o :
oso/9780198528692.001.0001.
[4] G. Apa icio-Es ems, A. Ga gallo-Pei ó, X. Roca, Defining me ic-awa e size-shape measu es o alida e and op imize cu ed high-o de meshes, Compu . Aided
Des. 168 (2024) 103667, h ps://doi.o g/10.1016/j.cad.2023.103667.
[5] G. Apa icio-Es ems, A. Ga gallo-Pei ó, X. Roca, Combining high-o de me ic in e pola ion and geome y implici iza ion o cu ed -adap ion, Compu . Aided
Des. 157 (2023) 103478, h ps://doi.o g/10.1016/j.cad.2023.103478.
[6] C.S. Peskin, Flow pa e ns a ound hea al es: a nume ical me hod, J. Compu . Phys. 10 (2) (1972) 252–271, h ps://doi.o g/10.1016/0021-9991(72)90065-4.
[7] R. Mi al, G. Iacca ino, Imme sed bounda y me hods, Annu. Re . Fluid Mech. 37 (2005) 239–261.
[8] F. So i opoulos, X. Yang, Imme sed bounda y me hods o simula ing fluid-s uc u e in e ac ion, P og. Ae osp. Sci. 65 (2014) 1–21.
[9] B.E. G iffi h, N.A. Pa anka , Imme sed me hods o fluid–s uc u e in e ac ion, Annu. Re . Fluid Mech. 52 (2020) 421–448.
[10] T. Ye, R. Mi al, H. Udaykuma , W. Shyy, An accu a e Ca esian g id me hod o iscous incomp essible flows wi h complex imme sed bounda ies, J. Compu .
Phys. 156 (2) (1999) 209–240.
[11] H. Udaykuma , R. Mi al, P. Rampunggoon, A. Khanna, A sha p in e ace Ca esian g id me hod o simula ing flows wi h complex mo ing bounda ies, J. Compu .
Phys. 174 (1) (2001) 345–380.
[12] P. Fu, G. K eiss, S. Zahedi, A bound p ese ing cu discon inuous Gale kin me hod o one dimensional hype bolic conse a ion laws, a Xi :2404.13936, h ps://
a xi .o g/abs/2404.13936, 2024.
[13] E. Fadlun, R. Ve zicco, P. O landi, J. Mohd-Yuso , Combined imme sed-bounda y fini e-diffe ence me hods o h ee-dimensional complex flow simula ions,
J. Compu . Phys. 161 (1) (2000) 35–60.
[14] N. Zhang, Z. Zheng, An imp o ed di ec - o cing imme sed-bounda y me hod o fini e diffe ence applica ions, J. Compu . Phys. 221 (1) (2007) 250–268, h ps://
doi.o g/10.1016/j.jcp.2006.06.012.
[15] H. Luo, H. Dai, P.J.F. de Sousa, B. Yin, On he nume ical oscilla ion o he di ec - o cing imme sed-bounda y me hod o mo ing bounda ies, Compu . Fluids 56
(2012) 61–76.
[16] F.-B. Tian, H. Dai, H. Luo, J.F. Doyle, B. Rousseau, Fluid–s uc u e in e ac ion in ol ing la ge de o ma ions: 3d simula ions and applica ions o biological sys ems,
J. Compu . Phys. 258 (2014) 451–469.
[17] S. Majumda , G. Iacca ino, P. Du bin, RANS sol e s wi h adap i e s uc u ed bounda y non-con o ming g ids, Annu. Res. B . 1 (2001).
[18] Y.-H. Tseng, J.H. Fe zige , A ghos -cell imme sed bounda y me hod o flow in complex geome y, J. Compu . Phys. 192 (2) (2003) 593–623, h ps://doi.o g/
10.1016/j.jcp.2003.07.024.
[19] P. Ango , C.-H. B uneau, P. Fab ie, A penaliza ion me hod o ake in o accoun obs acles in incomp essible iscous flows, Nume . Ma h. 81 (4) (1999) 497–520.
[20] E. B own-Dymkoski, N. Kasimo , O.V. Vasilye , A cha ac e is ic based olume penaliza ion me hod o gene al e olu ion p oblems applied o comp essible
iscous flows, J. Compu . Phys. 262 (2014) 344–357, h ps://doi.o g/10.1016/j.jcp.2013.12.060.
[21] R. Abg all, H. Beaugend e, C. Dob zynski, An imme sed bounda y me hod using uns uc u ed aniso opic mesh adap a ion combined wi h le el-se s and penal-
iza ion echniques, J. Compu . Phys. 257 (2014) 83–101.
[22] J. Kou, S. Joshi, A.H. de Mendoza, K. Pu i, C. Hi sch, E. Fe e , Imme sed bounda y me hod o high-o de flux econs uc ion based on olume penaliza ion,
J. Compu . Phys. 448 (2022) 110721, h ps://doi.o g/10.1016/j.jcp.2021.110721.
Jou nal o Compu a ional Physics 528 (2025) 113807
19
S. Colombo, G. Rubio, J. Kou e al.
[23] J. Kou, E. Fe e , A combined olume penaliza ion/selec i e equency damping app oach o imme sed bounda y me hods applied o high-o de schemes,
J. Compu . Phys. 472 (2023) 111678.
[24] J. Kou, A.H. de Mendoza, S. Joshi, S. Le Clainche, E. Fe e , Eigensolu ion analysis o imme sed bounda y me hod based on olume penaliza ion: applica ions o
high-o de schemes, J. Compu . Phys. 449 (2022) 110817, h ps://doi.o g/10.1016/j.jcp.2021.110817.
[25] V.J. Llo en e, J. Kou, E. Vale o, E. Fe e , A modified equa ion analysis o imme sed bounda y me hods based on olume penaliza ion: applica ions o linea
ad ec ion–diffusion equa ions and high-o de discon inuous Gale kin schemes, Compu . Fluids 257 (2023) 105869, h ps://doi.o g/10.1016/j.compfluid.2023.
105869.
[26] J. Kou, E. Fe e , A combined olume penaliza ion/selec i e equency damping app oach o imme sed bounda y me hods: applica ion o mo ing geome ies,
Phys. Fluids 35 (12) (2023) 121702, h ps://doi.o g/10.1063/5.0179779.
[27] T. Engels, D. Kolomenskiy, K. Schneide , J. Ses e henn, Nume ical simula ion o fluid–s uc u e in e ac ion wi h he olume penaliza ion me hod, J. Compu .
Phys. 281 (2015) 96–115, h ps://doi.o g/10.1016/j.jcp.2014.10.005.
[28] A. Massing, M. La son, A. Logg, M. Rognes, A Ni sche-based cu fini e elemen me hod o a fluid–s uc u e in e ac ion p oblem, Commun. Appl. Ma h. Compu .
Sci. 10 (11 2013), h ps://doi.o g/10.2140/camcos.2015.10.97.
[29] B. Liu, A Ni sche s abilized fini e elemen me hod: applica ion o hea and mass ans e and fluid–s uc u e in e ac ion, Compu . Me hods Appl. Mech. Eng.
386 (2021) 114101, h ps://doi.o g/10.1016/j.cma.2021.114101.
[30] E. Bu man, P. Hansbo, Fic i ious domain fini e elemen me hods using cu elemen s: II. a s abilized Ni sche me hod, Appl. Nume . Ma h. 62 (2012) 328–341.
[31] J. Kim, D. Kim, H. Choi, An imme sed-bounda y fini e- olume me hod o simula ions o flow in complex geome ies, J. Compu . Phys. 171 (1) (2001) 132–150,
h ps://doi.o g/10.1006/jcph.2001.6778.
[32] M. Uhlmann, An imme sed bounda y me hod wi h di ec o cing o he simula ion o pa icula e flows, J. Compu . Phys. 209 (2) (2005) 448–476, h ps://
doi.o g/10.1016/j.jcp.2005.03.017.
[33] F. So i opoulos, X. Yang, Imme sed bounda y me hods o simula ing fluid–s uc u e in e ac ion, P og. Ae osp. Sci. 65 (01 2013), h ps://doi.o g/10.1016/j.
pae osci.2013.09.003.
[34] A. Main, G. Sco azzi, The shi ed bounda y me hod o embedded domain compu a ions. Pa I: Poisson and S okes p oblems, J. Compu . Phys. 372 (2018)
972–995, h ps://doi.o g/10.1016/j.jcp.2017.10.026.
[35] N.M. A allah, C. Canu o, G. Sco azzi, Analysis o he shi ed bounda y me hod o he Poisson p oblem in gene al domains, a Xi :2006.00872, 2020.
[36] M. Funada, T. Imamu a, High-o de imme sed bounda y me hod o in iscid flows applied o flux econs uc ion me hod on a hie a chical Ca esian g id, Compu .
Fluids 265 (2023) 105986, h ps://doi.o g/10.1016/j.compfluid.2023.105986.
[37] S. Pé on, C. Benoi , T. Renaud, I. Ma y, An imme sed bounda y me hod on Ca esian adap i e g ids o he simula ion o comp essible flows a ound a bi a y
geome ies, Eng. Compu . 37 (07 2021), h ps://doi.o g/10.1007/s00366-020-00950-y.
[38] B. Cons an , S. Pé on, H. Beaugend e, C. Benoi , An imp o ed imme sed bounda y me hod o u bulen flow simula ions on Ca esian g ids, J. Compu . Phys.
435 (2021) 110240, h ps://doi.o g/10.1016/j.jcp.2021.110240.
[39] R. Codina, J. Baiges, App oxima e imposi ion o bounda y condi ions in imme sed bounda y me hods, In . J. Nume . Me hods Eng. 80 (2009) 1379–1405.
[40] E. Fe e , G. Rubio, G. N oukas, W. Laskowski, O. Ma iño, S. Colombo, A. Ma eo-Gabín, H. Ma bona, F.M. de La a, D. Hue go, J. Manzane o, A. Rueda-Ramí ez,
D. Kop i a, E. Vale o, A high-o de discon inuous Gale kin sol e o flow simula ions and mul i-physics applica ions, Compu . Phys. Commun. 287 (2023)
108700, h ps://doi.o g/10.1016/j.cpc.2023.108700.
[41] D.A. Kop i a, Implemen ing Spec al Me hods o Pa ial Diffe en ial Equa ions: Algo i hms o Scien is s and Enginee s, 1s edi ion, Sp inge Publishing Company,
Inco po a ed, 2009.
[42] G.J. Gassne , A.R. Win e s, F.J. Hindenlang, D.A. Kop i a, The b 1 scheme is s able o he comp essible Na ie –S okes equa ions, J. Sci. Compu . 77 (2018)
154–200.
[43] B. Ri ie e, Discon inuous Gale kin Me hods o Sol ing Ellip ic and Pa abolic Equa ions: Theo y and Implemen a ion, ol. 35, 2008.
[44] B. Cockbu n, G.E. Ka niadakis, C.-W. Shu, Discon inuous Gale kin Me hods: Theo y, Compu a ion and Applica ions, 1s edi ion, Sp inge Publishing Company,
Inco po a ed, 2011.
[45] T. Akenine-Mölle , E. Haines, N. Hoffman, A. Pesce, M. Iwanicki, S. Hillai e, Real-Time Rende ing, 4 h edi ion, A K Pe e s/CRC P ess, Boca Ra on, FL, USA, 2018.
[46] S. Dennis, G.-z. Chang, Nume ical solu ions o s eady flow pas a ci cula cylinde a Reynolds numbe s up o 100, J. Fluid Mech. 42 (3) (1970)
471–489, h ps://doi.o g/10.1017/S0022112070001428, ci ed by: 777, h ps://www.scopus.com/inwa d/ eco d.u i?eid=2-s2.0-0014938217&doi=10.1017%
2 S0022112070001428&pa ne ID=40&md5=adcab63a521b38ba40ad46d52736529e.
[47] B. Fo nbe g, A nume ical s udy o s eady iscous flow pas a ci cula cylinde , J. Fluid Mech. 98 (4) (1980) 819–855.
[48] J.-I. Choi, R.C. Obe oi, J.R. Edwa ds, J.A. Rosa i, An imme sed bounda y me hod o complex incomp essible flows, J. Compu . Phys. 224 (2) (2007) 757–784,
h ps://doi.o g/10.1016/j.jcp.2006.10.032.
[49] B. Fo nbe g, S eady iscous flow pas a sphe e a high Reynolds numbe s, J. Fluid Mech. 190 (1988) 471–489, h ps://doi.o g/10.1017/S0022112088001417.
[50] O. Boi on, G. Chia assa, R. Dona , A high- esolu ion penaliza ion me hod o la ge Mach numbe flows in he p esence o obs acles, Compu . Fluids 38 (3) (2009)
703–714, h ps://doi.o g/10.1016/j.compfluid.2008.07.003.
[51] P. Ango , C.-H. B uneau, P. Fab ie, A penaliza ion me hod o ake in o accoun obs acles in iscous flows, Nume . Ma h. 81 (1999) 497–520, h ps://doi.o g/10.
1007/s002110050401.
[52] F. Bassi, S. Rebay, A high-o de accu a e discon inuous fini e elemen me hod o he nume ical solu ion o he comp essible Na ie –S okes equa ions, J. Compu .
Phys. 131 (2) (1997) 267–279, h ps://doi.o g/10.1006/jcph.1996.5572.
Jou nal o Compu a ional Physics 528 (2025) 113807
20
