Keywords:

A Transport and shape calculus theorems 151 B Cut Cell method for Navier-Stokes equations 153 C Flower.jl package description 159 D Spurious currents study for static drops 163 E Numerical implementation of the super-slip model 169 F Kinematic transport of the contact angle 173 G Relaxation of the impermeability condition in the toy model 177

te tiens à remerierD tout d9ordD l9snstitut des ienes du glul et des honnées et son direteurD sl preyD de m9voir permis de réliser ette thèseD et fruno hesprés d9voir epté de fire prtie de ette ventureF wes remeriements vont églement à l9éole dotorle et son direteurD hjimedo uondoD insi qu9à l9snstitut ten le ond d9elemert dns son ensemleF te remerie ernud entkowikD wri elntiD henis ipp et qretr ryggvson pour leurs ommentires et questions qui ont été d9une grnde préision et ont ouvert l voie à de nomreuses extensions possiles de e trvilF te suis extrêmement reonnissnt ux rpporteurs de ette thèseD ghristophe tossernd et eter hmidD qui m9ont suivi tout u long de es trois dernières nnées et qui ont su iguiller mes reherhesD notmment pendnt les périodes de on(nement où tout semlit ien étrngeF weri à rneh ydi quiD dès notre première réunion dns son ureuD su diriger ve rio et mitieux projet de reherheF te sis que l (n de l thèse ne mrquer en uun s l (n de notre ollortion et j9en suis très heureuxF wes remeriements vont ussi à inent ve ghendeD qui été un exemple de jeune herheur pssionné et uprès de qui j9uris énormément pprisF weri à téphne leski de m9voir fit on(ne il y ientôt qutre ns et d9voir étéD depuis e momentElàD un mentor pour moiF te ne suris mnquer de remerier qustv emerg et les memres du groupe ur pour les préieuses disussions sur l ligne de onttF te tiens ussi à remerier tous mes professeurs etD en prtiulierD rélène humontetD gtherine eismn et tenEgmille ghssing pour leur soutien sns fille depuis le déut de mon ursusF te remerie shD wthisD elexndreD tous mes ollègues du pulse et de d9elemert plus lrgementD ve qui j9i trvillé et pssé euoup de ons moE mentsF weri à erenD glément et elejndro pour leur sympthieD leur onne humeur et leur mitiéD sns oulierD ien entenduD qonzlo qui m9 supporté penE dnt toutes es nnées en fe de luiF te remerie ploD vouisD vuile et tous mes mis de m9voir fit reltiviserD voire oulier mes petits prolèmes de dotorntF te remerie églement m fmille de ÷urD prisienneD qui se reonnîtr pr es lignesF weri à mes onlesD tntesD ousinsD ousines et à m grndEmère en qui m9ont toujours soutenu même si un hémisphère nous sépreF t9i une pensée pour mes grndsEprents wríEvuí et elerto à qui je dois euoupF weri à mes prentsD grolin et toséD et à m s÷urD imiliD d9voir toujours été à mes ôtésF te dédie ette thèse à rnnhD à quiD un seul meri ne su0rit ps à exprimer toute l grtitude et l9mour que je lui porteF Simulation et optimisation de phénomènes complexes dans des écoulements multiphasiques Résumé: hns ette thèseD nous étudions des phénomènes omplexes présents dns les éoulements multiphsiques en utilisnt diverses méthodes et modèles numériquesF v9ent est misD tout d9ordD sur le prolème diphsique de tefnD un prolème de hngement de phse où le mouvement de l9interfe est lié u sut des )ux ondutifs de hleurF ves formes omplexes de l9interfe dns l roissne ristllineD en présene d9e'ets d9nisotropieD sont prédites ve préision grâe à notre méthode d9ensemle de niveux ouplée à une nouvelle méthode de ellules oupées pour le trnsport di'usifF ne proédure d9optimistion de l forme de ette interfeD sée sur l dérivtion de l9djoint sous forme ontinueD ve une fontion de oût de type suivi d9interfeD est ensuite onstruite et testée sur di'érentes on(gurtions physiquesF in onsidérnt les équtions de xvierEtokes dns l9pproximtion de foussinesqD une étude sur l9pprition de l9instilité de yleighEfénrd en présene d9une frontière de fusion est églement réliséeF hns l deuxième prtieD nous étudions les prolèmes de lignes de ontt dynmiquesD où une interfe )uideE)uide est en ontt ve une frontière solideF lusieurs nouveux modèles sont mis en ple u sein d9une méthode du volume de )uideD tel que l ondition limite générlisée de xvierD où l vitesse de l ligne de ontt dépend de l ontrinte de oung non ompenséeD ou enore le modèleEjouetD où le )ux mssique existnt dns les méthodes à interfe di'use est mis en reltion ve l ourure de l9interfeD résultnt en une simple reltion d9ngle de ontt dynmiqueF Mots-clés: prolèmes de tefnD djoint ontinuD méthode d9ensemle de niveuxD méthode de ellules oupéesD lignes de onttD méthode du volume de )uide Simulation and optimization of complex phenomena in multiphase ows Abstract: sn this thesisD we study omplex phenomen in multiphse )ows using vrious new numeril methods nd modelsF he fous isD (rstD set on the twoEphse tefn prolemD phseEhnge prolem where the motion of the interfe is relted to the jump in ondutive het )uxesF he omplex interfil shpes in rystl growthD in the presene nisotropy e'etsD re urtely predited thnks to our level set method oupled with novel gut gell method for di'usive trnsportF e shpe optimiztion proedure with trkingEtype ost funtionl sed on ontinuous djoint derivtion is then uilt nd tested on vrious physil setupsD yielding fvorle resultsF gonsidering the xvierEtokes equtions in the foussinesq pproximtionD study of onset of yleighEfénrd instility in the presene of melting oundry is lso rried outF sn the seond prtD we investigte dynmi ontt line prolemsD where )uidE)uid interfe meets solid oundryF rious novel methods re implemented in olumeEyfEpluid frmeworkD suh s the generlized xvier oundry onditionD where the speed of the ontt line depends on the unompensted oung9s stressD or the toy modelD where the existing mss )ux in di'use interfe models is relted the interfe urvtureD resulting in simple dynmi ngle reltionF Introduction 1

I Shape optimization of melting and solidication processes 6

List of Figures IFI hemti of the liquid nd solid sudominsF F F F F F F F F F F F F F U IFP hotogrphy of rystl growing in n underEooled liquid thF F F W IFQ iew of the level set signed distne funtion ϕF F F F F F F F F F F F F IH IFR hemti of the twoEphse tefn prolem in the level set frmeworkF II IFS hotogrphy of n tmospheri surfe lyerF F F F F F F F F F F F F F IP IFT hemti of the melting oundry prolemF F F F F F F F F F F F F F F IQ PFI hemti of the interfe trking term in the ost funtionlF F F F IU QFI ixt nd pproximte representtion of ell ut y the interfeF F QP QFP voss of ury in meshEligned geometriesF F F F F F F F F F F F F F F QR QFQ inhned geometri informtion restoring uryF F F F F F F F F F QS QFR sntersetion of the )uid domin with grtesin elementsF F F F F F F F QV QFS hemti of the (rst kind pitiesF F F F F F F F F F F F F F F F F F F QW QFT hemti of the seond kind pitiesF F F F F F F F F F F F F F F F F RH QFU hemti of the di'erent ses in the mrhing squres lgorithmF F RI QFV gonvergene of the gut gell method when solving inside irleF F F RT QFW gonvergene of the gut gell method when solving outside irleF F RU QFIH gonvergene of the gut gell method when solving inside squreF F RV QFII gonvergene of the gut gell method when solving inside rystlF F RW RFI hemti of the tohnsenEgolell method for the grdient lultionF SR RFP gonvergene of the grdient omputtion in sttionry geometryF F SS RFQ ixtension of the veloity (eld for di'erent nrrow nd widthsF F F F ST RFR hemti of the dimondEell strtegy to ompute the rmiltoninF SV RFS gonvergene of the porwrdEfkwrd di'usion shemeF F F F F F F F TI RFT einitiliztion of n initil pertured signed distne funtionF F F F TQ RFU hemti of the tretment of the fresh nd ded ellsF F F F F F F F F TR SFI gonvergene of the temperture (eld in the plnr interfe seF F F TW SFP gonvergene of the rdius in the prnk sphere seF F F F F F F F F F F UH SFQ emperture (eld error mp in the prnk sphere seF F F F F F F F F UI SFR gonvergene of the temperture (eld in the prnk sphere seF F F F UP SFS xumeril instilities in the sene of surfe tension e'etsF F F F UQ SFT egulriztion introdued y the surfe tension e'etsF F F F F F F F UR SFU i'et of the surfe tension oe0ient on dendriti growthF F F F F F US SFV i'et of the grid re(nement on dendriti growthF F F F F F F F F F F F UT SFW i'et of the moleulr kinemti oe0ient on dendriti growthF F F UU SFIH enisotropy e'ets on rystl growth with α 0 = π/2F F F F F F F F F F UW SFII enisotropy e'ets on rystl growth with α 0 = π/4F F F F F F F F F F VH vi List of Figures SFIP i'etive yleigh numer for di'erent glol yleigh numersF F F VP SFIQ everge height for di'erent glol yleigh numersF F F F F F F F F F VQ SFIR imes series of the temperture (eld with Ra = 10 6 F F F F F F F F F F VR TFI gonvex ost funtionl s funtion of two prmeters a 1 nd b 1 F F F VT TFP yptimiztion results nd interfe shpes in the melting irle seF F VU TFQ gost funtionl vrition in the melting irle seF F F F F F F F F F F VV TFR yptimiztion results nd interfe shpes in the wullinsEekerk seF VW TFS gost funtionl vrition in the wullinsEekerk seF F F F F F F F F F WH TFT yptimiztion results nd interfe shpes in the growing rystls seF WI TFU gost funtionl vrition in the growing rystls seF F F F F F F F F WP UFI hensity nd visosity representtion in the oneE)uid formultionF F F WT VFI hemti of the stti ontt lineF F F F F F F F F F F F F F F F F F F F IHI WFI hemti of the Ph heightEfuntions ner the ontt line F F F F F F IHT WFP hemti of the qudtree disretiztion nd tree representtionF F F IHV WFQ hemti of geometril )ux estimtionF F F F F F F F F F F F F F F F F IHV WFR hemti of the xvier oundry onditionF F F F F F F F F F F F F F IHW WFS hemti of the urtin oting on(gurtionF F F F F F F F F F F F F IIP WFT ime series of urtin oting simultion with e a IPD g a IFQQF IIP WFU qrid resolution ner the ontt line with λ = 10 µmF F F F F F F F F F IIR WFV tedy stte solutions for the three )ow on(gurtionsF F F F F F F F F IIR WFW tility nlysis of the redued urtin oting systemF F F F F F F F IIS WFIH gonvergene of the ontt line nd in)etion point positionsF F F F F IIT WFII tility nlysis nd omprison with the flke experimentsF F F F F IIV WFIP tility nlysis nd omprison with the wrston experimentsF F F IIV WFIQ gomprison with the vvi experiments for spreding dropF F F F F F IPH IHFI dius s funtion of time for the superEslip modelF F F F F F F F F F IPQ IHFP vol pillry numer s funtion of time for the superEslip modelF IPQ IHFQ gontt ngle s funtion of time for the superEslip modelF F F F F F IPR IHFR her stress s funtion of time for the superEslip modelF F F F F F F IPR IHFS hemti of the ngle extrpoltion using the qxfgF F F F F F F F F IPU IHFT tedy stte menisus for Ca = 0.03 using the qxfgF F F F F F F F F IPV IHFU gontt line position s funtion of time using the qxfgF F F F F F IQI IHFV gontt ngle s funtion of time using the qxfgF F F F F F F F F F IQP IHFW her stress s funtion of time using the qxfgF F F F F F F F F F F IQQ IHFIH gomprison etween the noEslip model nd the qxfgF F F F F F F F IQR IHFII hemti of the p nd yp omputtionl dominsF F F F F F F F IRH IHFIP preding drop dimeter omprison of the p nd toy modelF F F F IRI IHFIQ preding drop interfes of the p nd toy modelF F F F F F F F F F IRP IHFIR oomed interfes of the p nd toy modelF F F F F F F F F F F F F F IRQ IHFIS ummry of the di'erent models for the dynmi ontt lineF F F F IRS Introduction he study of )uid )ows hs een lsting humn endevorD pursued to gin further understnding of xture nd industril proessesF prom mthemtil point of viewD these systems re modeled y prtil di'erentil equtions where the unknowns9 re oth spe nd time dependentF hue to the highly omplex nture of the phenomen rising from these )ows ! even in seemingly simple systems ! these equtions re prtilly unsolvle y nlytil tehniquesF xumeril nlysisD howeverD is n pproh tht is well suited to study these systemsD y using disrete pproximtions of the prtil di'erentil equtionsF sn the ontext of multiphse )ows ! s suset of )uid )ows orresponding to the simultneous )ow of mterils with two or more thermodynmi phses ! the min di0ulty omes from the high nonElinerity introdued y the interfe nd spei( numeril methods re required to trk the interfe with uryF sn this work we will study di'erent physil systems using severl originl numerE il methods nd modelsF sn the (rst prtD we will onsider phseEhnge prolems ! soElled tefn prolems ! involving therml energy exhnges etween solid nd liquid phseF sn the seond prtD we will investigte prolems involving dynmi ontt linesD where )uidE)uid interfe moves over solid sustrteF Two-phase Stefan problems tefn prolemsD nmed fter the eustrin mtheE mtiin toseph tefn tefn IVWID rler IWWS due to his sustntil ontriution to reserh on moving nd free oundriesD model trnsport nd trnsfer phenomenD in prtiulr solidEliquid phse hnge in evporting or hemilly reting )owsF uh phenomen govern the interfe motion in mny engineering relted prolems suh s dendriti solidi(tion ysher 8 ethin IWVVD turi 8 ryggvson IWWTD phse trnsformtion in metlli lloys egl et al. IWWVD nd solid fuel omusE tion rssn et al. PHPIF sn the pplitions of interest to this studyD the tefn ondition rises from the intertion of liquid nd solid phses ! oth onsidered inompressile ! resulting in moving liquidEsolid interfeD typilly freezing or melting frontF he speed of the front is diretly relted to the jump in the ondutive het )ux ross the interfeF sn one dimensionD this prolem hs een studied in depth using vrious numeril lgorithms frttkus 8 weiron IWWPD ose IWWQD tvierre et al. PHHTF sn higher dimensionsD vrious methods hve een used suh s the level set method ysher 8 ethin IWVVD vimre et al. PHPQ nd frontEtrking method turi 8 ryggvson IWWTF yne of the min hllenges ssoited with modeling suh prolems in multiple dimensions is due to the unstle dentriti pttern formtion wullins 8 ekerk IWTRD vnger IWVHD oods IWWPF sn rystl growthD for exmpleD underEooling triggers n instility mehnismD using the solid phse of the mteril to grow into the liquid phse in (ngerElike or dendriti fshionD resulting in omplex interfil shpesD whih re hllenging to predit Introduction numerillyF sn dditionD prmeters suh s nisotropy nd surfe tension or urE vture e'et @qisEhomsonA re shown to hve lrge impt on the dendriti shpe of the rystlD whih in turn need to e modeled urtely for the simultions to remin preditiveF sn this workD we present generl frmework for trking nd modeling rystl growth in the presene of urvture e'etsF his lgorithm then serves s vessel to mterilize the seondD nd minD ojetive of this study whih is extrting optimiztion strtegies to ontrol the resulting solidi(tion proessF he shpe of the interfe strongly e'ets the outome nd timeEfrme of the prodution proess in mny industril pplitions involving phse hngeF es resultD while prediting nd modeling the resulting solidi(tion proess remins t the forefront of mny reserh resD it is s desirle to extrt e0ient ontrol strtegies to mnipulte the motion of the interfeD for instneD y trking presried trjetoryF wo mjor types of optimiztion methods in use tody re @iA grdientEsedD nd @iiA derivtive free methodsF hile n e0ient lss of generi lgorithms ! eE longing to the lss of derivtive free methods ! sed on the surrogte mngement frmework wrsden et al. PHHV nd rti(il neurl networks ierret et al. PHHU hve een used for optimiztion in )uid mehnisD minly in the re of eroE dynmi shpe optimiztionD these often require mny funtion evlutionsD for trining purposes for exmpleF hen detiled simultions of interfil )ows re onernedD eh funtion evlution ommnds full ! potentilly unstedy ! gomE puttionl pluid hynmis @gphA simultionD using grdientEsed methods to e t n dvntgeF wost ommon methods in extrting the grdient informtionD on the other hndD re nlytil or use (nite di'erenesD neither of whih re suited to the on(gurtion of interest to this workF edjointEsed lgorithms present suitle lterntiveD s they llow the determintion of the grdient t ost omprle to single funtion evlution qiles 8 iere PHHHF he use of djoint methods for design nd optimiztion hs een n tive re of reserh whih strted with the pioneering work of ironneu ironneu IWUR with pplitions in )uid mehnisD nd hs een extensively used in eronutil shpe optimiztion y tmeson nd oEworkers tmeson IWVVD tmeson et al. IWWVF iver sine these groundreking studiesD djointEsed methods hve een widely used in )uid mehnis prtiulrly in the res of eroE nd thermoEoustis tuniper PHIHD vemke et al. PHIQF wore reently )ow regimes dominted y nonliner dynmisD suh s seprtion nd mixing hve lso een nlysed using djointEsed tehniques hurismy 8 elonso PHIPD hmidt et al. sn the ontext of tefn prolemsD vrious ontrol strtegies hve een emE ployed to trk the lotion of the interfeF sn oneEdimensionl settingD for exmpleD setEvlued (xed point equtions ro'mnn 8 prekels IWVP or linerE qudrti defet minimiztion unner IWVS hve een used to ontrol the loE tion of the frontF edjointEsed lgorithms hve lso een pplied to teE fn prolem used y heterogeneous retions on surfe of oneEdimensionl solid prtile rssn et al. PHPI to extrt sensitives with respet to vrious kiE neti prmetersF elterntivelyD in twoEdimensionl settingD djointEsed lgoE rithms hve een utilized previously together with (nite element nd (nite di'erE ene pprohes to trk nd ontrol the lotion of interfe y imposing het )ux @or tempertureA t the oundry in order to relize the desired interfe moE tion ung 8 rs IWWSD ng IWWUD rinze 8 iegenlg PHHUF sn prtiulrD fernuer nd rerzog fernuer 8 rerzog PHIID mking use of shpe lulus toolsD derived the set of djoint equtions to extrt ontrol strtegies for tefn prolems with shrp representtion of the interfeF imilr to the pproh of fernuer 8 rerzog PHII shpe lulus tools hve een employed in this work to extrt the orresponding djoint equtionsF roweverD ontrry to the previous studiesD ontrol strtegies re extrted here to suppress instilities ! s dentriti formtion ! of the solidi(tion proessF he ddition of urvture e'ets on the interfe ! qisEhomson reltion ! nd the omplex shpes enountered during the growth of the rystl require dedited numeril lgorithms ple of modeling oth the highly nonliner forwrd prolem ! orE ret representtion of the resulting interfe ! nd the resulting djoint equtionsF sn dditionD while previous studies mostly onentrted on tution y imposing temperture or het )ux t the oundriesD lterntive ontrol strtegies using the length of interfe re lso investigted hereD to identify the most relevnt nd e'etive numeril ontrol strtegies in the ontext of melting nd solidi(tion proesses pulln et al. PHPQF es n extension of the lssil twoEphse tefn prolemD we will onsider the onvetive motion indued y )ow in the )uid phse uirós odríguez et al. PHPP in di'erent setting thn the djointEsed optimiztion frmework de(ned previE ouslyF sn tht seD the xvierEtokes equtions in the foussinesq pproximtion ! de(ned lter ! will e solved in order to model the uoyny phenomenonD where the )uid exerts fore tht opposes the weight of solid phseF sn prtiulrD we will study the melting proess of solid where n initil )uid lyer is heted from elowF his motion will entrin the onset of yleighEfénrd instility ! where onvetion ells pper depending on the wve length of the lyer ! nd prmetri study on the e'et of the yleigh numerD similr to the ones of pvier pvier et al. PHIW nd vimre vimre et al. PHPQD will e rried outF Dynamic contact lines he ontt line formed t the intersetion of )uidE)uid nd solid oundry ! soElled the triple point ! is key element of mny ntuE rl nd severl tehnologil proessesD s desried in flke 8 hikhmurzev PHHPD ilson et al. PHHTD wrston et al. PHHWF he dynmi ontt line poses ! sine the erly work of ruh nd riven ruh 8 riven IWUI ! remrkle prolem Introduction euse of the ontrdition etween the noEslip ondition on the sustrte nd the motion of the ontt lineF e (xed ontt lineD for exmpleD is in ontrdition with simple oservtions of spreding drop on surfeF wthemtillyD this prdox introdues stress singulrity t the ontt line nd the numeril models to propE erly predit the motion of the triple point re still deted to this dy gox IWVTD hikhmurzev IWWRD flke et al. IWWWD iggers PHHRD iggers 8 ivns PHHSF sn the ontext of ontinuum mehnisD in shrpEinterfe modelsD the simE plest wy to relx the singulrity is to introdue xvier oundry ondition @xfgA ! y llowing the ontt line to slide over the sustrte ! hrterE ized y slip length of the order of the moleulr lengthF his model hs een proven to urtely reprodue the motion of the ontt line in vrious physil setups y either (tting the slip length with moleulr dynmis reE sults v© is et al. PHPHD v© is et al. PHPP or hoosing n pproprite slip relted to experimentl oservtions pulln et al. PHPHF enother essentil prmeter to onsider is the vrition of the ngle formed y the )uidE)uid interfe with the solid ! the dynmi ontt ngle ! s the ontt line movesF sn the lssil goxEoinov reltion oinov IWUTD gox IWVTD the dynmi ngle is funtion of the equilirium ngle ! property of the sustrte ! nd the ontt line speedF snspired y the work of vegendre nd wglio vegendre 8 wglio PHISD where they tested di'erent dynmi ngle modE elsD we will implement nd enhmrk them ginst experiments of spreding squlne drop vvi 8 wrmur PHHRF sn n ttempt to further understnd the physis involved in this prolemD we develop new nd originl numeril methods for dynmi ontt linesF he two previously desried models ! xfg nd dynmi ngleD lled lssil models in this disserttion ! will serve s sis for the implementtion of the dvned onesF sn the superEslip modelD further regulriztion of the ontt line is provided through seondEorder slip lengthF his modelD in nlogy with porous medi theory hevuhelle et al. PHHUD llows relxtion of the pillry pressure present in the xfgF e study on the e'et of the seondEorder slip in the se of spreding drop will e rried outF he seond dvned model is the generlized xvier oundry ondition @qxfgAF fy onsidering the e'et of the unompensted oung9s stress we llow n ngle devition with respet to the equilirium ngleF his oundry onditionD originlly derived for di'use interfe models in et al. PHHQD in et al. PHHT hs een proven to e wellEposed in the shrp interfe limit prike PHPHF he qxfg will e implemented in our olumeEyfEpluid frmework nd we will study its e'et on the onset of wetting filure ! where menisus trnsitions to liquid (lm ! in withdrwing plte setupF pinllyD we will use toy model to estlish reltion etween the mss )ux present in the hseEpield @pA model tqmin PHHHD emerg PHHQD grlson PHIP to urvture reltion t given distne of the ontt lineF hrough simple oundry onditions in yp modelD we will e le to reprodue the di'usionE driven wetting motion existing in the pD t muh lower ostF Original contributions in this work yne of the mjor ontriutions of our work is the derivtion of the ontinuous djoint twoEphse tefn prolem in the presene of urvture e'ets nd the onstrution of the orresponding djointEsed lgorithmF he numeril methods desried in ghpters Q nd R onstitute novE elty with respet to the existing onesD in prtiulrD the gut gell method for di'usive trnsportF he forementioned methods nd test ses re oded in Flower.jl @see ghpter S nd eppendix g for further detilsAD new tuli pkge developed s prt of this workF ell the ontt line models ! xfgD dynmi ontt ngleD superE slipD qxfgD toy model ! re implemented in BasiliskD n existing pltform for the solution of prtil di'erentil equtions @presented in ghpter WAF Outline he (rst prt of this work will e foused on the optimiztion of twoE phse tefn prolemsF sn ghpter ID the governing equtions in ontinuous form in level set frmework will e desriedF edditionllyD the )uid )ow extension will e desriedF sn ghpter PD we present the derivtion of the ontinuous djoint of the tefn prolem in presene of surfe tension e'ets nd the resulting djointE sed optimiztion lgorithmF sn ghpter QD the newly developed gut gell method for di'usive trnsport nd its numeril implementtion will e presented in depthF ghpter R presents the rest of the numeril methods used to solve oth the forwrd nd djoint tefn prolemsF sn ghpter SD numeril vlidtion of our sheme ginst nlytil solutions will e providedF woreoverD we will rry out studies on dendriti rystl grows nd on the onset of yleighEfénrd onvetion ells with melting oundryF pinllyD in ghpter TD we will present shpe optimiztion ses ! with vrile omplexity ! of melting nd solidi(tion proesses nd ompre the results otined with the djointEsed lgorithm with derivtiveEfree oneF sn the seond prtD fous will e set on the dynmi ontt lineF sn ghpE ter UD the oneE)uid formultion of the twoEphse xvierEtokes equtions nd the hseEpield extensions will e desriedF ghpter V disusses the prdoxes nd singulrities present in dynmi ontt line prolemsF sn ghpter WD the olumeE yfEpluid method nd the implementtion of the lssil models for ontt lines re desriedF heses models will e tested on di'erent physil systemsF pinllyD in ghpter IHD the newly developed dvned models for ontt lines will e detiled ndD ginD tested on vrious physil prolemsF Part I

Shape optimization of melting and solidication processes

Chapter 1

The two-phase Stefan problem e tefn prolem is spei( type of free oundry prolem9 hrterized y the temperture distriution of two di'erent phsesF he position of the interfe etween those two phses is a priori unknown nd will e prt of the solution of prtil di'erentil equtions desriing this systemF es (rst pproximtionD we will onsider the twoEphse tefn prolem in the presene of two immisile phses with mthing densities @one liquid nd the other solidAF e de(ne domin Ω prtitioned into two timeEdependent sudomins Ω 1 (t) nd Ω 2 (t) oupied y the liquid @1A nd solid @2A phsesD respetivelyF he externl oundry of the dominD denoted ∂ΩD is (xed wheres the interfe seprting oth phses Γ (t) = Ω 1 (t) ∩ Ω 2 (t) evolves in timeF e shemti of this on(gurtion is shown in pigure IFIF Liquid phase he energy trnsport mehnism in eh phse is the het produed y difE ferene in tempertureF he prtil di'erentil eqution governing this quntity is the het equtionF vet T i : (x, t) → R + denote the temperture (eld in either Chapter 1. The two-phase Stefan problem phse @i = 1, 2AD where x = (x, y) is the grtesin oordinte vetorF hen the densities ρ 1 nd ρ 2 re equlD onvetive het trnsfer vnishesF sn dditionD when the kground pressure is ssumed onstntD the het trnsport eqution simpli(es to

Ω 1 (t) Solid phase Ω 2 (t) Interface Γ(t)
∀i ∈ 1, 2 , ρ i c i ∂T i ∂t = ∇ • (k i ∇T i ) , t > 0, x ∈ Ω i , @IFIA
where ρ i D c i nd k i denote the densityD the spei( het pity t onstnt pressureD nd the therml ondutivityD for eh phseF es the phse hnge phenomen oursD there will pper ltent het whih either is sored or relesedF he ondition of het onservtion t given point on the moving interfe orresponds to the rte t whih het is generted t the oundry lned y the rte t whih this het )ows in either phseF hereforeD long the interfeD energy lne sttes tht

v S = -L 1 H k 1 ∂T 1 ∂n -k 2 ∂T 2 ∂n , x ∈ Γ, @IFPA
where L H denotes the ltent het of solidi(tion nd n is the outwrd norml unit vetor t the interfeF his jump is tken from phse I to phse P with ∂T i /∂n denoting the norml omponent of the temperture grdient in phse iF iqution IFP is ommonly referred to s the tefn onditionF sn lssil tefn prolemD the temperture is set to T D (x, t) = T M t the inE terfe where T M is onstnt equl to the melting temperture of the mterilF por prolems involving rystl growth however vnger IWVHD surfe tension e'ets must e dded to the thermodynmi oundry ondition y introduing depenE dene in the urvture t the frontF por exmpleD in the formtion of rystls @pigure IFPAD redution of the melting temperture for smll prtiles through the surfe tension e'ets provides the neessry stilizing fores llowing the retion of ptE ternsF woreoverD from moleulr kineti theoriesD we know tht the melting temperE ture might lso depend on the speed t whih the interfe movesF por tht purposeD we use the lssil qisEhomson reltionD s de(ned in ghen ghen et al. IWWU

T D (x, t) = T M -ε V v S -ε κ κ, x ∈ Γ, @IFQA
where κ denotes the urvture t the interfe @positive if the enter of urvture lies in the solid phseAD v S the veloity of the interfeD ε κ the surfe tension oe0ientD nd ε V the moleulr kineti oe0ientF nless stted otherwiseD in the rest of this disserttionD oth ε κ nd ε V re onsidered to e onstntsD nd the het pitiesD therml ondutivities nd ltent het re ll set to unityF sn dditionD to ese the nottionD the jump in grdient of temperture is denoted s [∇T ] 1 2 @lso de(ned in iqution IFPAF his set of equtions hrterize the twoEphse tefn prolemF hen solving the system numerilD the min di0ulty rises from the tretment of the interfe due to the nonEliner reltion etween the tefn ondition nd the two domins where the het equtions re solvedF 1.2 The level set method hen deling with the numeril pproximtion of interfil )owsD two lsses of methods re ommonly used to represent the interfeD nmely vgrngin or frontEtrking9 methodsD iulerin or frontEpturing9 methodsF he former uses prmeteristion of the interfe lotion @eFgF mrkers or moving meshesAD nd hs lredy een used in tefn prolems turi 8 ryggvson IWWT ut hs not een dopted in this work due to the inherent di0ulty of deriving the ontinuous djoint equtions with suh methods @see ghpter PAF he ltter n rodly speking e divided into two tegoriesX olumeEyfEpluid @ypA nd 10 Chapter 1. The two-phase Stefan problem vevel et methodsF he yp method will e desried in the seond prt of this disserttion when deling with moving ontt lines @see ghpter WAF sn djointEsed optimiztionD the yp method pikl et al. PHPH my led to numeril omplitions due to the pieeEwise liner reonstrution of the interfeF yn the other hndD the level set methodD where the interfe is impliitly de(ned s ontinuous funtionD hs een proven to e well suited in ontinuous djointEsed optimiztionD spei(lly for tefn prolems fernuer 8 rerzog PHIIF woreE overD this method hs een shown to urtely reprodue dendriti pttern forE mtion ghen et al. IWWUD vimre et al. PHPQF en impliit signed distne strtegy hs therefore een used here to represent the interfeF e level set funtion ϕ : (x, t) → R is onstrutedD suh thtD t ny time tD the front is equl to the zero level set of the funtion

Γ(t) = {x ∈ Ω : ϕ(x, t) = 0}. @IFRA
he level set funtion is initilly set to the signed distne funtionD with d the distne to the frontD suh tht sn prtieD the level set funtion is funtion of dimension D + 1D where D is the dimension of the prolem onsidered @D = 2 in this seAF he intersetion etween ϕ nd the HElevel plne orresponds to the interfe @pigure IFQAF yne of the dvntges of this method is the diret de(nition of the norml vetor nd the urvtureD quntities tht re neessry to ompute the tefn ondition nd the qisEhomson reltionF he norml unit vetor n is de(ned y

ϕ(x, 0) =    +d, x ∈ Ω 2 , 0, x ∈ Γ, -d, x ∈ Ω 1 . @IFSA
n = ∇ϕ(x, t) |∇ϕ(x, t)| , x ∈ Γ (t) , @IFTA nd the urvture κ y κ = ∇ • ∇ϕ(x, t) |∇ϕ(x, t)| , x ∈ Γ (t) . @IFUA
sn two dimensionsD the urvture κ t the front is omputed in nonEonservtive form using

κ = ϕ 2 y ϕ xx -2ϕ x ϕ y ϕ xy + ϕ 2 x ϕ yy ϕ 2 x + ϕ 2 y 3/2 , @IFVA
where ϕ x D ϕ xx D ϕ y D ϕ yy denote the (rst nd seond derivtive of ϕ with respet to x nd yF he interfe is moved through the level set dvetion equtionD where F is speed funtion de(ned everywhere

∂ϕ ∂t + F |∇ϕ| = 0. @IFWA ϕ(x, t) < 0 T 2 (x, t), x ∈ Ω 2 ϕ(x, t) > 0 T 1 (x, t), x ∈ Ω 1 ∂Ω v S = -[∇T i ] 1 2 T D = T M -ε V v S -ε κ κ ϕ(x, t) = 0 Γ Interface :
n sn the tefn prolemD the veloity is only de(ned t the interfe @iqution IFPAF he veloity is numerilly extended in the norml diretion of the interfeD round nrrow nd9 @de(ned in ghpter RAF fy omining iqutions IFP nd IFW nd using the norml de(nition @iqution IFTAD we otin the following eqution of motion

∂ϕ ∂t = -[∇T ] 1 2 • ∇ϕ. @IFIHA
iqution @IFIHA will move ϕ with the speed suh tht Γ remins lose to the zero level set of ϕF sn prtieD nd this is one of the min drwks of this methodD the level set funtion needs to e reinitilized suh tht it retins its signed distne funtion propertiesF he resulting twoEphse tefn prolem in the level set frmework is desried in pigure IFRF 1.3 Thermal convection and uid ow extension he model de(ned so fr ssumes tht the )ow in the )uid phse n e negletedF his is n pproprite hypothesis when onsidering equl densities in the )uid nd solid phsesF roweverD in some pplitionsD this ssumption is too restritiveF sn pigure IFSD n exmple of yleighEfénrd onvetion ells ppering in lyer moving upwrds is shownF o model suh phenomenD one hs to tke into ount oth onvetion nd the xvierEtokes equtions in the foussinesq pproximtion in the )uid phseF sn this workD we will fous on the sme setup s pvier et al. PHIWD where they onsider the evolution of horizontl lyer heted from elow nd omprised etween two wlls t distne H prtF he grvity is pointing downwrds g = -ge y nd the horizontl size of the domin is bH where b is the spet rtioF he temperture t the ottom wll is T = T 1 wall nd the one t the top wll is T = T 0 wall suh tht T 0 wall < T M < T 1 wall where T M is the melting tempertureF en shemti of this on(gurtion is shown in pigure IFTF he governing dimensionless equtions in the )uid phse re now

Pr -1 ∂u ∂t + u • ∇u = -∇P + Ra T e y + ∇ 2 u, ∇ • u = 0, ∂T ∂t + u • ∇T = ∇ 2 T, @IFIIA 1.3.
Ra = gα t ∆T H 3 νk 1 , @IFIQA
where g is the onstnt grvittionl elertionD α t the oe0ient of therml exE pnsion nd ∆T = T 1 wall -T 0 wall the di'erene in temperture etween the two pltesF sn the solid phseD when onsidering the nonEisotherml generl seD with T 0 wall ̸ = T M D the het eqution pplies

∂T ∂t = ∇ 2 T. @IFIRA
sn this on(gurtionD we will e le to study the onset of the yleigh onvetion ells nd the hnges in topology s the interfe melts through the solid @eE tion SFRAF sn the next hpterD we present the optimiztion proedure involving the derivtion of the ontinuous djoint of the twoEphse tefn prolem in the sene of )uid )owF

: Ω × [0, t f ] → R nd funtion ϕ : Ω × [0, t f ] → R suh tht                                                          ∂T 1 ∂t = ∇ 2 T 1 in Ω 1 (t) @pFA ∂T 2 ∂t = ∇ 2 T 2 in Ω 2 (t) @pFA T (x, 0) = T 0 (x) in Ω @pFA ∂T (x, t) ∂n = w (x, t) on ∂Ω @pFdA T (x, t) = T M -ε V v S -ε κ κ on Γ(t) @pFeA ∂ϕ ∂t = -[∇T i ] 1 2 • ∇ϕ on Γ(t) @pFfA ϕ(x, 0) = ϕ 0 (x)
in Ω @pFgA @pA rereD T i denotes the restritions of T to Ω i (t)D nd T 0 nd ϕ 0 re the initil onE ditions t t = 0 of the temperture (eld nd the level set funtionF he (rst two equtions pF nd pF orrespond to the het equtions in phses I nd P respeE tivelyF iqutions pF nd pFg orrespond to the initil temperture distriution in the whole domin nd the initil position of the front respetivelyF he equE tion pFd orrespond to the ontrol vrile wD the het )ux ting on the domin oundryF iqution pFe is the hirihlet oundry ondition t the interfe with the qisEhomson reltion tken into ountF pinllyD pFf orresponds to the dvetion eqution of the level funtion s de(ned in iqution IFIH where -[∇T i ] 1 2 is the tefn onditionF he next step is to de(ne the desired temperture (eld T d nd the desired position of the interfe ϕ d F he ontrol vrile w is then de(ned trough the optimiztion proedure to reh these desired quntitiesF hese onditions re known priori nd will drive our ost funtionl towrds its minimumF he following trkingEtype ost funtionl provides mthemtil desription of the ontrol gols stted ove he (rst term monitors the (nl temperture distriution nd is mostly used s n initilizer for the djoint temperture (eld for the djoint prolem @see etion PFPAF he seond term ontrols the reltive position of the level set funtion with respet to the desired oneF king dvntge of the signed distne funtion property of oth ϕ d nd ϕ f D the disrete form is simply the squre of the di'erene etween oth funtions omputed in the points elonging to the (nl interfe Γ f F en exmple of the omputtion is shown in pigure PFIF he weight β 2 ssoited to the interfe position is lwys lrger thn β 1 s the position of the interfe dittes not only the shpe itself ut lso the lotion of oth domins Ω 1 nd Ω 2 thus in)uening diretly the temperture distriutionF he third term ontrols the length of the interfeF st is used in ses where instilities use n inrese in interfe lengthD for exmple when ontrolling the wullinsEekerk instility @see ghpter TAF pinllyD the lst term penlizes the ontrol ost nd leds to the grdient eqution in the djoint prolemF Chapter 2. Adjoint-based optimization algorithm elgorithm I summrizes the steps needed to ompute J F e ssume tht eh hoie of the ontrol vrile w leds to unique sttes T (w) nd ϕ(w)F hereforeD the minimiztion prolem @wA reds min w J (T (w), ϕ(w), w) sujet to @pA. @wA Algorithm 1: gost funtionl omputtion input :

J (T, ϕ, w) = β 1 2 ˆΩ T f -T d 2 dx + β 2 2 ˆΓf ϕ f -ϕ d 2 ds + β 3 2 ˆΓf ds + β 4 2 ˆtf 0 ˆ∂Ω |w| 2 ds dt,
uD T f D T d ϕ f D ϕ d D β 1 D β 2 D β 3 D β 4 D t f output: J J = 0D X for x ∈ ∂Ω do for t ∈ [0, t f ] do J += β 4 w(x, t) 2 end end for x ∈ Ω do J += β 1 T f (x) 2 -T d (x) 2 if x ∈ Γ f then J += β 2 ϕ f (x) 2 -ϕ d (x) 2
tore the interfe entroid position in X end end gompute the totl length of the interfe L = arclength(X)

J += β 3 L

Continuous adjoint derivation

Lagrange functional

sn this work we hose the 4di'erentiteEthenEdisretize4 pproh to solve the minE imiztion prolem @wAF e derive the orresponding djoint prolem in ontinE uous wy nd disretize fterwrdsF vet Θ e the djoint temperture nd ψ the djoint level set funtionF sn order to ompute the grdient of J (w)D we introdue the following vgrnge funtionl L @dxD dt nd ds re omitted for revityA

L (T, ϕ, w, Θ, Θ D , Θ I , ψ) = J (T, ϕ, w) - ˆtf 0 ˆΩ1 ∂T 1 ∂t -∇ 2 T 1 Θ 1 - ˆtf 0 ˆΩ2 ∂T 2 ∂t -∇ 2 T 2 Θ 2 - ˆtf 0 ˆ∂Ω (T -w) Θ D - ˆtf 0 ˆΓ (T -(T M -ε V v S -ε κ κ)) Θ I - ˆtf 0 ˆΓ ∂ϕ ∂t + [∇T i ] 1 2 • ∇ϕ ψ.
@PFIA imilrly to the temperture (eldD Θ i : (x, t) → R + denotes the djoint temperture (eld in either phse @i = 1, 2AF he vgrnge multipliers Θ D nd Θ I re used for the oundry onditions on the domin nd the interfe respetivelyF he djoint system is otined y setting to zero the derivtives of the vgrnge funtionl with respet to

T nd ϕ X L T (•) = L ϕ (•) = 0F etting the initil onditions to T (x, 0) = T 0 (x),
ϕ(x, 0) = ϕ 0 (x), @PFPA leds to δT (x, 0) = δϕ(x, 0) = 0 when lulting the derivtives in the diretion δT nd δϕF sn the following setionsD we detil the proedure of extrting the djoint tefn prolem @eA using the vgrnge funtion L @iqution PFIAF 2.2.2 Adjoint temperature e strt y deriving the djoint temperture equtions y setting L T (•) = 0

L T δT = -ˆΩ β 1 T f -T d δT f - ˆtf 0 ˆΩ1 ∂δT 1 ∂t -∇ 2 δT 1 Θ 1 - ˆtf 0 ˆΩ2 ∂δT 2 ∂t -∇ 2 δT 2 Θ 2 - ˆtf 0 ˆ∂Ω δT Θ D - ˆtf 0 ˆΓ δT Θ I - ˆtf 0 ˆΓ ψ [∇δT i ] 1 2 • ∇ϕ.

20

Chapter 2. Adjoint-based optimization algorithm es expetedD the terms in the vgrngin tht do not depend on the temperture vnishF e now move the sptil nd temporl derivtives towrds the djoint stte ΘF e pply integrtion y prtsD one with respet to time @using the eynolds trnsport theorem eA nd twie with respet to spe @using qreen9s formulA

L T δT = -ˆΩ β 1 T f -T d δT f - ˆΩf 1 δT f Θ f 1 + ˆΩ0 1 δT 0 Θ 0 1 - ˆtf 0 ˆΩ1 ∂Θ 1 ∂t δT + ˆtf 0 ˆ∂Ω 1 δT Θ 1 v S1 + ˆtf 0 ˆ∂Ω 1 ∂δT ∂n Θ 1 - ˆtf 0 ˆΩ1 ∇δT ∇Θ 1 - ˆΩf 2 δT f Θ f 2 + ˆΩ0 2 δT 0 Θ 0 2 - ˆtf 0 ˆΩ2 ∂Θ 2 ∂t δT + ˆtf 0 ˆ∂Ω 2 δT Θ 2 v S2 + ˆtf 0 ˆ∂Ω 2 ∂δT ∂n Θ 2 - ˆtf 0 ˆΩ2 ∇δT ∇Θ 2 - ˆtf 0 ˆ∂Ω δT Θ D - ˆtf 0 ˆΓ δT Θ I - ˆtf 0 ˆΓ ψ [∇δT i ] 1 2 • ∇ϕ ⇐⇒ L T δT = -ˆΩ β 1 T f -T d δT f - ˆΩf 1 δT f Θ f 1 + ˆΩ0 1 δT 0 Θ 0 1 - ˆtf 0 ˆΩ1 ∂Θ 1 ∂t δT + ˆtf 0 ˆ∂Ω 1 δT Θ 1 v S1 + ˆtf 0 ˆ∂Ω 1 ∂δT ∂n Θ 1 - ˆtf 0 ˆ∂Ω 1 ∂Θ 1 ∂n δT + ˆtf 0 ˆΩ1 δT ∇ 2 Θ 1 - ˆΩf 2 δT f Θ f 2 + ˆΩ0 2 δT 0 Θ 0 2 - ˆtf 0 ˆΩ2 ∂Θ 2 ∂t δT + ˆtf 0 ˆ∂Ω 2 δT Θ 2 v S2 + ˆtf 0 ˆ∂Ω 2 ∂δT ∂n Θ 2 - ˆtf 0 ˆ∂Ω 2 ∂Θ 2 ∂n δT + ˆtf 0 ˆΩ2 δT ∇ 2 Θ 2 - ˆtf 0 ˆ∂Ω δT Θ D - ˆtf 0 ˆΓ δT Θ I - ˆtf 0 ˆΓ ψ [∇δT i ] 1 2 • ∇ϕ, with v S1 nd v S2 the veloities of the ontrol volumes Ω 1 nd Ω 2 respetively with v S1 = ⃗ v S • n nd v S2 = -⃗ v S
• nF hese two terms re only nonEzero on ΓF orting the terms y their domin of integrtion nd setting δT 0 = 0D we otin

L T δT = ˆtf 0 ˆΩ1 ∂Θ 1 ∂t + ∇ 2 Θ 1 δT + ˆtf 0 ˆΩ2 ∂Θ 2 ∂t + ∇ 2 Θ 2 δT - ˆΩ1 Θ f 1 -β 1 T f -T d δT f - ˆΩ2 Θ f 2 -β 1 T f -T d δT f + ˆtf 0 ˆ∂Ω - ∂Θ ∂n δT + ∂δT ∂n (Θ -Θ D ) + ˆtf 0 ˆΓ Θ 1 ⃗ v S • n -Θ 2 ⃗ v S • n - ∂Θ 1 ∂n + ∂Θ 2 ∂n -Θ I δT + ˆtf 0 ˆΓ (Θ -ψ|∇ϕ|) [∇δT ] 1 2 • n ∀δT = 0,
where the seond to lst termD orresponds to the norml jump in djoint temperE tureD whih gives

Θ I = -[∇Θ] 1 2 • n euse Θ 1 = Θ 2
on ΓF fy ltering the diretions of vritions nd eliminting ertin termsD we otin the djoint temperture proE lem

                                               - ∂Θ 1 ∂t = ∇ 2 Θ 1 in Ω 1 (t) - ∂Θ 2 ∂t = ∇ 2 Θ 2 in Ω 2 (t) Θ(x, t f ) = β 1 (T f -T d ) in Ω ∂Θ(x, t) ∂n = 0 on ∂Ω Θ(x, t) = ψ|∇ϕ| on Γ(t) Θ I = -[∇Θ] 1 2 • n on Γ(t)
@eA he (rst two equtions of the djoint temperture prolem @eA re the het equE tions in reverse timeF he third eqution is the initil ondition for the djoint temperture (eldD tht depends on the desired temperture distriution

L ϕ δϕ = D J (T, ϕ, w); δϕ +D ˆtf 0 ˆΩ ∂T ∂t -∇ 2 T Θ; δϕ -D ˆtf 0 ˆ∂Ω ∂T ∂n -u Θ D ; δϕ -D ˆtf 0 ˆΓ (T -T M -ε V v S -ε κ κ) Θ I ; δϕ -D ˆtf 0 ˆΓ ∂ϕ ∂t + [∇T i ] 1 2 • ∇ϕ ψ; δϕ .
e divide the ontriutions of the djoint level set termEyEtermF he ontriution from the ost funtionlD removing the terms tht do not depend on ϕD simpli(es to

D J (T, ϕ, w); δϕ = D β 2 2 ˆΓf |ϕ f | 2 ; δϕ + D β 3 2 ˆΓf 1; δϕ .
fy using the theorem on the derivtive of oundry integrl @heorem fAD we otin

D J (T, ϕ, w); δϕ = - β 2 2 ˆΓf δϕ f |∇ϕ| ∂|ϕ f | 2 ∂n + κ|ϕ f | 2 - β 3 2 ˆΓf δϕ f |∇ϕ| κ = - β 2 2 ˆΓf δϕ f |∇ϕ| ∂|ϕ f | 2 ∂n + κ |ϕ t f | 2 + β 3 β 2 .
fy speifying tht δϕ = 0 on δΩD the seond nd third terms re then equl to zeroF por the fourth termD we need to ssume tht the djoint level set ϕ nd the multiplier Θ I re de(ned on ll of ΩF nder suh ssumptionsD we n pply the sme oundry integrl theorem @heorem fA

-D ˆtf 0 ˆΓ (T -T M -ε V v S -ε κ κ) Θ I ; δϕ = - ˆtf 0 ˆΓ δϕ |∇ϕ| ∂ (T -T M -ε V v S -ε κ κ) ∂n + κ (T -T M -ε V v S -ε κ κ) Θ I .
he vlue of the temperture on Γ is lwys equl to T -T M -ε V v S -ε κ κF woreoverD using the ft tht T M is onstnt long the interfeD the ove eqution simpli(es to

-D ˆtf 0 ˆΓ (T -T M -ε V v S -ε κ κ) Θ I ; δϕ = - ˆtf 0 ˆΓ δϕ |∇ϕ| ∂ (T -ε V v S -ε κ κ) ∂n Θ I = - ˆtf 0 ˆΓ δϕ |∇ϕ| ∂T ∂n -ε V ∂v S ∂n -ε κ ∂κ ∂n Θ I .
vooking t the resulting terms in the integrlD we hve ∂T /∂n orresponding to the norml derivtive of the tempertureD whih is well de(ned in ll ΩD ∂v S /∂n tht is onstnt s the veloity of the interfe will e extended nuE merilly in the norml diretion y proedure de(ned in ghpter RD ∂κ/∂nD the norml derivtive of the urvture exists in the level set frmework nd n e omputed in ll of Ω ut hs no true 4physil4 mening hene we ignore this term when dveting the vrilesD leding to n inomplete djoint derivtionF xeverthelessD the optimiztion results otined in ghpter TD in ses where ε κ is nonEzeroD on(rm tht this ssumption is eptleF pinllyD y pplying the sme oundry integrl theorem together with the hin ruleD the lst term simpli(es to

-D ˆtf 0 ˆΓ ∂ϕ ∂t + [∇T i ] 1 2 • ∇ϕ ψ; δϕ = - ˆtf 0 ˆΓ δϕ |∇ϕ| ∂ ∂n ∂ϕ ∂t + [∇T i ] 2 1 • ∇ϕ ψ + κ ∂ϕ ∂t + [∇T i ] 2 1 • ∇ϕ ψ + ˆtf 0 ˆΓ ∂δϕ ∂t + [∇T i ] 2 1 • ∇δϕ ψ = - ˆtf 0 ˆΓ δϕ |∇ϕ| ∂ ∂n ∂ϕ ∂t + [∇T i ] 2 1 • ∇ϕ ψ + ˆtf 0 ˆΓ ∂δϕ ∂t + [∇T i ] 2 1 • ∇δϕ ψ.
fy using the ssumption of onstnt veloity v S extended in the norml diretion on ll of ΩD we otin tht

-D ˆtf 0 ˆΓ ∂ϕ ∂t + [∇T i ] 1 2 • ∇ϕ ψ; δϕ = - ˆtf 0 ˆΓ δϕ |∇ϕ| ∂δϕ ∂t + [∇T i ] 2 1 • ∇δϕ ψ.
sn order to otin the djoint tefn onditionD we need to move the derivtives from δϕ to ψF essuming tht ψ is de(ned on ll of ΩD we n use the orollry g Chapter 2. Adjoint-based optimization algorithm on integrtion y prts in time on moving surfe

- ˆtf 0 ˆΓ δϕ |∇ϕ| ∂ϕ ∂t + [∇T i ] 2 1 • ∇ϕ ψ = - ˆΓf δϕ f ψ f + ˆΓ0 δϕ 0 ψ 0 - ˆtf 0 ˆΓ δϕ ∂ψ ∂t + ∇(δϕ ψ) • v S + δϕ ψ div Γ v S + ψ[∇T i ] 2 1 • ∇δϕ = - ˆΓf δϕ f ψ f + ˆΓ0 δϕ 0 ψ 0 - ˆtf 0 ˆΓ δϕ ∂ψ ∂t - ˆtf 0 ˆΓ δϕ (∇(ψ) • v S + ψ div Γ v S ) + ˆtf 0 ˆΓ ψ[∇T i ] 2 1 • ∇δϕ -ψ∇δϕ • v S .
he lst term nels due to the ft tht [∇T i ] 2 1 = v S on ΓF he resulting terms of the djoint level setD using δϕ 0 = 0D give

L ϕ δϕ = - β 2 2 ˆΓf δϕ f |∇ϕ| ∂|ϕ f | 2 ∂n + κ |ϕ f | 2 + β 3 β 2 - ˆtf 0 ˆΓ δϕ |∇ϕ| ∂T ∂n -ε V ∂v S ∂n Θ I - ˆΓf δϕ f ψ f - ˆtf 0 ˆΓ δϕ ∂ψ ∂t - ˆtf 0 ˆΓ δϕ (∇(ψ) • v S + ψ div Γ v S ) .
fy regrouping the terms y their domins of integrtion nd requiring L ϕ δϕ = 0 ∀ δϕD we otin tht

ψ f = - 1 |∇ϕ| β 2 2 ∂|ϕ f | 2 ∂n + κ |ϕ f | 2 + β 3 β 2 - ∂ψ ∂t = ∇ψ • v S + ψ div Γ v S + 1 |∇ϕ| ∂T ∂n -ε V ∂v S ∂n Θ I .
elling tht Θ I ws de(ned in the djoint temperture nd tht v S is onstnt in the norml diretionD we hve

ψ f = - 1 |∇ϕ| β 2 2 ∂|ϕ f | 2 ∂n + κ |ϕ f | 2 + β 3 β 2 ∂ψ ∂t + div(ψv S ) = 1 |∇ϕ| ∂T ∂n -ε V ∂v S ∂n [∇Θ]
L w δw = β 4 2 ˆtf 0 ˆ∂Ω |w| 2 ; δw - ˆtf 0 ˆ∂Ω ∂T ∂n -w Θ D ; δw = ˆtf 0 ˆ∂Ω (β 4 |w| -Θ) δw,
thus otining t the grdient eqution 0 = β 4 w + Θ on ∂ΩF rereD we used the previously de(ned multiplier Θ D tht is identilly equl to Θ on ∂Ω @see djoint temperture derivtion in etion PFPFPAF e now gther the djoint temperture nd level set with the grdient eqution to onstrut the djoint tefn prolem @eA elowF

pind funtion Θ : Ω × [t f , 0] → R nd funtion ψ X Ω × [t f , 0] → R suh tht                                                                    - ∂Θ 1 ∂t = ∇ 2 Θ 1 in Ω 1 (t) @eFA - ∂Θ 2 ∂t = ∇ 2 Θ 2 in Ω 2 (t) @eFA Θ(x, t f ) = β 1 (T (t f ) -T t f ) in Ω @eFA ∂Θ(x, t) ∂n = 0 on ∂Ω @eFdA Θ(x, t) = ψ|∇ϕ| on Γ(t) @eFeA ∂ψ ∂t + div(ψ⃗ v S ) = 1 |∇ϕ| ∂T ∂n -ε V ∂v S ∂n [∇Θ i ] 1 2 • n on Γ(t) @eFfA ψ(x, t f ) = - β 2 2 ∂ ∂n |ϕ t f | 2 + κ |ϕ t f | 2 + β 3 β 2 in Ω @eFgA 0 = β 4 w + Θ on ∂Ω @eFhA @eA
iqutions eF nd eF orrespond to the het equtions in reverse timeF iqution eF is used to initilize the djoint temperture (eld in oth phsesF sn iqution eFeD the inEhomogeneous xeumnn oundry ondition imposed on ∂Ω existing in the forwrd prolem @pA is mpped to n homogeneous lterntiveF he hirihlet oundry ondition t the interfe is de(ned in iqution eFf where the djoint temperture is equl to the djoint level set vlue ψ ugE mented y the rmiltonin of the level set funtion |ϕ|F Chapter 2. Adjoint-based optimization algorithm sn iqution eFfD the djoint level set dvetion eqution is de(nedF he soure term orresponds to the djoint tefn ondition nd tkes into ount the norml jump in grdient of djoint temperture ross the interfeF he djoint level set funtion is initilized using iqution eFgF he lst iqution @eFhA represents the optimlity onditionD where the leftE hnd side represents the grdient used to updte the ontrol vrile wF 2.3 Gradient-based optimization procedure sn optimiztion prolems where the informtion on the grdient n e extrted through the djointD the onvergene is strongly improved nd this will e shown when ompring our djointEsed proedure with derivtiveEfree oneF sn the fmE ily of the djointEsed methodsD the grdientEsed ones re the methods of hoie to prove tht the optiml ontrol pproh is resonle due to their stright forwrd implementtionF hereforeD we hose to solve the minimiztion prolem @wA y using the limited memory fpq @vEfpqA methodD qusiExewton method origE inlly desried in viuF 8 xoedl IWVWF he min hrteristi of this method is tht it determines the desent diretion y preonditioning the grdient with n pproximtion of the ressin mtrixF his informtion is otined using pst pE proximtions @the numer of pproximtions is determined y the memory length prmeter whih is set to m = 10 in ll our optimiztion test sesA s well s the grdientF es n initil guess for the initil ressinD we use the sled identity mtrix s desried in right 8 xoedl PHHTF elgorithm P summrizes the vE fpq method used in our numeril exmplesF he lgorithm is stopped t given itertion n if one of the following riteri is ful(lled he reltive di'erene in ontrol vrile w n-1 -w n w n < 10 -8 F he reltive di'erene in ost funtionl J n-1 -J n J n < 10 -8 F his riterion n e relxed to llow temporry inrese of the ost funtionlD for exmple to espe9 lol minimumF he norm of the grdient |∇J n | < 10 -6 F sn the next ghpterD we will present novel gut gell method to solve di'usive trnsportF he het equtions in oth the forwrd nd djoint prolem will e solved y this method tht modi(es the stndrd (nite di'erene formul to tke into ount the oundry onditions of given interfe in the underlying grtesin gridF 2.3. Gradient-based optimization procedure 27 Algorithm 2: yptimiztion proedure using the vEfpq method input :

w 0 D m = 10 output: wD T D ϕD ΘD ψ k ← 0D l ← 0
while not converged do olve the forwrd tefn prolem @pA for T k nd ϕ k olve the djoint tefn prolem @eA for Θ k nd ψ k gompute the grdientX

∇J k = β 4 w k + Θ k if k ≥ 1 then s k-1 = σ k-1 d k-1 g k-1 = ∇J k -∇J k-1 if (s k-1 ) T g k-1 ≤ 0 then l ← 0 else if (s k-1 ) T g k-1 > 0 then l → l + 1 if l > m then
emove {s l-m , g l-m } end edd {s l-m , s l-m } end end ghoose n initil pproximtion to the inverse of the ressin H k 0 gonstrut the diretion d k = -H k ∇J k hetermine σ k using vine erh lgorithm with ktrking where utEell method on stggered gridsF he ommon trdeEo' in utEell methods is etween expliit hndling of mny speil ses in whih omputtionl ell n e interseted y ojet oundries nd mking use of lkEox lestEsqures methods to onstrut derivtivesF he seond method n e shown to lk the desired onservtion propertiesF he level set funtionD de(ned in ghpter ID is used there to ompute geoE metri momentsD suh s the entroid oordintes of wet volumes or the re of wet fesF yne dvntge of the proposed method is tht these geometri (elds re the only informtion required to modify lssil (nite di'erenes formuls in the viinity of oundriesD while reovering wellEknown seondEorder formuls wy from them nd ommodting ny strethingF he de(nition of these geE ometri (elds nd their numer is determined from ury onsidertionsF st will e shownD in prtiulrD tht the proposed opertors degenerte to lssiE l formuls for the meshEligned oundriesF he disrete lulus of worinE ishi worinishi et al. IWWVD worinishi PHIH is leverged to provide onise expresE sions for the disrete opertorsD whether for hirihlet or xeumnn oundry ondiE tionsF sn this workD we will only fous on the hirihlet seF lidtion re provided tht ssess the sheme9s seond order uryD stility nd ility to urtely represent trnsfersF 3.1.2 Governing principles e suset of the opertors de(ned y worinishi et al. IWWV is reemployedD nmely the seond order di'erentitionD interpoltionD nd permnent produt opertors on unitry grids de(ned respetively s

σ k = argmin J (w k + σ k d k ) pdte w k+1 = w k + σ k d k k → k +
δϕ δξ x i+ 1 /2,j = ϕ i+1,j -ϕ ij , @QFPA ϕ x i+ 1 /2,j = ϕ i+1,j + ϕ ij 2 , @QFQA nd ϕψ x i+ 1 /2,j = ϕ i+1,j ψ ij + ψ i+1,j ϕ ij 2 .
@QFRA hese de(nitions re strightforwrd to extend to either stggered (ϕ i+ 1 /2,j ) nd

(ϕ i,j+ 1 /2 ) or nodl (Φ i+ 1 /2,j+ 1 /2 ) (eldsF vikewiseD opertors ting upon the seond diretionD δ(•)/δξ y D (•) y nd (•, •)
y D re de(ned in the sme mnnerF pinllyD the restritions to the oneEdimensionl se is otined y dropping the lst indexF hese opertors will e used in the susequent formul when deriving the disrete opertorsF gut gell wethods re (rmly grounded in the pinite olume wethodD whih de(nes the primry disrete vriles s ellEwise verges over mesh elements @s opposed to pointEwise vlues in the pinite hi'erene wethodD for exmpleAF he design of the pinite olume opertors is then sed on the pplition of tokes9 theoremF por exmpleD given slr (eld T D this theorem sttes tht in grtesin oordinte systemD the x omponent of the grdient q ≡ ∇T verged over ell Ω my e omputed s

|Ω| q x = ˆΩ ∂T ∂x dV = ˛∂Ω T e x • dS @QFSA
where |•| denotes the mesure opertorD dS the outwrdEpointing surfe elementD e x the unit vetor long the x diretion nd (∂•) the ontour opertorF por the ske of presenttionD the se displyed in pigure QFI is onsideredD where Ω onsists of the intersetion of phse domin nd omputtionl ell @ right hexhedronAF he ontour ∂Ω then onsists of the union of the three plnr fes A - x D A - x nd A - y s well s the oundry surfe ΓF e pieeEwise liner pproximtion of ΓD denoted ΓD of length Γ nd unit norml (n x , n y )D n e de(ned s done in pigure QFIF epplying iqution QFS to Ω with T = 1 then yields whih highlights the existene of fundmentl reltionD sometimes referred to s urfe gonservtion vw @gvA

ˆ Ω ∂1 ∂x dV = A + x -A - x + n x Γ = 0,
A + x -A - x = -n x Γ . @QFTA
sn other wordsD the knowledge of (|A α |) α∈{x,y,z} impliitly de(nes pieeEwise liner pproximtion to the oundryF es onsequeneD this surfe informtionD heneforth referred to s the surfe pitiesD my serve to pproximte the rightE hnd side of iqution QFSF sf the unknowns T ±

x/y re de(ned s verges over the wet res A ±

x/y D the formul

˛∂ Ω T e x • dS = A + x T + -A - x T --A + x -A - x D,
is ext provided D is the hirihlet ondition verged over the pproximte oundE ry ΓF o omplete the de(nition of the verged xEomponent of the grdientD the volume pity V ≡ |Ω| is lso requiredD whih results in the following tenttive grdient opertor

q v1 x ≃ A + x T + -A - x T --A + x -A - x D /V.
st is worth stressing tht the use of the gv @iqution QFTA in q v1

x gurntees tht the disrete grdient vnishes when the solution nd oundry vlues re mthing onstnts @T + = T -= DAF his nottion n e generlized to ritrry dimensions for ny oundry geometry using the di'erentition opertor δ • /δξ α D α ∈ {x, y} @de(ned in iqution QFPAD s follows

∀α ∈ {x, y} , grad v1 α (T α , D) = 1 V δA α T α δξ α - δA α δξ α D , @QFUA
where ll omponents of the disrete vetor (eld Q = (q α )D α ∈ {x, y, z}D re olloE ted with DF sn this (rst version of the grdient opertorD the primry unknowns T α re olloted with the surfe res A α D wheres the hirihlet oundry ondiE tion D is stggered inEetweenF his onstrutionD referred to s erkw iEgrid erkw 8 vm IWUUD relies on the de(nition of multiple temperture (eldsF uh grid on(gurtion is not the one dopted y the weg pproh rrlow 8 elh IWTSD whih fvors the gEgrid tht de(nes single temperture (eld olloted with the D (eld hereF e gEgrid however mens tht the temperture unknowns T nd surfe pities A α re stggeredD in whih se the ltter together with V should e interpolted s follows

∀α ∈ {x, y} , grad v2 α (T, D) = 1 V α δA α T δξ α - δA α δξ α D α , @QFVA
whih uses the interpoltion opertor • α D α ∈ {x, y, z} @de(ned in iqution QFQAF

Loss of accuracy with interpolation

st should e pointed out tht formuls other thn iqution QFV n lso e writE ten without interpoltion of the geometri pitiesD for exmple y olloting ll surfe pities (A α ) with the primry vrile T F rowever in the ontext of seond order opertor suh s the slr vplin @iqution QFIAD the need for interpoltion will resurfe in the pproximtion of the divergene opertorF his setion therefore fouses on the limittions of the seond tenttive formul @iqution QFVAD more spei(lly its filure to revert to lssil (rst order pproxE imtion of the seond order derivtive in the limit where the oundry is orthogonl to the diretion of interestF his is the entrl point of the proposed utEell methodD nmely the enhnement of the geometri desription of the oundry y mens of dditionl volume nd surfe pities to revise the grdient nd divergene operE tors so s to hieve (rstEorder ury in the viinity of meshEligned oundriesF o illustrte the limittion of the tenttive grdient formul @iqution QFVAD the disretiztion of the seondEorder derivtive long x in the meshEligned twoE dimensionl on(gurtion displyed in pigure QFP is onsideredD where the )uid oupies the rightmost ellsF his on(gurtion is hrterized y

V 0 = 0D V 1 = (h x -g)h y D V 2 = h x h y D A 1/2 = 0 nd A 3/2 = A 5/2 = h y @hereD A stnds for A x sine only the x ontriution is onsideredAF sing these expressionsD iqution QFV simpli(es to q -1/2 = 0D q 1/2 = T 1 -D 1 g , @QFWA nd q 3/2 = T 2 -(T 1 + D 1 ) /2 (g + h x ) /2 . @QFIHA Figure 3
.2: Insucient geometric information resulting in loss of accuracy in mesh-aligned geometries.

his pproximtion of the grdient is prolemti for two resonsF et the oundryD the xEgrdient vlue @q 1/2 A is under predited y ftor of 2D sine the denomintor of the rightEhnd side of iqution QFW stnds t g when it should mth the distne etween the points where D 1 nd T 1 re de(nedD g/2F ewy from the oundryD the xEgrdient vlue @q 3/2 A depends on the oundry ondition D 1 D when one would simply expet the di'erene (T 2 -T 1 ) to pper in the numertor of the rightEhnd side of iqution QFIHF his simple exerise highlights the loss of ury ssoited with the interpoE ltion of the geometri pitiesF his n e ssoited with the ft tht the pities re de(ned s volume nd surfe integrls of the hrteristi funtion of the )uid domin

Ω f ⊂ ΩD de(ned s ∀x ∈ Ω, H f (x) ≡ ˆy∈Ω f δ (x -y) dV, @QFIIA
where Ω denotes the omputtionl domin nd δ the multiEdimensionl hir delt funtionF H f is not di'erentile in the lssil senseD nd one should tred reE fully not to interpolte or di'erentite its surfeE or volumeEverged vluesF 3.1.4 Additional geometric information to restore accuracy en intuitive ide to llevite the interpoltions in iqution QFV is to dd new inforE mtion where the volume @ellEentered nd denoted V A nd surfe @feEentered nd denoted (A α )A pities were previously interpoltedF hese new quntitiesD referred to s seondEkind pitiesD omplement the lredy used (rstEkind piE ties V nd (A α )D lso ome in volume @feEentered nd denoted (W α )A nd surfe @ellEentered nd denoted (B α )A forms for ritrry geometriesF his dditionl quntities yield the (nl grdient formul

∀α ∈ {x, y} , grad α (T, D) = 1 W α δB α T δξ α + δ(A α α -B α )D δξ α - δA α δξ α D α , @QFIPA tht supersedes grad (v1) α nd grad (v2) α F
Figure 3.3: Enhanced geometric information restoring accuracy in mesh-aligned geometries.

o show how the ddition of the seondEkind pity restores (rstEorder ury in the grdient omputtionD the on(gurtion displyed in pigure QFQ is onsideredF ine only x derivtives re onsideredD A gin will stnds for A x D wheres W nd B will respetively stnd for W x nd B x F fering this in mindD the on(gurtion under study is hrterized y

V 0 = 0D V 1 = 2f h y nd V 2 = h x h y D A -1/2 = A 1/2 = 0D A 3/2 = A 5/2 = h y D B 0 = 0D B 1 = B 2 = h y nd (nlly W -1/2 = 0D W 1/2 = f h y D W 3/2 = gh y nd W 5/2 = h x h y F sing these expressionsD iqution QFIP simpli(es to q -1/2 = 0 nd q 1/2 = T 1 -D 1 f , q 3/2 = T 2 -T 1 g nd q 5/2 = T 3 -T 2 h x .
T 0 does not pper sine it is outside of the )uid dominD nd the oundry ondition @DA ppers only in the fes djent to the oundryF he formuls otined from iqution QFIP re lssil formuls sine f D g nd h x re the distnes over whih the di'erenes T 1 -D 1 D T 2 -T 1 nd T 3 -T 2 re de(nedF pinllyD in the )uid domin nd wy from the oundriesD iqution QFIP simply reverts to the lssil grdient

formul q x,i+1/2 = T i+1 -T i x i+1 -x i .
sn ftD the ddition of the seondEkind pities is lso su0ient to de(ne the @ellEenteredA volumeEweighted divergene opertorD whih onsists of the sum of the following ontriutions ∀α ∈ {x, y} where N α denotes the oundry vlue of q α F sf one sets the divergene to the produt of the volume V with the lol vlue of the soure term σ s in the originl oisson prolem @iqution QFIAD the on(gurtion displyed in pigure QFQ yields the trivil eqution 0 = 0 in the (rst ellD nd

h y q 3/2 -N 1/2 = 2f h y σ 1 , h y q 5/2 -q 3/2 = h x h y σ 2 ,
in the othersF egin lssil formuls tht degenerte to the lssil formul

h y q i+1/2 -q i-1/2 = h y x i+1/2 -x i-1/2 σ i ,
in the )uid domin wy from the oundryF pinllyD the unknown N = (N α ) n e eliminted y sustituting the grdient formul @iqution QFIPA in the divergene formul de(ned elow

div α q, N = δA α Q α δξ α + δ(B α -A α )N α δξ α - δB α δξ N α α .
@QFIQA he oundry ontriution @the lst two terms in the rightEhnd side of iquE tion QFIQA re set to

α δ(B α -A α )N α δξ α - δB α δξ N α α = α δ(B α -A α )Q α δξ α - δB α δξ Q α α ,
whih mount to identifying the het )ow through the oundry to the norml omE ponent of the temperture grdientF sn the on(gurtion displyed in pigure QFQ this yields only one nonEtrivil equtionD N 1 = q 1 F utting it ll togetherD the proE posed grdient nd divergene opertorsD de(ned for ritrry oundry geometries in iqutions QFIP nd QFIQD disretize the oisson prolem @iqution QFIA in the on(gurtion displyed in pigure QFQ s 0 = 0D

h y T 2 -T 1 g - T 1 -D f = 2f h y σ 1 , nd h y T 3 -T 2 h x - T 2 -T 1 g = h x h y σ 2 ,
in the three ells displyedD while reverting to the lssil formul

h y T i+1 -T i x i+1 -x i - T i -T i-1 x i -x i-1 = h x h y σ i ,
in the )uid domin wy from the oundryF es onsequeneD formuls iquE tions QFIP nd QFIQ n e interpreted s generliztions of the lssil seondEorder formuls to ommodte the presene of ritrry oundries while preserving (rstE order ury in the presene of meshEligned sesF he next setion lri(es the nottion employed thus frD in prtiulr the de(nition of volume nd surfe E pities of the (rst nd seond kindsF 3.2 Volume and surface capacities 3.2.1 Mesh and geometry input es fr s the grtesin mesh is onernedD retiliner mesh with n x × n y ells is de(ned y speifying the following sets of userEde(ned sisss

x1 /2 < x3 /2 < • • • < x nx+ 1 /2 , y1 /2 < y3 /2 < • • • < y ny+ 1 /2 .
smportntlyD the mesh need not e uniformF eny given ell Ω ij D identi(ed y multiEindex ijD (i, j) ∈ 1, n x × 1, n y D orresponds to the set of points (x, y) tht simultneously stisfy x i-1 /2 < x < x i+ 1 /2 nd y j-1 /2 < y < y j+ 1 /2 F egrding the oundry desriptionD there exists wide rnge of tehniques to de(ne )uid doE minD suh s simpliil meshes or gonstrutive olid qeometry @gqA primitives nd opertionsF smpliit representtions y mens of void frtion or distne funtion @vevel etA re lso ommonly usedF egrdless of the method employedD the ssemly of the ut ell opertors requires the omputtion of res nd volE umes tht orrespond to the intersetion of the )uid domin with grtesin elements @fes or ellsAD s displyed in pigure QFRF sn this workD these omputtions re perE formed using wrhing qures lgorithmD tht only requires the level set funtion vlues t the omputtionl ells @detiled in the following setionAF pollowing the omputtion of the pitiesD the geometry input is disrdedF 3.2.2 Capacities of the rst kind gonsider the grtesin mesh displyed in pigure QFRD prtitioned into )uid @Ω f A nd solid @Ω s A domins seprted y oundry @ΓAF sn (nite volume settingD the primry vriles Φ ij onsist of verges of ny given ontinuous (eld (x, y) → ϕ (x, y) over the intersetion of the )uid domin with ny given hexhedrl ellD de(ned s follows

V ij ≡ ˆΩij ϕ (x) H f (x) d 2 x,
ϕ ij V ij ≡ ˆΩij ϕ (x) H f (x) d 2 x,
where H f is the )uid hrteristi funtion de(ned in iqution QFIIF he set V ≡ (V ij ) is referred to s the volume pities of the (rst kindF hen the (eld under onsidertion is linerD these verges oinide with the vlues t the )uid enter of mssD displyed in pigure QFSD de(ned for ells fully or prtilly oupied y the )uid @referred to s mixed ellsAF elthough it does not pper expliitly in the gut gell opertorsD the oordintes of the )uid enter of mss @displyed with rosses in pigure QFSA re still required to de(ne the seond kind pitiesD nd re therefore temporrily storedF hey re denoted s X nd Y nd de(ned for ny ell

Ω ij s X ij Y ij V ij ≡ ˆΩij x y H f (x) d 2 x.
he seond step onsists in omputing the re of eh of the fes wet y the )uidF feuse the mesh is grtesinD the fes djent to eh ell re lelled sed on the diretion they re orthogonl toF hese quntitiesD referred to s surfe pitiesD re stggered nd denoted (A α ) @α ∈ {x, y}AF hey re de(ned s

A x i+ 1 /2,j ≡ ˆyj+1/2 y j-1/2 H f x i+ 1 /2 , y dy, nd A y i,j+ 1 /2 ≡ ˆxi+1/2 x i-1/2 H f x, y j+ 1 /2 dx.

Volume and surface capacities 39

Figure 3.5: First kind capacities: V (lled areas), A 1 (dashed vertical lines), A 2 (dashed horizontal lines) and X and Y (crosses).

Capacities of the second kind

he oordintes of the )uid enter of mss re used s followsF por eh diretionD the volume informtion is enrihed y mesuring how muh )uid lies etween eh enter of mssF his yields s mny sets of stggered volumesD nmely

W x i+ 1 /2,j ≡ ˆyj+1/2 y j-1/2 ˆXi+1,j X ij H f (x) d 2 x, nd W y i,j+ 1 /2 ≡ ˆYi,j+1 Y ij ˆxi+1/2 x i-1/2 H f (x) d 2 x,
tht will e denoted (W α )D α ∈ {x, y}D nd referred to s volume pities of the seond kindF he pities W x nd W y re represented s olored res in the on(gurtion displyed in pigure QFT nd QFTD respetivelyF vikewiseD the re wet y the )uid for the meshEligned fes tht interept the )uid enter of mss will e required in eh ellF his yields n dditionl set of ellEentered quntitiesD

B x ij = ˆyj+1/2 y j-1/2 H f (X ij , y) dy, nd B y ij = ˆxi+1/2 x i-1/2 H f (x, Y ij ) dx,
referred to surfe pities of the seond kindF he pities B x nd B y re represented s olored dshed lines in the on(gurtions displyed in pigure QFT nd QFTD respetivelyF 40 Chapter 3. A novel Cut Cell method for diusive transport (a) W1 (lled areas), B1 (dashed vertical lines) and X and Y (crosses).

(b) W2 (lled areas), B2 (dashed horizontal lines) and X and Y (crosses).

Figure 3.6: Schematic of the second kind capacities.

Marching squares algorithm

he omputtion of res nd volumes is rried out y using mrhing squres lgorithm using the vlues of the level setF his lgorithm is used to determine one of the IT di'erent possile ses in two dimension @pigure QFUAD thus llowing the omputtion of the rossing points etween the interfe nd the edgesF es preliminry stepD we wnt to lote the mixed ells in our dominF e look t vlues of the level set funtion whih re of di'erent sign thn the neighors in eh diretionF ith this onditionD we ensure tht we re lose9 to the interfeF xote tht this is not su0ient to determine if the ell is mixed or not s one ould tg the ells djent to the mixed onesF yne this (rst suEset of ells is knownD we n pply the mrhing squres lgorithmF he lgorithm works s follows IF he vlues of the orners of ellD hrterized y these R rdinl positions @southEwestD southEestD northEestD northEwestAD re omputed using iE qudrti interpoltionF PF qiven the signs of the ornersD the unique isovlue of the ell is omputed y the formul iso = SW +2 SE +4 NE +8 NW . @QFIRA QF hepending on the isovlue @seAD the intersetion points re loted reEusing the iEqudrti interpoltionF o determine the vlues of the ornersD we perform iEqudrti interpoltion on the 3 × 3 stenil entered on the ell of interestF por this exeriseD we ssume onstnt sping of the grtesin grid in ll dimensionsF e wnt to lulte the interpoltion mtrix M given y

ϕ(x, y) = 2 i=0 2 j=0 m ij x i y j = x 2 x 1   m 2,2 m 2,1 m 2,0 m 1,2 m 1,1 m 1,0 m 0,2 m 0,1 m 0,0     y 2 y 1   = XM Y T ,
with ϕ(x, y) the set of known dt points of the level set funtionD yielding

Φ = GM G T   ϕ 1, 1 ϕ 1,0 ϕ 1,1 ϕ 0, 1 ϕ 0,0 ϕ 0,1 ϕ 1, 1 ϕ 1,0 ϕ 1,1   Φ =   ( 1) 2 1 1 0 2 0 1 1 2 1 1   G   m 2,2 m 2,1 m 2,0 m 1,2 m 1,1 m 1,0 m 0,2 m 0,1 m 0,0   M   ( 1) 2 (1) 2 1 1 1 1   G T .
Chapter 3. A novel Cut Cell method for diusive transport he interpoltion mtrix M n now diretly e omputed y

M = G 1 Φ G T 1 = G 1 Φ G 1 T ,
with the mtrix G 1 de(ned s

G 1 =   0.5 1 0.5 0.5 0 0.5 0 1 0   .
yne the interpoltion mtrix is determinedD we hve the interpolted vlues of the level set funtion t the orner nodes ϕ i±1/2,j±1/2 @ieF t the points @x = ±0.5D y = ±0.5AAF he isovlue orresponding to the ell is then determined y the reltion desried in iqution QFIRF xote tht ses S nd IH do not existD s we ssume tht the interfe is rossing the underlying grtesin grid only one per ellF he sme interpoltion mtrix is then used to (nd the lotion of the intersetion point y (nding the zeros of the iEqudrti interpoltion long the edges of interestF ith the rossing pointsD we n determine the entroid of the wetted re y using the following formulF qiven polygon de(ned y n verties

(x 0 , y 0 ), (x 1 , y 1 ), ..., (x n-1 , y n-1 )D the point C = (C x , C y )
is the entroid where

C x = 1 6A n-1 i=0 (x i + x i+1 )(x i y i+1 -x i+1 y i ), C y = 1 6A n-1 i=0 (y i + y i+1 )(x i y i+1 -x i+1 y i ),
nd A is the signed re de(ned y the shoele formul

A = 1 2 n-1 i=0 (x i y i+1 -x i+1 y i ).
he lotion of the entroid of the ell will not only e used in the omputtion of the seond kind pities ut lso in the initiliztion of the temperture (eld in the fresh ells when the interfe unovers new ells s the front moves @see ghpter RAF woreoverD we de(ne the interfe entroid s the mid point of the segment rossing the ell whih will e used in the omputtion of the tefn onditionF elgorithm Q summrizes the steps presented in this setionF he purpose of this lgorithm is twoEfolded X @iA it onstruts the disrete vplE in opertors using the fes nd volume pitiesD nd @iiA it lotes the interfe in our grtesin grid s well s the ell nd interfe entroidsY informtions whih will e used in most of the numeril steps of solving the tefn prolemF Algorithm 3: gut gell method pseudoEode input : ϕD N @size of the dominA output: L @vplin opertorA for (i, j)

∈ (2 : N 1, 2 : N 1) do if ϕ i,j ϕ i±1,j±1 < 0 then gonstrut the 3 × 3 stenil Φ gompute the interpoltion mtrix M = G 1 Φ G 1 T
hetermine the sign of the verties loted t ( 0.5, 0.5) SW , (0.5, 0.5) SE , (0.5, 0.5) NE , ( 0.5, 0.5) NW gompute the isovlue iso = SW +2 SE +4 NE +8 NW if iso ̸ = 0 and iso ̸ = 15 then tore the indies (i, j) of the mixed ell vote the intersetion points depending on the se gompute nd store the ell entroid on oth sides of the interfe gompute nd store the interfe entroid gompute the fes nd volume pities of (rst kind

A 1 , A 2 , V nd seond kind B 1 , B 2 , W 1 , W 2 end end end
gonstrut the vplin opertor L

Validation on stationary geometries

ith the disrete vplin opertors @one for eh phseA now onstrutedD we n solve the het equtions on oth domins

∂T ∂t

= LT, @QFISA on oth dominsF rereD L is the disrete liner opertorF e ouple the gut gell spe disretiztion with grnkExiolson time disretiztionD where τ is the time stepD ∆ the uniform grid sping nd n the urrent itertionD resulting in the following disrete system

T n -T n-1 τ = 1 2 1 ∆ 2 LT n + 1 ∆ 2 LT n-1 , @QFITA
whih requires the solution of liner system forming pentdigonl mtrixF e vlidte the method in di'erent sttionry setupsF e onvergene study is rried out for these sesD where the referene solution is tken s the simultion with the highest numer of points per dimensionF e look t the L 2 norm of the error in mixed ellsD full ells @ells tht re not mixedA nd in ll the ellsF sn eh seD the initil temperture (eld is set to zero nd we impose hirihlet oundry ondition t the interfeF he rtio τ /∆ 2 = 0.5 is kept onstnt s we inrese the numer of pointsF e onsider the following test sesX

IF e solid irle of rdius R = 0.85D initilized in 2 × 2 dominF he level set funtion is de(ned s ϕ(x, y) = x 2 + y 2 -R.
he hirihlet T D = 1 oundry ondition is imposed t the interfeF e solve only for the phse inside of the irle until (nl time t f = 0.03125F he simultion is rried out for di'erent resolutions x a ITD QPD TRD IPV orresponding to RD ITD TRD PST itertionsD respetivelyF he referene solution is tken for x a PSTF he results re summrized in pigure QFVF es expetedD the order of onvergene of the error in full ells is lose to P while the order of onvergene in mixed ells is slightly less thn PF his drop in order in mixed ells is due to the ssumption of pieeEwise liner interfe pproximtion s well s the umultion of errors of the iEqudrti interpoltionF he mximl errors re lolized in ells where the wetted re is smll @typilly smller thn 5% of ∆ 2 AF xeverthelessD the glol order of onvergene in ll ells is extly PF PF imilr to the (rst seD we initilize solid irle of rdius R = 0.75F he domin size re the sme s well s the onsidered grid resolutions thn in the previousF his timeD we solve outside of the irle with hirihlet oundry ondition T D = 1 t the interfe nd insulted oundry onditions t the domin oundryF he resultsD presented in pigure QFW re similr to se I with slight drop in solute errorF his is due to the ft tht there re less points per dimeter thn previouslyF his se vlidtes the implementtion of the xeumnn oundry ondition imposed t the domin oundriesF QF sn this third seD we initilize squre of re 1.6 × 1.6 in 2 × 2 dominF he level set funtion is de(ned s ϕ(x, y) = max((x -0.8), (x + 0.8), (y -0.8), (y + 0.8).

e impose hirihlet oundry ondition T D = 1 nd solve inside of the squre until the sme (nl time t f = 0.03125 with the sme resolutions onE sidered previouslyF sn pigure QFIHD we n see tht the mximl errors re loted t the ornersF he order of onvergene for full ells is similr to the irle ses s well s for ll ellsF his se exhiits the roustness of the method when deling with mesh ligned geometries s explined in etion QFPF RF pinllyD in the lst seD we onsider rystl in 2 × 2 domin where the level set funtion is de(ned s

ϕ(x, y) = x 2 + y 2 -R -0.2 cos (6α) ,
where α is the ngle of the interfe with respet to the x xis nd R = 0.7F

Validation on stationary geometries 45

et the interfeD we impose the qisEhomson reltion

T D = ε κ κ,
with κ the urvture nd ε κ = 0.01F he resulting temperture (eld will now depend on the sign nd mplitude of κF e solve only for the phse inside of the irle until (nl time t f = 0.0078125F he simultion is performed for di'erent resolutions x a QPD TRD IPVD PST orresponding to RD ITD TRD PST itertions respetivelyF he referene solution is tken for x a SIPF sn pigure QFIID we n oserve drop in order of onvergene for full ells with respet to the ses where T D ws onstntF his is explined y the ury of the urvture omputtion @iqution IFVAF he error is mximl in regions where the rdius of urvture is lrgeD where the interfe is qusi ligned with the gridF ith these vlidtion sesD we lose the hpter on the gut gell method for di'usive trnsportF sn the next hpterD we desrie the rest of the numeril steps of the twoEphse tefn prolemF st should e noted tht in ses where the )ow in the liquid phse is onsideredD the lgorithm is modi(ed suh tht the xvierEtokes equtions in the foussinesq pproximtion re solved prior to the onvetionEdi'usion eqution @iqution IFIIAY thus efore the seond stepF he solution of the xvierEtokes equtions in the gut gell frmework re detiled in eppendix fF ih of these steps nd their orresponding onvergene properties re presented in the following setions nd summrized in elgorithm RF hile solving the forwrd prolemD we store the temperture (eldD the interfe position nd the phseEhnge veloity for use in the djoint lgorithm@see elgorithm TAF Algorithm 4: qeneri lgorithm for the solution of the forwrd prolem input : T 0 , T D , w, ϕ 0 , u 0 , P 0 , t f output: es stted previouslyD the motion of the interfe is solely dependent on the jump in the temperture grdientF st is therefore importnt to ompute the norml grdient in temperture of eh phse urtelyD nd to this endD the tohnsenEgolell method tohnsen 8 golell IWWV is used

∇T | Γ = 1 d B -d A d B d A (T D -T * A ) - d A d B (T D -T * B ) , @RFIA
with T D the hirihlet vlue imposed t the interfeD T * A nd T * B the interpolted vlues of the temperture (eld on points A nd B respetivelyD nd d A nd d B the distnes etween the interfe entroid to A nd B respetivelyF he lgorithm for omputing the grdient in one phse is s followsX IF e shifted 3 × 3 stenil is hosenD s shown in pigure RFIF PF e line from the interfe entroid is st in the norml diretion nF QF he rossing points A nd B of this line nd the vertil @or horizontlD deE pending on the norml orienttionA segments of the neighoring Q points re identi(ed nd the distnes

d A nd d B re omputedF RF he vlues T * A nd T * B re interpolted using T 1 A , T 2 A , T 2 A nd T 1 B , T 2 B , T 2 B D reE spetively
yne the norml grdient is omputed in eh phseD using iqutionF RFID the jump is omputed s

[∇T ] 1 2 = ∇T 1 | Γ -∇T 2 | Γ . @RFPA
he disrete veloities of the front in the prtil ells re initilized with this jump nd will e used s the oundry ondition in the veloity extension lgorithmF o vlidte the method within our gut gell frmeworkD we onsider sttionry irle of rdius R = 0.5 in 1 × 1 domin nd we initilize the temperture (eld with similrity solution of the het eqution

T (r) =    T ∞ 1 - f (r) f (R) , r > R 0, r < R @RFQA
with T ∞ = -0.5 given underEooling tempertureD nd

f (r) = E 1 (1/4 r 2 )
where where F is the extended veloity (eld equl to v S t the frontD t ⋆ denotes pseudoE time nd S(ϕ) is the signture funtion

E 1 (t) = ˆ∞ x e -t t dt T 1 A T 1 B T 2 A T 2 B T 3 A T 3 B n T D T * A T * B d B d A Solid Liquid
S(ϕ) =    -1 if ϕ < 0 0 if ϕ = 0 +1 if ϕ > 0 @RFSA
iqution RFR is then disretized using (rst order upwind sheme nd integrted in time y forwrd iuler method until stedy stteF king n s the norml vetor de(ned s n

= (n x , n y ) = ϕ x / (ϕ 2 x + ϕ 2 y ), ϕ y / (ϕ 2 x + ϕ 2 y ) , @RFTA
the disretistion leds to

F n+1 ij =F n ij -τ ⋆ S ij n x ij + F ij -F i-1j ∆ + S ij n x ij -F i+1j -F ij ∆ + S ij n y ij + F ij -F ij-1 ∆ + S ij n y ij -F ij+1 -F ij ∆ @RFUA
where ∆ is the uniform grid spingD (x) + = max(0, x) nd (x) -= min(0, x)D nd the time step τ ⋆ is hosen so tht τ ⋆ /∆ 2 = 0.45F he pseudoEtime spwn in the veloity extension lgorithm is purely (titious nd the numer of itertions in iqution RFU orresponds to the width of the nrrowEnd @xfA round the HE level set where the veloities re initilizedF pigure RFQ shows n exmple of n initil veloity (eld for di'erent nrrowEnd widthsD fter one itertion of the het 56 Chapter 4. Numerical methods for forward and adjoint problems equtionD where the temperture t the interfe T D = ε κ κ depends only on the urvtureF he omputed disrete veloities re positive in the kinks nd negtive in the tips driving the initil rystl towrds irulr shpeF

∂ϕ ∂t + F ∇ϕ |∇ϕ| • ∇ϕ = 0, @RFVA
whih is then divided into onservtive nd nonEonservtive terms

∂ϕ ∂t + ∇ • F ϕ ∇ϕ |∇ϕ| a -ϕ∇ • F ∇ϕ |∇ϕ| b = 0, @RFWA
resulting in seond order prtil di'erentil eqution kin to weighted di'usion equtionF he (rst term @aA hs di'usion oe0ient F ϕ tht depends on the solution nd represents nonliner urvture )ow wheres in the seond term @bA the solution is multiplied y the urvture of its levelEsetsF he min ide ehind this sheme is to distinguish two sesX if the produt F ϕ is negtive @positiveD respetivelyA then a represents forwrd @kwrdD respetivelyA di'usion nd b represents kwrd @forwrdD respetivelyA di'usionF he forwrd di'usion is treted impliitly while the kwrd di'usion is treted expliitly leding to semiEimpliit sheme with di'usive gpv numerF sn order to disretize iqution RFWD we use the sme nottion s

F φpq |∇ϕ pq | (ϕ q -ϕ p ) - q∈N (p) F φp |∇ϕ pq | (ϕ q -ϕ p ) = 0, leding to ˆp ∂ϕ ∂t dx + q∈N (p) F ( φp -φpq ) |∇ϕ pq | (ϕ p -ϕ q ) = 0. @RFIIA
vooking t the term F ( φpφpq )D we n distinguish two sesX sf the term is positiveD we hve forwrd di'usion9 or in)ow towrds the ellF sf the term is negtiveD we hve kwrd di'usion9 or out)ow from the ellF e therefore de(ne the di'usion oe0ient a pq s

a pq = F ( φp -φpq ) |∇ϕ pq | , @RFIPA
nd the relted dominnt forwrd nd kwrd di'usion prts s a f pq = max(a pq , 0), a b pq = min(a pq , 0). @RFIQA sing kwrd iuler time disretiztionD tking the forwrd ontriution exE pliitly nd the kwrd ontriution impliitlyD iqution RFII gives the following 58 Chapter 4. Numerical methods for forward and adjoint problems liner system

ϕ n p + τ ∆ 2 q∈N (p) a f pq (ϕ n p -ϕ n q ) implicit = ϕ n 1 p + τ ∆ 2 q∈N (p) a b pq (ϕ n 1 p -ϕ n 1 q ) explicit @RFIRA
where τ is the time stepD ∆ the uniform grid sping nd n given time stepF he question remining is how to tkle the reonstrution of φp D φpq nd |∇ϕ pq |F he ltter is onstruted using the dimondEell strtegy desried in wikul 8 yhlerger PHIHF gonsider ϕ 1 pq nd ϕ 2 pq the vlues of ϕ loted t the end points x 1 pq nd x 2 pq of the edge e pq s shown in pigure RFRF ϕ 1 pq nd ϕ 2 pq re de(ned s the verge of the R djent vlues of ϕF fering this in mindD we n de(ne the reonstruted rmiltonin s

|∇ϕ pq | = ϕ p -ϕ q ∆ 2 + ϕ 1 pq -ϕ 2 pq ∆ 1/2 @RFISA x p x q x 2 pq x 1 pq n pq a f
S f p = q∈N (p) a f pq , S b p = q∈N (p) a b pq . @RFIUA
QF sf the forwrd di'usion is dominnt @ieF S f p > S b p A then no further steps re needed s the disretiztion will len towrds the impliit prtF RF yn the other hndD if the kwrd di'usion is dominnt @ieF S f p < S b p AD we need to smooth the reonstruted solution for stilityD using the following formul

φp = 1 4 q∈N (p) (ϕ p + ϕ q + ϕ 1 pq + ϕ 2 pq ), φpq = 1 4 q∈N (p)
ϕ pq . @RFIVA sn se of step RD the di'usion oe0ients re reomputed using the new reonstruE tionsF elgorithm S summrizes the steps to onstrut the liner system t eh time stepF he resulting mtrix is pentdigonl mtrix s in the gut gell method oupled with grnkExiolson sheme @see etionQFRAF

Numerical validation

his sheme llows us the relxtion of the gpv onditionD usully present in most of level set methodsF o hek the roustness of the method for τ /∆ 2 rtio exeeding the usul gpv onditionD we onsider n initil level set funtion ϕ 0 in 2 × 2 dominD given y ϕ 0 (x, y) = x 2 + y 2 -R, @RFIWA with R = 0.8 the rdius of the HElevel setF he veloity (eld is set to F = -1 in the whole domin nd we run the simultion until the (nl time t f = 0.3625 for di'erent resolutions x a ITD QPD TR nd di'erent gpv numers rnging from I to ITF he L 2 norm of the error with respet to the nlytil solution is omputed in the whole dominF esults in pigure RFS show seondEorder ury for ny given gpv numer for the retrting irle seF 60 Chapter 4. Numerical methods for forward and adjoint problems Algorithm 5: porwrdEfkwrd di'usion sheme input :

ϕ n 1 D F output: ϕ n et the oundry onditions for ϕ n 1 foreach element x p ∈ Ω\∂Ω do se the reonstrutions φp = ϕ p , φpq = 1 2 (ϕ p + ϕ q )
gompute the rmiltonin using the dimondEell strtegy

|∇ϕ pq | = ϕ p -ϕ q ∆ 2 + ϕ 1 pq -ϕ 2 pq ∆ 1/2
gompute the forwrd nd kwrd di'usion oe0ients

a pq = F ( φp -φpq ) |∇ϕ pq | a f pq = max(a pq , 0), a b pq = min(a pq , 0) gompute the dominnt ontriution S f p = q∈N (p) a f pq , S b p = q∈N (p) a b pq if S f p < S b p then se the reonstrutions φp = 1 4 q∈N (p) (ϕ p + ϕ q + ϕ 1 pq + ϕ 2 pq ), φpq = 1 4 q∈N (p)
ϕ pq eompute the forwrd nd kwrd di'usion oe0ients

a pq = F ( φp -φpq ) |∇ϕ pq | a f pq = max(a pq , 0), a b pq = min(a pq , 0)
end end olve the liner system hepending on the seD usully fter one or more itertions of the time dvnement shemeD the level set funtion will ese to e n ext signed distne funtionD neessitting reinitiliztion step to enfore this riterion t the front Γ (t)F qiven funtion ϕ 0 D whih is not signed distne funtionD it n e evolved into funtion ϕ y solving the iikonl eqution

ϕ n p + τ ∆ 2 q∈N (p) a f pq (ϕ n p -ϕ n q ) = ϕ n 1 p + τ ∆ 2 q∈N (p) a b pq (ϕ n 1 p -ϕ n 1 q )
   ∂ϕ ∂t ⋆ = S (ϕ 0 ) (1 -|∇ϕ|) in Ω, ϕ(x, τ ) = 0 on Γ (t) , @RFPHA
where t ⋆ gin denotes pseudoEtimeD S(ϕ 0 ) is the signture funtion @iqution RFSAD nd ϕ the signed distne funtion @one stedyEstte is rehedAF here exist mny numeril methods for solving the iikonl eqution ndD here we hve dopted tht of win win PHIHD lso reently used in the work of vimre vimre et al. PHPQF his pproh relies on seond order ixy sptil disretiztion with suEell resolution ner the interfeF he oneEsided ixy (nite di'erene @in the x diretion onlyA yields

D + x ϕ ij = ϕ i+1,j -ϕ ij ∆ - ∆ 2 minmod (D xx ϕ ij , D xx ϕ i+1,j ) , D - x ϕ ij = ϕ i,j -ϕ i-1,j ∆ + ∆ 2 minmod (D xx ϕ ij , D xx ϕ i-1,j ) , @RFPIA
where

D xx ϕ ij = (ϕ i-1,j -2ϕ ij + ϕ i+1,j
) /∆ 2 is the seond order derivtive of ϕ ij nd the minmod9 limiter is zero when oth rguments hve opposite signsD nd tkes the rgument of smllest solute vlue otherwiseF he numeril rmiltonin |ϕ| is omputed s follows

|∇ϕ| ij ≃ H G D + x ϕ ij , D - x ϕ ij , D + y ϕ ij , D - y ϕ ij , @RFPPA
where funtion H g is given y

H G (a, b, c, d) =        max (a -) 2 , (b + ) 2 + max (c -) 2 , (d + ) 2 when S ϕ 0 ⩾ 0, max (a + ) 2 , (b -) 2 + max (c + ) 2 , (d -) 2
when S ϕ 0 < 0.

@RFPQA xer the interfeD the (nite di'erenes need to e modi(ed in order to impose ϕ = 0 where ϕ 0 = 0F o this endD qudrti ixy polynomil interpoltion gives

D + x ϕ ij = 0 -ϕ ij ∆ + - ∆ + 2 minmod (D xx ϕ ij , D xx ϕ i+1,j ) , @RFPRA
where

∆ + =              ∆   1 2 + ϕ 0 ij -ϕ 0 i+1,j -S ϕ 0 ij -ϕ 0 i+1,j √ D ϕ 0 xx   where ϕ 0 xx > ε, ∆ ϕ 0 ij ϕ 0 ij -ϕ 0 i+1,j elsewhere, @RFPSA with ϕ 0 xx = minmod ϕ 0 i-1j -2ϕ 0 ij + ϕ 0 i+1,j , ϕ 0 ij -2ϕ 0 i-1,j + ϕ 0 i-2j , nd D = ϕ 0 xx /2 -ϕ 0 ij -ϕ 0 i-1,j 2 -4ϕ 0 ij ϕ 0 i-1,j .
he negtive oneEsided ixy di'erene D - x is otined similrlyF e forwrd iuler sheme is then used for time dvnementF his rmiltonin extrtion will lso e used for hekEpointing in the djoint prolem @iqutions eFe nd eFfA when neededF e vlidte the method y initilizing pertured solutionD similr to the test se found in win PHIHD in 4 × 4 domin Stencil at time t n X full cells X mixed cells X empty cells n RF he vlues T * A nd T * B re interpolted using

T 1 B T 1 A T 2 B T 2 A T 3 B T 3 A n T new T B T A Stencil at time t n+1
T 1 A , T 2 A , T 2 A nd T 1 B , T 2 B , T 2 B respeE tivelyF
SF he oordintes o = (x new , y new ) of the ryenter of the ell to initilize is lotedF he lst step is done y using the disrete fe pities de(ned in etion QFQF en orthonorml oordinte system R = (o, (x ′ , n)) similr to the prolE(tted urvture found in opinet PHHW is then de(nedF pinllyD the new temperture vlueD T new D is linerly extrpolted y solving for a 0 nd a 1 in the resulting system of equtions

T * A = a 0 x ′ A + a 1 , T * B = a 0 x ′ B + a 1 , @RFPUA
where a 1 = T new F 4.5 Algorithm for solving the adjoint Stefan problem e present here the lgorithm for solving the djoint prolem @iqution eAF he djoint vriles Θ nd ψ re solved kwrd in time in order to ompute the grdient eqution @iqution eFhAF hese vriles re initilized with the forwrd vriles t (nl timeF es mentioned in ghpter PD the djoint level set ψ is no longer signed distne funtion ut n uxiliry vrile tht enters in the hirihlet 4.5. Algorithm for solving the adjoint Stefan problem 65 oundry ondition of the djoint temperture Θ t the interfeF he lotion of the interfe t given time step n is known thnks to ϕ n D n informtion tht is stored during the forwrd resolution of the prolem @elgorithm RAF woreoverD the onstruted vplin t the sme time step nD omputed thnks to elogrithm QD is reEused to solve the djoint het equtionsF he steps required to solve the djoint prolem re desried in elgorithm TF Algorithm 6: elgorithm for the solution of the djoint prolem input : T, ϕ, ε V , t f output: Θ 0 snitilize the djoint temperture (eld Θ f @iqution eFA snitilize the djoint level set funtion ψ f @iqution eFgA et the oundry onditions for Θ @iqution eFdA repeat 1. ghekpoint the fesD volumes nd entroids with ϕ n 2. gompute the djoint tefn ondition

[∇ Θn ] 1 2 ← Θ n
3. olve the djoint dvetion eqution

ψ n-1 ← ψ n , T n , ϕ n , v n s , [∇ Θn ] 1 2 , ε V @iqution eFfA
4. olve the djoint het equtions Θn-1 ← Θn , ψ n-1 , ϕ n 5. glen or snitilize ded or fresh ells Θ n-1 ← Θn-1 until initial time 0 ← t f Y pinllyD we fous on the solution of the (rstEorder onservtion lw @iqution eFfA on ΓF olving this (rstEorder eqution on time time moving surfe n e numeril hllengeF portuntelyD ll the derivtives n e extended on ll of ΩF veding to the following equtionD

∂ψ ∂t + div(ψ⃗ v S ) = rhs in Ω, @RFPVA
where the right hnd side

rhs = 1 |∇ϕ| ∂T ∂n -ε V ∂v S ∂n [∇Θ i ] 1 2 • n,
is omputed eforehndF xumerillyD iqution RFPV is solved using the sme porwrdEfkwrd di'usion sheme presented in etion RFQFIF ith thisD we lose the ghpter on numeril methods for oth the forwrd nd djoint tefn prolemsF sn the next ghpterD we present numeril results of the forwrd prolem s well s vlidtions on nlytil solutions of the twoEphse tefn prolemF BroadcastingX pplies given funtion elementEwise on multiEdimensionl rry without using extr memoryF Shared-memory parallelizationX esy implementtion multiEthreding in the min loopsF 68 Chapter 5. Numerical simulations of two-phase Stefan problems he numeril shemes presented in ghpter R re oded in n inEhouse pkge Flower.jl @see eppendix g for more detilsAF ih time step of the forwrd or djoint twoEphse tefn prolem requires the solution of 2 + 1 liner systemsX IF he two het equtions @one for eh phseA using the gut gell method oupled to the semiEimpliit grnkExiolson sheme @iqution QFISAF PF he level set dvetion eqution using the porwrdEfkwrd di'usion shemeD either in the forwrd prolem @iqution RFIRA or in the djoint se @iquE tion RFPVAF he liner systems re solved using the gonjugte qrdient itertive method imE plemented in the Iterative.jl pkgeF sn the se of )uid )owD the visous term in the xvier tokes equtions @see eppendix fA is solved using i @ortleD ixtensile oolkit for ienti( gomputtionA fly et al. IWWUF his g lirry is wrpped in the tuli lnguge in PETSc.jl F he optimiztion ses presented in ghpter T re rried out using the vEfpq method @djointEsedA or the rtiule wrm lgorithm @derivtiveEfreeA impleE mented in the Optim.jl pkge wogensen 8 iseth PHIVF 5.2 Analytical solutions of two-dimensional problems 5.2.1 Planar interface e vlidte our numeril methods for tefn prolems ginst nlytil solutionsF es (rst enhmrkD we onsider moving plnr interfe ounded y two wllsF his se is tken from vimre et al. PHPQD where they enhmrk their hyrid emeddedElevel set oundry methodD nd from grnk IWVU @originl test seAF en initil ie lyer is melting from the top where we impose positive temperture T 1 wall F he melting temperture is set to T M = 0 nd we only solve for the liquid phseF he tefn numer is set to

St = c L (T 1 wall -T M ) L H = 2.85,
where c L is the liquid het pity nd L H the ltent hetF he initil temperture in the liquid phseD only depending in the y oordinteD is set to

T L (y, t 0 ) = erf 1 -y 2 √ t 0 erf(l) , @SFIA
where t 0 is the initil time t whih the simultion is strted nd l = 0.9F he nlytil position of the interfe is given y the following formul he simultion is strted t t 0 = 0.03 suh tht there re enough ells to ompute the tefn ondition @iqution RFIAF he simultions is run until (nl time t f = 0.1 for di'erent numer of points per dimension x a ITD QP nd IPVF he rtio τ /∆ 2 = 0.5 is kept onstnt s we inrese the sptil resolutionF sn pigure SFID we show the L 2 nd L ∞ error norm of the temperture (eld with respet to the nlytil solution @iqution SFIAF he method shows n seond order onvergene in oth normsF his se vlidtes the initiliztion of fresh ells sineD s the ie lyer meltsD it unveils previously nonEexisting ells in the liquid phse whih need to e properly initilized @freshEell prolemAF 5.2.2 Growing Frank's spheres he growth of n ie sphere surrounded y n underEooled liquid ws initilly studied y prnk prnk IWSH where he showed tht the initil sphere rdius evolved in selfEsimilr mnner s the squreEroot of timeF his prolem n e used to vlidte the ury of the numeril sheme y ompring the results with n nlytil solutionF sn this on(gurtionD the temperture (eld is given y

y = 1 -2l √ t. @SFPA
T (r, t) = T (s) =    T ∞ 1 - F (s) F (S) , s > S, 0, s < S, @SFQA
where r = x 2 + y 2 nd s = r/t 1/2 F T ∞ is given underEooling tempertureD nd F (s) denotes the similrity solution of the het eqution

F (s) = E 1 (1/4 s 2 ).
es the sphere expndsD numeril errors my led to n unwnted ltertion of the initil shpe due to the wullinsEekerk instility wullins 8 ekerk IWTRD where pertured solution n led to unstle dendriti growth in the se of zero melting temperture @T M = 0AF e therefore test the roustness of our method using the initil prmeters reommended in elmgren elmgren IWWQF he initil level set funtion is irle of rdius R 0 = 1.56 in 8 × 8 domin surrounded y n initil negtive temperture (eld @iqution SFQ with T ∞ = -0.5AF he initil time is set to I nd the simultion is dvned until (nl time t f = 2F e rry out the simultions with di'erent grids x a QPD TR nd IPVF he results of this study re summrized in the following (guresF pigure SFP shows the onvergene of the rdius towrds the nlytil solutionF pigure SFQ shows the (nl error mp in temperture (eld s well s the interfe positions t di'erent times t = 1, 1.5, 2F e n lerly see tht the initil sphere tends towrds the nlytil (nl shpe s we inrese the numer of pointsF he intrinsi regulriztion of our method llows the level set funtion to retin its initil irulr shpeD voiding spurious interfe osilltionsF sn pigure SFRD we present the error norm in temperture (eld for the di'erent type of ellsF he order of onvergene is lose to IFS for ells tht re not mixedF he following test ses will fous on the dendriti growth of n initil rystl shpe in uniform underEooled temperture (eldF 5.3 Crystal growth in an under-cooled liquid bath 5.3.1 Surface tension eects grystl growth is n unstle phenomenon tht ours spontneously in ntureF sts pperne is the result of ompetition etween the nturl growth due to the tefn ondition @iqution IFPAD whih tends to streth the tipsD nd the surfe tension e'et present in the qisEhomson reltion @iqution IFQAD whih tends to restore the )tness of the interfe @see vnger IWVH for further detilsAF xumeril reprodution of suh ptterns is hllengeD sine the dendriti growth depends on the smllest resolved length sleF woreoverD in the sene of surfe tensionD the numeril noise nd the wullinsEekerk instility will led to n unphysil growth of the initil rystlF en exmple of rystl growth without surfe tension e'ets is shown in pigure SFSF he initil REfolded rystl expnds in n unordered fshion with numerilly produed dendritesF yn the other hndD when tking into ount the surfe tension e'etsD the solution is regulrized @pigure SFTAF sn this test seD we sses the e'et of the moleulr kineti oe0ientF e tke the sme initil ondition thn in etion SFQFID we (x the surfe tension ε κ = 0.001 nd we run the simultion until (nl time t f = 0.4F pigure SFW shows the e'et of the Chapter 5. Numerical simulations of two-phase Stefan problems kinemti oe0ient s we vry ε V from 0 to 0.0016F es the oe0ient is inresedD the tip of the dendrites move fster from the initil shpeF et the tipD where the initil veloity is mximlD due to the tefn onditionD the melting temperture is redued y n ε V v S ftor thus further entrining the growing motionF 5.3.4 Anisotropy eects sn order to urtely reprodue the rystl shpes produed y dendriti growthD we introdue vrile surfe tension oe0ientD similr to n 8 rs PHHT nd vimre et al. PHPQD in the qisEhomson reltion @iqution IFQA to ount for nisotropy e'ets

εκ (α) = ε κ 1 + A 8 3 sin 4 1 2 M (α -α 0 ) -1 , @SFSA
where ε κ is the surfe tension oe0ientD A represents the weight of the nisotropy e'etD M the mode numerD α the ngle of the interfe with respet to the xExis nd α 0 the presried ngle of symmetry on whih the dendrites will growF o vlidte the nisotropy e'etsD we initilize sixEfold rystl in n underE ooled liquid with T ∞ = -0.8 nd (xed veloity oe0ient ε V = 0F he nisotropi weight is (xed t A = 0.4 nd mode numer M = 6F he simultions re rried out for two di'erent presried ngles α 0 = π/2, π/4 until (nl time t f = 0.09 in 2 × 2 domin with N = 300F pigure SFIH shows the rystl growing in the diretion of the presried ngle α 0 = π/2D when the ngle of symmetry is ligned with the initil rnhesF he miniml surfe tension e'et long the symmetry lines llows the the primry rnhes to grow further wy from the initil rystlF e n lso oserve seE ondry rnhes ppering nd propgtingF hen α 0 = π/4 @di'erent from the initil orienttion of the rnhesA the rystl rottes towrds the presried diE retion @pigure SFIIAF hese results vlidte our implementtion of the nisotropy e'etsF

Ra e = Ra(1 -T M ) h3 , @SFTA
where h is the verge )uid height de(ned s

h(t) = 1 b ˆb 0 h(x, t)dx. @SFUA
hen the e'etive yleigh numer rehes the wellEknown ritil vlue Ra c = 1707.76D the initil di'usionEdriven motion is trnsformed into onvetionEdriven oneF e rry out simultions for di'erent Ra = 10 3 , 10 4 , 10 5 nd 10 6 with tefn numer set to St = 0.25F he spet rtio is set to b = 8 nd the grid used in our simultion is 512 × 64F he simultion re run until the mximum height rehes 0.9HF pigure SFIP shows the e'etive yleigh numer s funtion of time for the di'erent sesF hen Ra = 10 3 nd Ra = 10 4 the ritil yleigh numer is not rehed nd the motion of the )uid lyer is not 'etedD thus remining di'usionEdriven oneD similrly to the plnr motion desried in etion SFPFIF sn the Ra = 10 5 nd Ra = 10 6 sesD howeverD we see Ra e rossing the threshold inditing the onset of the instilityF his onset is visile in pigure SFIQD where the verge height strongly inreses when the onvetionEdriven motion tkes overF pigure SFIR shows time series of the temperture (eld nd interfe position for the Ra = 10 6 seF sn tht seD the ritil yleigh numer is rehed t t ⋆ = 0.1034 nd the (rst ifurtion ppersD reting the onvetion ellsF he size of the onvetion ells will then vry with the seondry ifurtions mehnismF hen the verged height h mthes the hrteristi wvelength of the onvetion rollsD the onvetion ells hve su0ient time to merge nd then stilizeF e lso note tht the interfe is deformed ording to the shpe of the ellsF

Actuation function parametrization

e rnge of optimiztion ses for the tefn prolem with vrying omplexity re presented hereF he results nd performne of the djointEsed lgorithm @elgoE rithm PA re nlyzed nd ompred to derivteEfree lgorithm @with no grdient requirementA for these di'erent setupsF sn ll of the ses the desired temperture (eldD T f D nd desired level setD ϕ f D re omputed eforehndD nd re used to drive the ontrol vrile w tht ts s xeumnn @or hirihletA oundry ondition on the domin oundryF he tutor w is prmetrized with pourier sis w = n p=1 a p cos (pπx) + n p=1 b p sin (pπx) , @TFIA where x orresponds to the ounds of the domin nd a nd b to the sis oe0E ientsF he numer of oe0ients to (t is 2nD with n eing the order of eh sis funtionF hrough the optimiztion proessD the mplitude of eh sis funtion is determined using the grdient eqution @iqution eFhAF fy opting for prmeterized distriution we ensure the smoothness of the tution funtionF hue to the high sensitivity of the ost funtionl with respet to the sis onsidered ! too mny prmeters will rete multiple lol minim ! the numer of prmeters re kept t low enough vlue to ensure the onvexity of the prolem while llowing sptil vrition of the tution funtionF en exmple of ost funtionl J s funtion of the two oe0ients tht prmetrize the tutor is shown in pigure TFIF sn tht simple seD the onvexity of the ost funtionl is ler nd the grdient eqution will drive J towrds the glol minimumF woreoverD in ll of the optimiztion ses onsideredD the initil guess ! the strting point in our ost funtionl spe ! is set to w = 0 in order not to dd is regrding the initil desentF 6.2 Adjoint-based optimization cases 6.2.1 Asymmetric melting of a solid circle sn the (rst seD we onsider n symmetri expnsion nd retrtion of n initil ie irle of rdius R = 0.75 in 2×2 dominF he surfe tension oe0ient is onstnt nd set to ε κ = 0.002 while ε V = 0F he nonEzero surfe tension oe0ient is dded to further regulrize the level set funtionF he numer of points per dimension is N = 64 nd the (nl time is set to t f = 0.1F he ontrol w ts on the whole domin oundry nd is of hirihlet typeF he djoint oundry onditionD on the other hndD remins homogeneous xeumnn oundry onditionF he tutor w is prmeterized s u = 2 p=1 a p cos (pπx) + 2 p=1 b p sin (pπx) . @TFPA 6.2. Adjoint-based optimization cases 87 he oe0ients in the ost funtionl @iqution PFIA re set to β 1 = 1D β 2 = 10 nd β 4 = 10 3 F sn this seD the term ontrolling the interfe length β 3 is set to zeroF pigure TFP shows the level set of the forwrd prolem evolving towrds the desired shpe fter IT itertions orresponding to TT solutions of the forwrd prolem @pA nd IT solutions of the djoint prolem @eAF he lgorithm is le to reover oth the retrtion nd the expnsion of the initil irleF Figure 6.2: Iterations 0, 3, 10 and 16 of the optimization procedure for the circle test case. The blue curve represents the desired shape and the red one the nal position of the interface at a given iteration. The nal time is t f = 0.1 and the interface is plotted with a time step of 0.01. The color map corresponds to the error in temperature eld. The inset shows the actuator w at a given iteration (red) and the desired one (blue).

Chapter 6. Shape optimization results

pigure TFQ shows the evolution of the normlized ost funtionl J /J 0 s funtion of the itertionsF et itertion TD we see the desent slowing down nd lmost stoppingD inditing the presene of lol minimumF xeverthelessD the lgorithm mnges to overome it nd ontinue its desentF

Controlling the Mullins-Sekerka instability

sn this seond seD we onsider pertured initil plnr shpe in 2 × 2 domin suh tht the initil level set funtion is given yD ϕ 0 (x, y) = y + 0.6 + A cos(2πx), @TFQA nd the temperture (eld yD

T 0 (x, y) = -1 + e -T∞ϕ 0 (x,y) , ϕ 0 > 0, 0, ϕ 0 < 0, @TFRA
with n mplitude of the perturtion A = 0.05 nd n underEooling temperture T ∞ = 1.2F he surfe tension nd veloity oe0ients re set to zero in this seF es desried in ghpter SD n initil perturtion suh s tht presried in this prolem leds to wullinsEekerk type instility wullins 8 ekerk IWTRD ghen et al. IWWU hrterized y unstle dendriti growthF he purpose of this optimiztion test se is to redue the instility y imposing n optiml tution on the top oundry ounterEting the motion of the tipsF he ontrol vrile w is of xeumnn type nd we set the (nl time t f = 0.5 with N = 64F imilr to the previous se the tution is prmetrized s sn pigure TFRD we n oserve how the initil tip splitting is redued s we go through the optimiztion proedureF he minimum is ttined fter IH itertions of the optimiztion lgorithm orresponding to SR solutions of the forwrd prolem @pA nd IH solutions of the djoint prolem @eAF Figure 6.5: Normalized cost functional as a function of the iterations of the L-BFGS optimization procedure for the Mullins-Sekerka case.

he vrition of the ost funtionl @pigure TFSA is similr to the previous oneD with n initil steep desent fter few itertions tht eventully plteusF sn this seD the optimiztion lgorithm is stopped y the riterion on the grdient |∇J | < 10 6 F 6.2.3 Growing crystals with topology changes sn this lst seD we onsider three rystls symmetrilly disposed in n underE ooled liquidF he rystls will grow nd eventully mergeF he ojetive of this optimiztion is to drive the (nl shpe towrds the desired one y ting on the oundries of the whole dominD nd therey to suppress the nisotropy e'etsF he surfe tension oe0ient is set to ε κ = 0.0005 nd the veloity oe0ient to ε V = 0.002F e hoose reltively smll surfe tension oe0ient with respet to the veE loity one in order to llow for strong dendriti formtion nd to exmine the optiE miztion lgorithm in hllenging se where the topology remins omplexF he nisotropy e'ets re dded y setting α 0 = π/2 nd M = 4 in iqution SFSF he het )ux indued through tution will hve to ompete with these e'ets in order to drive the interfe towrds the desired oneF 6.2. Adjoint-based optimization cases 91 he simultions re run in 4 × 4 domin with N = 100 until (nl time t f = 0.45F he underEooling initil temperture is set to T ∞ = -0.6 nd the ontrol vrile w is of xeumnn typeF sn ontrry to the previous seD we need to restrit the sis used for the tutor @y (tting only two prmetersA in order to desend towrds the glol minimumF he following prmetriztion is used u(x, p) = p 1 ((1 + cos(π/8x))/2) 4 + p 2 ((1 + sin(π/8x))/2) 4 . @TFTA Figure 6.6: Iterations 0, 15, 17 and 20 of the optimization procedure for the growing crystals test case. The blue curve represents the desired shape and the red one the nal position of the interface at a given iteration. The nal time is t f = 0.45 and the interface is plotted with a time step of 0.045. The color map corresponds to the error in the temperature eld. The inset shows the actuator w at a given iteration (red) and the desired one (blue).

92

Chapter 6. Shape optimization results

pigure TFT shows the evolution of the surfeF he minimum is ttined fter PH itertions of the optimiztion lgorithm orresponding to PS solutions of the forwrd prolem @pA nd PH solutions of the djoint prolem @eAF st n e seen tht initilly the rystl is driven towrds the domin orners y the nisotropy prmetersF hrough the optimiztion proedure howeverD the (nl rystl shpe tends towrds the desired oneF woreoverD the topology hnges re impliitly tken into ount thnks to the level set methodF vooking t the vrition of the ost funtionl through the optimiztion proeE dure @pigure TFUAD we n notie n initil plteu until itertion ISF his is used y the lrge di'erene in the initil shpe ! t itertion H with w = 0 ! nd the desired oneF fy strting elsewhere in the prmeter spe of J D the desent towrds the minimum might our fsterF le TFI summrizes the results otined for eh seF gse yptimiztion prmeters esults sn the previous setionD we showed tht the djointEsed optimiztion proedure ppers to e roust lgorithm to ontrol the shpe of melting or solidi(tion frontD even in the presene of dendriti instilities nd nisotropi e'etsF e now ompre the results to derivtiveEfree method @where no informtion on the grdient is requiredAD the rtile wrm optimiztion lgorithm hn et al. PHHWF his method ttempts to improve glol onvergene y swithing etween four evolutionry sttesX explortionD exploittionD onvergeneD nd jumping outF sn the jumping out stte it tries to tke the est prtile nd move it wy from its lol optimumD to improve the ility to (nd glol oneF e test the derivtiveEfree optimiztion proedure on the sme three ses presented previously using the sme prmeters β 1 to β 4 for the ost funtionlF le TFP summrizes the omprison etween oth methodsF es expetedD the derivtiveEfree one ttins lower minimum t ost of muh higher numer of funtion evlution @one order of mgnitude higherAF gse vEfpq rtile wrm J lls ∇J lls 7.1 Two-phase ows in the one-uid formulation sn the (rst prt of this workD we studied melting nd solidi(tion proessesD modeled y the twoEphse tefn prolem @ghpter IAD where the motion of the interfe ws solely funtion of the het trnsfer ross the interfe etween )uid nd solid phseF sn tht seD when onsidering the )ow in the )uid phseD we solved the xvierEtokes equtions in the foussinesq pproximtionF sn this seond prt of the disserttionD we will study physil prolems governed y the inompressile twoEphse xvierEtokes equtionsF plows of immisile )uids re uiquitous in xture nd in everydy lifeY for exE mple in se wvesD rin drops or in ules rising in glss of sprkling wterF he irEwter )ows re hrterized y lrge density nd visosity rtios s well s high surfe tensionF his interfil fore is n energy per re tht resists the retion of new interfeF hese disprities in the mteril properties ross the interfe will tend to generte omplex phenomen involving di'erent sptil nd tempoE rl slesF he prtiulr se of the motion of )uidE)uid interfe over solid sustrte hikhmurzev IWWRD flke et al. IWWWD

β 1 β 2 β 3 β 4 J lls ∇J lls J f inal /J 0 welting
J f inal /J 0 J lls J f inal /
ρ (∂ t u + u • ∇u) = -∇p + ∇ • (2µD) + σκδ s n + g, @UFIA
supplemented y the divergeneEfree ondition ∇ • u = 0, @UFPA nd trnsport eqution for the density

∂ t ρ + ∇ • (ρu) = 0, @UFQA
with u the )uid veloityD ρ the )uid densityD µ the )uid visosity nd D the deforE mtion tensor de(ned s D ij = (∂ i u j + ∂ j u i )/2D σ the surfe tension oe0ientD κ the urvtureD δ s the hir distriution expressing the ft tht the surfe tension term is onentrted t the interfeD n the unit norml to the interfe nd g the elertion of grvityF sn the oneE)uid formultionD olor funtion c(x, t) is used to hrterize the interfe lotion nd the )uid properties @density nd visosityAF he volume frE tion c(x, t) will e hosen s the volume frtion oupied y referene phseF he density nd visosity re now de(ned y

Fluid 1 Fluid 2 ρ = ρ 1 µ = µ 1 ρ = ρ 2 µ = µ 2 c = 1 c = 0 c = 0.5 ρ = 0.5ρ 1 + 0.5ρ 2 µ = 0.5µ 1 + 0.5µ 2
ρ(c) ≡ cρ 1 + (1 -c)ρ 2 , µ(c) ≡ cµ 1 + (1 -c)µ 2 , @UFRA
with ρ 1 D ρ 2 nd µ 1 D µ 2 the densities nd visosities of the (rst nd seond )uids respetivelyF pigure UFI shows n exmple of the disrete representtion of the density nd visosityF he dvetion of the density @iqution UFQA is then nturlly repled y the dvetion eqution for the olor funtion

∂ t c + ∇ • (cu) = 0. @UFSA
sn prt I of this workD this eqution ws solved using n impliitEexpliit sheme to dvet level set funtion @ghpter RAF sn the ontext of the twoEphse xvierE tokes equtionsD when onsidering the shrpEinterfe model ! s opposed to the di'use interfe model desried in the following setion ! iqution UFS will e solved using onservtiveD nonEdi'usive geometri olumeEyfEpluid @ypA shemeF he yp method nd the tretment of the oundry onditions will e detiled in etions WFI nd WFPF 7.2 Phase-Field extensions sn the hseEpield methodD the governing equtions re derived from the thermoE dynmi potentils of the systemD together with the ssumption of surfe energy ssoited with n interfe emerg PHHQF st isD thereforeD possile to onsider difE ferent physil situtions with reltive eseF st is lso strightforwrd to implement numerillyD sine interfes re not trked expliitlyF snstedD vrile is introE dued tht hs di'erent onstnt vlues in the two phses with steep trnsition etween the two in the di'use interfe tqmin PHHHD grlson PHIPF yne of the mjor drwks of this methodD howeverD is tht width of the interfe ε must e smll to mth the proper interfe dynmisF imilrly to the oneE)uid formultionD the hseEpield model introdues phse vrile C(x, t) rnging from 1 to -1F his phseE(eld vrile C(x, t) is governed y onvetionEdi'usion eqution

∂C ∂t = F d -F c , @UFTA
where F d is the di'usive )ux nd F c is the onvetive )uxF he ltterD in the ontext of inompressile )owsD tkes the simple form

F c = u • ∇C, @UFUA nd the di'usive )ux is F d = -M ∇ϕ, @UFVA
where M is proportionlity oe0ient lled the hseEpield moility nd ϕ the hemil potentil de(ned s

ϕ = βΨ ′ (C) -α∇ 2 C. @UFWA
sn the hemil potentilD we hve two prmeters α nd βD whih re relted to the surfe tension σ nd the hrteristi thikness of the di'use interfe ε s

σ = 2 3 2αβ, ε = α β . @UFIHA 98
Chapter 7. Two-phase Navier-Stokes equations and Phase-Field extensions sn dditionD it ontins the derivtive of the stndrd douleEwell potentil

Ψ (C) = (C + 1) 2 (C -1) 2 4 . @UFIIA
ith the de(nitions ove @iqutions UFT to UFIIAD the onvetionEdi'usion equE tionD referred to s the ghnErillird eqution ghn 8 rillird IWSVD n e writE ten s

∂C ∂t = ∇ • M ∇ βΨ ′ (C) -α∇ 2 C -u • ∇C. @UFIPA
iqution UFIP is the hseEpield nlog of iqution UFSF he oundry onditions for this fourthEorder prtil di'erentil eqution will e detiled in our study of dynmi wetting y the hseEpield method @etion IHFQFIAF king the sme nottion s in iqution UFID the momentum eqution in the xvierEtokes equtions eomes ρ

(∂ t u + u • ∇u) = -∇p + ∇ • (2µD) -C ∇ϕ + g, @UFIQA
where -C ∇ϕ orresponds to the surfe tension fore nd ts over the di'use interfe regionF his form of the surfe tension foring is the so lled potentil form tqmin IWWWF woreoverD the density nd visosity re now de(ned through the hseEpield vrile C s

ρ(C) ≡ ρ 1 C + 1 2 + ρ 2 C -1 2 , µ(C) ≡ µ 1 C + 1 2 + µ 2 C -1 2 .
@UFIRA e now stte the oundry onditions in the presene of ontt line @see ghpter VA on solidD hrterized y its equilirium ngle @the ngle etween the )uidE)uid interfe nd the wllAF he onvetionEdi'usion equE tion @iqution UFIPA is fourthEorder prtil di'erentil eqution nd requires two oundry onditionsF pirstD we impose nonEequilirium wetting ondiE tion tqmin PHHHD in et al. PHHQ on the solid wllD

-µ f ε ∂C ∂t + u • ∇C = α∇C • n -σ cos θ e g ′ (C) , @UFISA
where µ f is ontt line frition prmeterD hving the sme units s ulk dynmi visosityF rereD θ e is the equilirium ontt ngle nd g (C) = 0.5 -0.75C + 0.25C 3 , @UFITA is swithing funtion desriing smooth trnsition etween oth phsesF he unit norml vetor n is direted from the )uid to the surrounding solidF sf one sets µ f = 0D the ontt ngle is lwys enfored to the equilirium ngle θ e F his will e the se in the toy model with ompring p to yp simultions @etion IHFQAF xonEzero ontt line frition llows the dynmi ontt ngle to evolve nturlly s funtion of ontt line speedF he seond oundry ondition for the phse funtion is zero di'usive )ux of hemil potentil through the oundries ∇ϕ • n = 0. @UFIUA sn the following hpterD we will present the min topi of our studyX the dynmi ontt linesD the singulrities involved nd the models to urtely simulte the motion of the )uidE)uidEsolid intersetion pointF Chapter 8

Singularities in ows with contact lines sn the stti seD for droplet t equilirium @pigure VFIAD the equilirium ontt ngle is well known to e determined y oung9s lw

σ SG -σ SL = σ LG cos(θ e ) @VFIA
where σ SG D σ SL nd σ LG re the solidEgsD solidEliquid nd liquidEgs surfe tensions respetivelyF xumerillyD while still eing sujet spurious urrents depending on the hrteristis of the system @see eppendix hAD the stti se does not pose ny prtiulr issuesF 8.2 Paradoxes and singularities in the dynamic case he di0ulties rise when onsidering dynmi ontt lineF sn ftD the motion of the ontt line posesD sine ruh nd riven ruh 8 riven IWUI remrkle prolem euse of the ontrdition etween the noEslip ondition on the sustrte nd the motion of the ontt lineF nder these ssumptionsD one nturlly ends up with n immoile ontt line tht is inonsistent with simple oservtions @for exmple wter drop sliding on lefAF his noEslip prdox indues stress singulrity t the ontt lineD summrized in the words of ruh nd rivenX xot even rerules ould sink solid FFF his prdox does not exist in di'use interfe modelsD suh s hseEpield @eE tion UFPAD s the ontt line dvnes through di'usion even in the noEslip senrioF rowever speil tretment of the oundry ondition t the wll is required for shrp interfe modelsD suh s olumeEyfEpluid @etion UFIAF yne the (rst methE ods tht ws proposed to relieve the fore singulrity t the ontt line ws to simply relx the noEslip onditionF sn the slip length theoryD expressed s xvier oundry ondition @xfgAD termD proportionl to the sher veloity pro(leD is dded in the oundry ondition suh tht

u x | y=0 -λ ∂u x ∂y | y=0 = 0, @VFPA
where u x | y=0 is the tngentil veloity t the wll @loted t y = 0A nd λ the slip lengthF rnslting this oundry ondition into reltion for the ontt line speed yields

U CL ∼ λ ∂u ∂y | y=0 , @VFQA
where U CL is the ontt line speedF his simple slip model hs proven to e su0ient to solve the ontt line motion in some hydrodyE nmi ses ilson et al.

U CL ∼ λ ∂u ∂y | y=0 + σ µ (cos θ d -cos θ e ) , @VFWA
where σ is the surfe tension nd µ the visosity of the primry )uidF his oundry ondition will e tested in withdrwing plte setup where n urte modeling of the dynmi ontt line is ruil to predit the onset wetting filureF sn the following hpters we will present the implementtion of the models nd their pplition to physil setupsF sn ghpter WD the lssil slip nd dynmi ngle models will e desried in depth nd tested on di'erent prolemsF sn ghpE ter IHD the superEslip nd qxfg models will e presented s well s the toy modelD where we will relte the mss )ux present in the hseEpield model to simple slip nd dynmi ngle reltion in the yp methodF

ρ (∂ t u + u • ∇u) = -∇p + ∇ • (2µD) + σκδ s n + g, ∂ t ρ + ∇ • (ρu) = 0, ∇ • u = 0 @WFIA
with u the )uid veloityD ρ the )uid densityD µ the )uid visosity nd D the deforE mtion tensor de(ned s D ij = (∂ i u j + ∂ j u i )/2D σ the surfe tension oe0ientD κ the urvtureD δ s the hir distriution funtion used for the shrp interfe modelD 106 Chapter 9. Numerical methods for contact lines n the unit norml to the interfe nd g the elertion of grvityF por twoEphse )owD the volume frtion c(x, t) is de(ned s the integrl of the (rst )uid9s hrteristi funtion in the ontrol volumeF he volume frtion c(x, t) is used to de(ne the density nd visosity in the ontrol volume

ρ(c) ≡ cρ l + (1 -c)ρ g , µ(c) ≡ cµ l + (1 -c)µ g , @WFPA
with ρ l D ρ g nd µ l D µ g the densities nd visosities of the liquid phse nd the gs phse respetivelyF he dvetion eqution for the density is then repled y the eqution for the volume frtion

∂ t c + ∇ • (cu) = 0
he projetion method is used to solve the inompressile xvierEtokes equtions omined with fellEgollelEqlz dvetion sheme nd yp method for interE fe trkingF he resolution of the surfe tension term is diretly dependent on the ury of the urvture lultionF he reightEpuntion methodology is ypEsed tehE nique for lulting interfe normls nd urvtures efkhmi 8 fussmnn PHHVD efkhmi 8 fussmnn PHHWF eout eh interfe ellD )uid heights9 re lulted y summing )uid volume in the grid diretion losest to the norml of the interfeF sn two dimensionsD 7 × 3 stenil round n interfe ell is onstruted nd the heights re evluted y summing volume frtions horizontlly @pigure WFIA

h j = k=i+3 k=i-3 c j,k ∆,
with c j,k the volume frtion nd ∆ the grid spingF he heights re then used to 9.1. The Volume-Of-Fluid method 107 ompute the the interfe norml n nd the urvture κ n = (h x , -1),

κ = h xx (1 + h 2 x ) 3/2 , @WFQA
where h x nd h xx re disretized using seondEorder entrl di'erenesF st is importnt to note tht numeril spei(tion of the ontt ngle 'ets the overll )ow lultion in two wysX st de(nes the orienttion of the yp reonstrution in ells tht ontin the ontt lineF st in)uenes the lultion of the surfe tension term y 'eting the urE vture omputed in ells t nd ner the ontt lineF he orienttion of the interfeD hrterized y the ontt ngle ! the ngle eE tween the norml to the interfe t the ontt line nd the norml to the solid oundry ! is imposed in the ontt line ellF e now present some of the hrteristis of the olumeEofEpluid xvierEtokes solverF por further detilsD we refer the reder to rdovelli 8 leski IWWWD opinet 8 leski IWWWD opinet PHHWD opinet PHISD opinet PHIVF e stggered in time disretiztion of the volumeEfrtionGdensity nd pressure omined with timeEsplitting projetion method leds to the following time disretiztion

ρ n+ 1 2 u * -u n ∆t + u n+ 1 2 • ∇u n+ 1 2 = ∇ • µ n+ 1 2 (D n + D * ) + (σκδ s n) n+ 1 2 , c n+ 1 2 -c n-1 2 ∆t + ∇ • (c n u n ) = 0, u n+1 = u ⋆ - ∆t ρ n+ 1 2 ∇p n+ 1 2 ,
@WFRA whih requires the solution of the oisson eqution

∇ • ∆t ρ n+ 1 2 ∇p n+ 1 2 = ∇ • u ⋆ . @WFSA he momentum eqution n e rewritten s ρ n+ 1 2 ∆t u ⋆ -∇ • µ n+ 1 2 D ⋆ = ∇ • µ n+ 1 2 D n + (σκδ s n) n+ 1 2 + ρ n+ 1 2 u n ∆t -, u n+ 1 2 • ∇u n+ 1 2 , @WFTA
where the rightEhnd side depends only on vlues t time n nd n + 1/2F his equE tion is solved using multilevel oisson solverF he veloity dvetion term u ∇u n+ 1 2 is estimted using the fellEgollelEqlz fell et al. IWVWD opinet PHHQ seondEorder upwind shemeF pe is disretized using qudtree prtitioning in Ph @pigure WFPAF ell the vriles re olloted t the enter of eh squre disretiztion volumeF gonsistently with (niteEvolume formultionD the vriles re interpreted s the volumeEverged vlues for the orresponding disretiztion volumeF e projetion method is used for the sptil disretiztion of the pressure orretion eqution nd the ssoited divergene in the oisson equtionF o solve the dvetion eqution the geometril yp sheme is used nd proeeds in two stepsX IF snterfe reonstrutionF PF qeometril of )ux estimtion nd interfe dvetionF he reonstrution is pieewise liner interfe lultion9 @vsgAD followed y vgrngin dvetionF sn the vsg tehniqueD given volume frtion c(x, t) nd n pproximte norml vetor nD liner interfe is onstruted within eh interfe ellD whih orresponds extly to c(x, t) nd nF sn pigure WFQD we illustrte the priniple of geometril )ux estimtionF he totl volume whih will e )uxed to No-Slip

u i+ 1 /2j ∆t u i+ 1 /2j C C a
u x = 0 λ u x = λ ∂u x ∂y Slip u i j u i j+1 u i+1 j+1 u i+1 j u i-1 j+1 u i-1 j u i-1 j-1 u i j-1 u i+1 j-1
u ghost values are imposed Chapter 9. Numerical methods for contact lines tenils loted lose enough to the oundries of the domin will extend eyond it @pigure WFRAF he stenil vlues outside the domin @ghost vluesA need to e initilizedF hese vlues tht depend on the vlues inside the domin re set in order to provide the disrete equivlent of the xfg s follows

u x [ghost] + u x [ ] 2 + λ u x [ghost] -u x [ ] ∆ = U, ⇐⇒ u x [ghost] = 2∆ 2λ + ∆ U + 2λ -∆ 2λ + ∆ u x [ ],
with u x [ghost] the tngentil veloity t the ghost ellD u x [ ] the tngentil veloity of the ell inside the dominD ∆ the grid sping nd λ the slip lengthF es u is solution of the oisson equtionD we lso need to de(ne the disrete homoE geneous ounterprt of the xfg

u h x [ghost] = 2λ -∆ 2λ + ∆ u x [ ].
e will show tht the xvier oundry ondition oupled with onstnt ontt ngle is su0ient to model the ontt line motionF sn etion WFQFID we study the onset of wetting filure in urtin oting systemD where liquid flls on moving plte nd strts oting the surfeF e re le to reprodue the nonEmonotoni ehviour of the ritil veloity s the liquid )owErte inresesF woreoverD in v© is et al. PHPHD we enhmrked our yp simultions with phseE (eld @pA nd moleulr dynmis @whA simultions in nnosopi shered droplet setupF sn tht se we used the slip length nd the onstnt ontt ngle s (tting prmetersF he slip hd to e lolized ner the ontt line pointF he oundry ondition for the tngentil omponent in the xfg ws modi(ed s follows

u x | y=0 -f x ε λ ∂u x ∂y | y=0 = U, @WFVA
where f x ε is ell funtion de(ned s

f d ε =    1 + cos (πd/ε) 2 2 |d| < ε, 0 |d| ≥ ε, @WFWA
where d is the distne to the ontt line point nd ε the width of the ell funtion @relted to interfe width in the phseE(eld modelAF

Ca CL = µ l U CL σ , @WFIIA
where µ l is the )uid visosity nd σ the surfe tensionF hen the surrounding )uid is of lower visosityD the funtion g(θ) simpli(es to

g(θ) = ˆθ 0 x -sin x cos x 2 sin x dx. @WFIPA
sn etion WFQFPD we ompre our model with experiments of spreding squlne dropF 9.3 Application to physical systems 9.3.1 Onset of wetting failure in a curtain coating setup he (rst system onsidered is the urtin oting pulln et al. PHPHF sn the urE tin oting system @pigure WFSAD liquid is flling with veloity V on plte moving t veloity U F hen the liquid rehes the solid sustrteD it strts oting the free surfeD s shown in the time series exmple @pigure WFTAF e stedyEstte solution is only otined for given sets of physil prmeters nd the onset of wetting filure n e predited y studying rnge of pillry nd eynolds numers y vrying U nd V F e will mke the sme ssumption s viu viu et al. PHITX onstnt ontt ngle oupled with xvier oundry ondition is su0ient to model the dynmi wetting systemD provided tht the ir stresses re tken into ountF sn our modelD these stresses re diretly tken into ount in the twoEphse xvierEtokes solverF e onstnt ontt ngle θ m is imposed through height funtions efkhmi 8 fussmnn PHHVD efkhmi 8 fussmnn PHHW nd we use the xvier oundry ondition for the tngentil veloity t the wllF he multiEsle nture of the urtin oting on(gurtion n led to numeril di0ultiesD in prE tiulr on the resolution of the smllest length sleD the slip length λF h c d c ρ l , µ l ρ g , µ g θ m s Figure 9.5: Schematic of the curtain coating conguration. The system parameters are: h c the curtain height, d c the curtain width, ρ l , ρ g and µ l , µ g the densities and viscosities of the liquid phase and the gas phase respectively, U the substrate velocity, V the feed ow velocity and θ m the imposed contact angle. The inection point noted IP corresponds to the point at which the curvature of the interface is zero. The distance from the triple point to the IP will be used as a control quantity. woreoverD the regulriztion of the ontt line prdox is diretly dependent on the ury of the solution ner the ontt lineF he dptive mesh re(nement llows good resolution of the interfe t the triple pointF e will study the onvergene of the solution of the viu viu et al. PHIT on(gurtion s we inrese the numer of grid points per slip length λ/∆F he system will then e extended to ompre with experiments of flke et al. IWWW nd wrston et al. PHHWF es (rst stepD to vlidte our modelD we reprodue the urtin oting on(gurtion desried in viu et al. PHITF e onsider smll urtin height h c a 10 -2 m nd smll urtin width d c a 10 -3 m with lrge slip length λ a 10 -5 mF es the ury of the interfe reonstrution depends on the resolution of the smllest length sleD these previous onsidertions drstilly derese the omputtionl ost of this multiEsle prolemF he )uid propertiesD orresponding to glyerinted wter nd irD re the following

ρ l = 1000kgFm 3 , ρ g = 1.2kgFm 3 , µ l = 25mFs, µ g = 0.018mFs,
with the surfe tension set to σ = 70mxFm 1 F he visosity rtio µ g /µ l = 7.2 10 -4 is kept onstnt in our omputtionsF he sustrte veloity U is vried from HFI to IH mFs -1 nd the feed )ow veloity V from HFI to I mFs -1 F he dimensionless numers governing the )ow reX the pillry numer

Ca = µ l U σ ,
vrying from HFI to PFSD the eynolds numer

Re = ρ l V d c µ l ,
vrying from I to RH nd the fond numer

Bo = ρ l g σ d c V U 2 ,
vrying from 10 -3 to 10 -1 F he ontt ngle is kept onstntD θ m = 90°F o determine whether the simultion with given set of physil prmeters rehes the stedyEstte solutionD we set very lrge (nl time nd omputeD t eh time stepD the di'erene on the veloity (eld etween two susequent time stepsF sf the di'erene is lower thn given thresholdD we n onlude tht the )ow hs rehed stedy stte nd tht there is no wetting filureF sn pigure WFUD we show n exmple of yp simultion for grid ∆ = 0.156 µmD orresponding to TR grid spings per slip lengthF he dptive mesh re(nement llows good resolution of the interfe t the triple pointF sn this prtiulr seD for e a QH nd g a PFTD stedy stte solution is rehedF 114 Chapter 9. Numerical methods for contact lines 9.3. Application to physical systems 115 fy vrying the sustrte veloity nd the feed )ow veloityD we re le to reover the sme qulittive )ow on(gurtionsX ed pullingD right eneth the liquid urtin nd heel formtion @pigure WFVA nd similr oting window @pigure WFWA s in viu et al. PHITF he hydrodynmi ssist hs the most impt @ieF the moving plte veloity U is mximumA when the ontt line is eneth the liquid urtinF his on(gurtion llows stronger pressure due to the liquid inerti t the triple point preventing the formtion of ules nd therefore preventing the wetting filure from ourringF o determine the oting windowD we look for the (rst unstedy solution while inresing the g numer for given e numerF he error r reltes to the di'erene in g vlues etween the lst stedy solution nd the (rst unstedy oneF e onvergene study of the resolution of the interfe is onduted for this onE (gurtionF es the mximl level of re(nement is inresedD the resolution of the interfe t the ontt line is improvedF he mirosopi ontt ngle tends to the presried one of 90°s the smllest ell size is deresed from IH µm to HFIST µmD orresponding to numer of grid points per slip length inresing from I to TRF e hoose the TR grid points per slip length solution s the referene solution for the ontt line positionF sn pigure WFIHD we plotted the reltive error of the ontt line position nd the reltive error of the distne from the ontt line position to the in)etion point s funtion of ∆ for the e a QHD g a PFT seD whih is lose to the stility limitF he results otined for oth quntities onsidered show seondEorder onvergene of the yp methodF woreoverD the distne from the ontt line position to the in)etion point of the referene solution is SH µmF his result ompres fvorly with the experimentl tehniques used to ompute the ontt ngle flke et al. IWWWD flke 8 hikhmurzev PHHPF e hve demonE strted tht for su0ient resolution of the interfeD our yp model with xvier oundry ondition on the moving sustrte oupled to n imposed ontt ngle is le to reprodue the min stility results of the urtin oting on(gurtion of viu viu et al. PHITF woreoverD this is done with fr fewer points per slip length thn in the (nite element omputtions suggesting higher ury nd roustness of the yp method with height funtionsF

Squalane drop spreading with dynamic contact angle

sn this studyD we use the model hynQ9 desried in vegendre 8 wglio PHISF he slip length λ is (xed to ∆/2 nd l is smll length sle of the order of the nnometerF et eh time stepD the dynmi ontt ngle θ d is omputed nd imposed through the height funtionsF iqution WFIH gives

θ d = g 1 g (θ S ) + Ca CL log L l .
@WFIQA sn prtieD g nd g 1 n e pproximted y (tting polynomil s desried in hupont 8 vegendre PHIH

g(x) = 1 9
x 3 -0.00183985 x 4.5 + 1.845823 × 10 6 x 12.258487 , g 1 (x) = 1 4.33

x 1/3 + 0.0727387 x -0.0515388 x 2 + 0.00341336 x 3 .

e ompre our yp dynmi ngle model with experiments from vvi nd wrE mur vvi 8 wrmur PHHRF hey onsider the spreding of squlne drop on solid sustrteF he drop of volume V = 4πR 0 /3 is initilized on wll with no initil veloityF he yhnesorge numer hrterizing the )ow is

Oh = µ l ρ l σR 0 = 0.21,
nd stti ontt ngle is θ S = 41.5°F e rry out simultions using the sme nnosopi length l = 10 9 m s in vegendre 8 wglio PHISF sn pigure WFIQD we present the omprison etween our simultions nd the experimentl dt for di'erent gridsF he level of re(nement orresponds to the smllest grid size in the qudtree prtitioning of the omputE tionl dominF vevels R to T orresponds to VD IT nd QP grid points per dimeterF por low level of re(nementD the ontt line eventully gets pinned nd the motion is stoppedF et level TD the (nl rdius mthes the experimentl oneF woreoverD in trnsient stteD we n see good greement etween othD vlidting our impleE menttion of the dynmi ngle model in yp frmeworkF sn the next hpterD we present the newly developed in models for moving ontt linesX es shown in ghpter WD the xvier oundry ondition !oupled to onstnt or dynmi ontt ngle! is su0ient to regulrize the solution t the ontt line in some sesF he noEslip prdox is relxed nd the numeril slip in simultions n e used s (tting prmeterF xeverthelessD hevuhelleD tossernd nd leski hevuhelle et al. PHHU pointed out @in eppendix fA tht the xfg fils to suppress ll the singulrities t the ontt lineF sn prtiulrD the pillry pressure diverges t the triple pointF fy studying the interfil slip in porous mediD they showed tht from the frinkmn equtions nd the hry lw in the se of thin oundry lyerD one n derive the following oundry ondition for the tngentil veloity t the wll

u x -λ ∂u x ∂y -a λ 2 ∂ 2 u x ∂y 2 = U, @IHFIA
where u x is the xEomponent of the veloityD λ the slip lengthD a positive oe0ient of order one nd U the wll veloityF fy introduing two seprte slip lengths λ 1 122 Chapter 10. Advances in models for contact lines nd λ 2 D iqution IHFI eomes

u x -λ 1 ∂u x ∂y -λ 2 ∂ 2 u x ∂y 2 = U. @IHFPA
his formultion llows us to independently vry the (rst orderEslip λ 1 nd seondE order slip λ 2 to investigte the e'et of eh termF he numeril implementtion nd vlidtion of the superEslip model in the yp frmework re presented in epE pendix iF 10.1.2 Application to a spreading drop o study the e'et of the superEslip model on the dynmi ontt lineD we onsider drop initilized with θ 0 = 90 °ontt ngle on stti sustrte @U = 0AF e impose n onstnt ngle θ e = 60 °tht will determine the equilirium shpeF he e'et of λ 1 nd λ 2 will only pper in the trnsient stteF he visosities nd densities in oth phses @µ l D µ g D ρ l D ρ g AD nd the surfe tension σ re equl nd set to oneF he initil rdius is R 0 = 0.5 nd ll the simultions re run until (nl time t f = 10 with (xed grid sping ∆ = 0.0078125 orresponding to IPV points per initil dimeterF e onsider four di'erent ses IF λ 1 = 0 nd λ 2 = 0 orresponding to the homogeneous hirihlet oundry ondition u x = 0 @noEslipA t the wllF PF λ 1 = 0.1 nd λ 2 = 0 orresponding to xvier oundry ondition s deE sried in etion WFPFIF QF λ 1 = 0 nd λ 2 = 0.01 where only the seondEorder slip modi(es the spreding motionF RF λ 1 = 0.1 nd λ 2 = 0.01 where λ 2 = λ 2 1 s in iqution IHFIF sn the following (guresD we look t the drop displementD the ontt line speedD the ontt ngle nd the sher stress s funtion of time for the four sesF sn pigure IHFID we plotted the rdius s funtion of timeF he (rst notle e'et is the ft tht ses PDQ nd R overshoot the (nl stti solutionF sn those sesD the ontt line position relxes k to the equilirium one wheres in se I @noEslipAD it lwys remins elow the ltterF his is one of the nturl hrteristis of the slip modelsY llowing the ontt line to slip over the sustrte results in fster motion tht my pss over the equilirium positionF snterestinglyD it ppers tht the seondEorder slip is the leding order termF he ump9 is lrger in se Q @seond order slipA thn in se P @(rst order slipAF woreoverD we n see tht the e'et of oth slips is dded s the ump is mximl in the fourth seF sn ll the ses the ontt line eventully relxes to its equilirium position s expetedF pigure IHFP shows the lol pillry numer @iqution WFIIA omputed in the ontt line ell s the drop is spredingF he ontt line speed is initilly zeroD inreses until mximl vlue efore tending to zero s the stti shpe is rehedF sn ses PD Q nd R the speed eomes negtive fter overshootE ing the equilirium position lerly inditing tht the ontt line motion is reversed to reh the equiliriumF sn the (rst seD howeverD Ca CL is lwys positive or equl to zeroF he highest Ca CL is ttined in se RD inditing ginD the ddition of the e'ets of λ 1 nd λ 2 F sn the noEslip seD we n oserve osilltions of the speed s it is tending to zeroF hese osilltions re relted to grid e'etsF es the ontt line rosses ell oundryD the newly omputed volume frtion is very smllF his indues loss of ury in the omputtion of the surfe tension termF he osilltionsD lthough smllerD re lso present in the slip se wheres they re lmost entirely removed in ses Q nd R @where λ 2 ̸ = 0AF his lst oservtion illustrtes the regulrizing e'et of the seondEorder slipF sn pigure IHFQD we show the vrition of the ontt ngle s funtion of timeF he ngle is omputed t eh time step in the ontt line ellF he (rst notiele thing to oserve is the orreltion etween osilltions in ontt line speed nd ontt ngleF his is in greement with the previous remrks on the loss of ury when the ontt line rosses ell oundryF sndeedD s detiled in etion WFID the ontt ngle diretly in)uenes the lultion of the urvture term therefore modifying the vlue of the veloity when solving the twoEphse xvierEtokes equtionsF e lso note tht the (nl ontt ngle is slightly di'erent from the imposed equilirium ngle θ e = 60 °F here is di'erene etween the imposed9 ngle through height funtions nd the extrted9 ngleF his ∆θ num D tht exists in ll of our ontt line simultionD seems to depend on the )ow hrteristis nd the grid resolutionD nd is still eing investigtedF pinllyD in pigure IHFRD we look t the sher stress ∂u x /∂y t the ontt line s funtion of timeF his lst omprison etween the ses is essentil to understnding the regulrizing e'et of the seondEorder slipF sn the noE slip simultionD the sher stress strongly osilltes nd hnges sign s the drop rehes the equilirium shpeF he hnge in sign indites strong modi(tion of the )ow round the ontt line positionF hese osilltions ! lthough with no hnge in sign ! re still present in se P wheres in ses Q nd R @λ 2 ̸ = 0AD the sher stress is smooth outF ell these onsidertions demonstrte the regulrizing e'et of the seondEorder slipF xonethelessD the superEslip model still needs to e tested on physil systemsF he newly onsidered systems would need to llow @veryA high numer of grid points per slip

λ 1 s λ 2 = λ 2 1 nd λ 1 ≪ 1F 126
u x -λ ∂u x ∂y = U + f x ε σ µ l (cos(θ d ) -cos(θ e )) , @IHFQA
where µ l is the liquid visosityD σ the surfe tensionD θ d the dynmi ontt ngleD θ e the presried ontt ngle nd f x ε ell funtion de(ned s

f x ε = 1 -tanh 2 x ε ε . @IHFRA
xote tht this ell funtion is di'erent thn the one de(ned in etion WFPFIF he width ε is now relted to the di'use interfe modelF he reltion etween the ontt ngle nd the veloity of the ontt line is now Ca CL = f (θ d )F his mens tht the ontt ngle should e solution of the )ow insted of eing presriedD s it ws the se until nowF sn the yp method suh nonEimposition9 of the ngle is not possile t (rst glneF sndeedD even if we were not to impose n ngle through the height funtionsD the homogeneous xeumnn oundry ondition on the volume frtion t the wll

∂c ∂y | y=0 = 0
would still e requiredF his symmetry oundry ondition imposes ! y defult ! 90 °ontt ngle for the interfeF portuntelyD from kinemti reltions on the dynmi ontt lineD the uthors of prike et al. PHIVD prike et al. PHIW derived reltion etween the numeril ontt ngle nd n pprent9 ngleD spei( distne wy from the wllF he numeril ngle n e extrpolted using the following formul

θ ext = θ app + 3 2 ∆ κ 1 + h y sin(θ app ) @IHFSA
where θ ext is the extrpolted ngleD θ app the pprent ngleD ∆ the grid spingD κ the urvture nd h y the (rst order derivtive of the height funtion in the y diretion @norml to the wllAF e shemti of this extrpoltion is presented in 10.2. The generalized Navier boundary condition 127 pigure IHFSF yne the extrpolted ngle is omputedD we impose it through height funtionsD s if we would for regulr onstnt ontt ngle

θ d = θ ext
his method llows us to hve free9 dynmi ontt ngle in the sense tht it is now solution of the )owF e vlidtion of the extrpoltion method is provided in eppendix pF

h 0 h 1 h 2 θ d ⃗ n ⃗ n y x 3 2 ∆ θ app Figure 10
.5: Extrapolation of the contact angle using the apparent angle located 3/2 ∆ away from the wall.

Forced dewetting and lm transition

e study the e'et of the qxfg on the onset of (lm trnsition in fored dewetting setup efkhmi et al. PHIUF e onsider solid sustrte tht is withdrwn from visous liquid poolF he interfe my either sustin sttionry stte menisusD if elow ritil pillry numerD Ca c D or ontinue to move up the sustrte until depositing thin (lm to ritrry heightsF en exmple of stedy stte menisus is shown in pigure IHFTF sn efkhmi et al. PHIUD the uthors used hirihlet oundry ondition for the tngentil veloity nd showed tht there exists n impliit gridEdependent slip ontrolling the trnsition to liquid (lmF sn n ttempt to demonstrte tht the qxfg is gridEindependent ! provided tht the width ε of the ell funtion is well resolved ! we will rry out similr simultions for di'erent widthsD gridsD pillry numers nd equilirium nglesF he solid plte on the left is withdrwn t speed U nd the omputtionl domin is L × L ox where L = 9l c with l c = 1 the pillry lengthF he liquid density is set to ρ l = 1 nd the gs density to ρ g = 0.2F he liquid nd gs visosities re equl nd set to µ l = µ g = U D suh thtD with σ = 1D we hve the reltion Study of the eect of the width ε es (rst stepD in order to investigte the e'et of the width ε nd the rtio ε/∆D we (x the pillry numer to Ca = 0.1 nd n width ε 0 = 0.078125 tht will serve s seF e run the simultion for IV di'erent sesD orresponding to P equilirium nglesD Q widths nd Q grid spingsX

U = √ Ca
θ e = 80 °, 90 °D ε = ε 0 , 2ε 0 , 4ε 0 D ∆ = ε 0 /8, ε 0 /4, ε 0 /2F
he numer of grid points per width ε/∆ will therefore rnge from P to V in the se ε = ε 0 D from R to IT in the se ε = 2ε 0 nd from V to QP in the se ε = 4ε 0 F edditionllyD the slip length is set λ = ε in order to hve unique mirosopi length sle to ontrolF sn the following resultsD we will investigte the ontt line height nd the ontt ngle s funtion of time s well s the sher stress t (nl timeF pigure IHFU shows the normlized ontt line height (hh 0 )/l c s funtion of timeF sn ll of the sesD the ontt line rehes stedy stte position implying tht Ca = 0.1 < Ca c F sn the θ e = 80 °ses the (nl height is higher thn for the θ e = 90 °ses whih is expeted s more hydrophoi sustrte @θ e = 80 °A will fore the ontt line towrds the top of the dominF hen deresing εD the e'et of the unompensted oung stress is further lolized round the ontt line nd the slip length is deresedF he (nl position of the ontt line is mximl for the smllest εF vooking t the rtio ε/∆D we n see ler onvergene ptternF he height t highest resolution onverges towrds ertin vlue omprised etween the P other resolutionsF roweverD looking the ε = ε 0 with P grid points per width ses @dotted lk linesAD we n oserve signi(nt drop in the (nl position with respet to the R nd V grid points per width sesF his lerly shows tht P points per width is not enough to resolve the smllest length sleF his ssertion is in ordne with the previous study on the xvier oundry ondition @etion WFQFIAF pigure IHFV shows the ontt ngle s funtion of time in logrithmi sleF xote tht this ngle is not the one we impose through the extrpoltion method @iqution IHFSAF st is the extrted ngle t the ontt line ell @similrly to pigure IHFQAF gonsider (rst the θ e = 80 °sesF he ngle strts from 90 °D rosses the the equilirium ngle nd eventully settles towrds (nite vlueF hen the dynmi ngle rosses this thresholdD the sign of the unompensted oung stress hngesD s cos(θ d ) -cos(θ e ) is negtive when θ d > θ e nd positive when θ d < θ e F his trnsition ppers reltively fstD mening tht the e'et the stress pushes further the ontt line in the diE retion of the moving plteF sn the 80 °sesD the ngle is lwys elow the threshold exept for ε = ε 0 with P grid points per widthF sn tht spei( seD one n oserve signi(nt osilltions of the ngle tht re in ordne with the previous remrks on the (nl heightF sndeedD the (nl vlues of the ngles re lmost the sme regrdless of the numer points per width mening tht these osilltionsD only ppering in trnsient stteD still in)uene the stedy stte solutionF sf we now onsider the di'erenes indued y modifying the solute vlue of εD we n oserve tht the (nl ngle seems to onverge towrds ertin vlue s the e'et of the oung stress is further lolizedF pinllyD in pigure IHFW we look t the sher stress t (nl time s funtion of the y position of the wllF e n lerly see the e'et of the redution of the widthF por lrge εD the sher stressD while still rehing mximum t the ontt line positionD is further spred out thn for smll εF yne interesting oservtion is the nonEmonotoni vrition of the stress in the smllest ε sesF here is n in)exion ppering t ertin distne wy from the ontt lineF his e'et still needs to e investigtedF 130 Chapter 10. Advances in models for contact lines prom these resultsD we n drw two onlusions IF e smll vlue of ε with respet to the mrosopi length sle ! the pillry length in this se ! ppers s the nturl wy of treting the qxfg in the shrp interfe limitF his is understndle s ε is relted to the interfe width in the di'use interfe modelF PF ε needs to e wellEresolved in order to otin onverged resultsF es expetedD the smllest length sle needs to e muh lrger thn the grid sping ∆ ≪ εF fering these onlusions in mindD we n now ompre our qxfg simultions to the noEslip simultions in order to investigte the gridEdependene ! or independene ! of the modelsF

Comparison between the no-slip and GNBC models e onsider the sme setupD this time with Ca = 0.03D two di'erent equilirium ngles θ e = 60 °, 90 °nd (xed width ε = 0.078125F he grids onsidered re ∆ = ε/4 = 0.01953125 @orse gridAD ∆ = ε/8 = 0.009765625 @medium gridA nd ∆ = ε/16 = 0.004882813 @(ne gridAFF pigure IHFIH shows the results otined with the noEslip model @u x = U AF sn the 90 °sesD the liquid pool rehes stedy stte nd the trnsition to (lm does not ourF here isD s expetedD gridEdependene of the modelD s the (nl height is inresed when ∆ is deresedF his dependene is mpli(ed in the 60 °sesF he trnsition now ours in the (ne grid se wheres it does not in the other sesF es the grid is deresedD the numeril slip ppering in the noEslip model is deresed s wellD therefore tending to true noEslip9 motion where the trnsition would exist for smll vlues of CaF hese oservtions re in greement with the ones of efkhmi et al. PHIUF pigure IHFIH shows the results otined with our qxfg modelF st is lerD from the results thtD in this prtiulr seD tht the qxfg model exhiits true gridE independent ehviorF sn the θ e = 60 °se there is no trnsition to (lm formtion even with (ne gridF prom this omprison nd the previous study on the e'et of the widthD we n onlude thtD in this prtiulr se with Ca = 0.03D our qxfg model is n imE provement with respet to the noEslip modelF roweverD we need to further investigte this ehvior y simulting wide rnge of pillry numers nd looking t the e'et of the unompensted oung stress independently of the slip @λ = 0AF sn the hseEpield @pA methodD the ghnErillird equtions re formulted from the thermodynmis of n immisile twoEomponent mixture nd phse funtion is used to represent the moving interfeF he veloity (eld n stisfy noEslip ondition t the ontt lineD nd the ontt line moves due the di'usive mss trnsport tht is present in twoEomponent system ue 8 peng PHIIF hile the p method hs the dvntge of needing less empiril (tting thn the yp methodD it is onsiderly more expensiveF sn this studyD we show how to otin results tht re essentilly equivlent to results from the p methodD y doing muh heper yp simultionsF e introdue omputtionl oundry t spei(ed short distne from the wll in the yp simultionF here we pply xfg for tngentil veloityD nd n dynmi ontt ngle reltion tht ounts for the interfe urvture t the wllF he sis for the model is theoretil reltion etween totl di'usive mss trnsport ross the interfe nd urvture t the wllF e ompre full phse (eld simultions with yp simultions using these oundry onditionsF he input prmeters for the yp simultions re the sme s for the full p solution nd the yp solutions reprodue the p solutionsD t muh lower ostF e strt e tking the dimensionless inompressile twoEphse xvierEtokes equtions in the hseEpield pproximtion @de(ned in ghpter UAF es referene quntities we hoose U D LD L/U nd µ l U/L for veloityD lengthD time nd pressureD with L denoting mrosopi lengthD for exmple the droplet rdiusF he dimensionless numers tht pper re the eynolds numer

Re = ρ l U L µ l , the elet numer Pe = U L α ,
where α = M σ/ε is the mss di'usion oe0ientF he ghn numer ! or dimenE sionless interfe width ! is de(ned s Cn = ε L .

Ca F = µ F U/σ is the line frition pillry numerD sed on the line frition pE rmeter µ F F sn the rest of this studyD we onsider the line frition oe0ient to e zero @µ F = 0AD in order to enfore the dynmi ontt ngle to the equilirium oneF he dimensionless equtions therefore red

ρ l ∂u ∂t + u • ∇u = - 1 Re ∇P + 1 Re ∇ 2 u + 1 Re Cn Ca ϕ∇C, @IHFTA ∂C ∂t + u • ∇C = 1 Pe ∇ 2 ϕ, @IHFUA with ϕ = -Cn 2 ∇ 2 C + Ψ ′ (C) , @IHFVA nd Ψ ′ (C) = (C 2 -1)C.
he wetting oundry onditionD with µ F = 0D is now

3 2 √ 2 Cn ∇C • n = cos(θ e )g ′ (C), with g(C) = 2 + 3C -C 3 4 .
essume length sle δ for di'usionD whih is lso oupled to the )ow (eldF essume visous lne in the momentum equtionD nd tht the veloity is properly sled @dimensionless |u| ∼ 1)F istimtes of the relevnt terms iqutions IHFT to IHFV give

C = C 0 + Cn •C 1 + O(Cn 2 ), ϕ = ϕ 0 + Cn •ϕ 1 + O(Cn 2 ).
he prolems for the zeroth nd (rst order re now

ϕ 0 = - ∂ 2 C 0 ∂p 2 + Ψ ′ (C 0 ) , @IHFIIA nd ϕ 1 = - ∂ 2 C 1 ∂p 2 - 1 R ∂C 0 ∂p + C 1 Ψ ′′ (C 0 ) . @IHFIPA
he zeroth order eqution IHFII gives the equilirium solution s C 0 = -tanh(p/ √ 2)D with ϕ 0 = 0F wultiplying the (rst order eqution IHFIP y ∂C 0 /∂p nd prtilly integrtingD solvility ondition is otined for the (rstEorder vlue of the hemil potentil

ˆϕ1 C 0,p dp = -ˆ(C 0,p C 1,pp + C 1 (Ψ ′ (C 0 )) ,p )dp - 1 R ˆ(C 0,p ) 2 dp, @IHFIQA
where , p denotes derivtive with respet to pF he hemil potentil is ontinuous t the interfeD so the (rst order hemil potentil ϕ 1 is onstnt in equE tion IHFIPF fy prtilly integrting the (rst term on the right hnd side tht is seen to e zeroF sing tht ˆϕ1 C 0,p dp = -2,

nd ˆ(C 0,p ) 2 dp = 2 √ 2 3
, the result gives the vlue of the hemil potentil t the interfe

ϕ = Cn ϕ 1 = 1 2 2 √ 2 3 Cn R .
@IHFIRA e will now pply the sme perturtion pproximtion to iqution IHFUF eginD the zeroth order gives the equilirium onentrtion pro(leD ut the (rst order eE omes

-u • Pe •C 0,p = (-C 1,pp + C 1 Ψ ′′ (C 0 )) ,pp - 1 R C 0,ppp . @IHFISA
sn order to derive reltion etween the mss )ux on the left hnd side nd the urvtureD we multiply this eqution y

F = ˆp 0 C 0 dp. @IHFITA
fy suessive prtil integrtionsD the (rst term in the right hnd side is shown to vnishD nd the end result is ˆu

• Pe •F F pp dp = 1 R ˆF F pppp dp, @IHFIUA with ˆF F pp dp = 2 (2)(1 -ln 2) = 0.8679, nd ˆF F pppp dp = 2 √ 2 3 = 0.9428.
sn iqution IHFIUD ftor δD the p distne to the stgntion pointD is introdued on the left hnd side to ount for n integrtion in the vertil diretionF he veloity pro(le on the left hnd side goes to zero eyond δD ut the right hnd side omes from the di'usion (eld round the ontt line whih n hve ontriutions from further wyF e now onsider the xvierEtokes equtions @iqution IHFTA to ouple this to the veloity (eldF essuming lne etween the visous term nd the pillry termD we estimte the relevnt terms in the xvierEtokes eqution s 1 Re

a u δ ∼ 1 R 2 √ 2 3 1 Ca Cn Re .
rere dimensionless prmeter a is introdued to ount for the visous stress τ s τ = au δ .

utting everything together yields the following reltionsX the hemil potentil ϕ is relted to the rdius of urvture R t equilirium y

ϕ = 1 2 2 √ 2 3 Cn R , @IHFIVA
nd the lol urvture is funtion of the )ow prmeters 10.3.2 Phase-Field simulations of dynamic wetting by the Volumeof-Fluid method e onsider drop spreding over noEslip motionless sustrte ! loted t y = 0 ! where the equilirium ngle θ e governs the (nl shpe of the dropF sing the reltions desried oveD we ompre full p simultions ! in the sense tht the full domin is solved ! with yp simultions using the toy model pproximtionF sn the toy modelD the ottom oundry is now omputtionl oundry tht is loted distne +δ wy from the wllF he domin omprised etween the wll y = 0 nd the omputtionl oundry y = δ is not simultedF hereD we pply the dynmi ngle reltion @iqution IHFPHAD where the pprent ngle θ app is now dynmi ontt ngle θ d

R = 2 √ 2 3 1 √ 0.
θ d = cos 1 cos θ e + 3 2 √ 2 Ca CL a , @IHFPIA
with Ca CL the ontt line pillry numerF woreoverD in order to simulte true noEslip ondition t the wllD we use δ s the slip length for the tngentil veloity @gin loted t the omputtionl oundryA

u x | y=δ -δ ∂u x ∂y | y=δ = 0. @IHFPPA
hen onsidering spreding drop with this modelD we need to relx the imperE meility ondition in the xfg y llowing mss )ux through the omputtionl oundryF his sinking veloity v sink is relted to the mount of mss tht we lose ! in the se of θ e < θ 0 ! or gin ! in the se of θ e > θ 0 F he oundry ondition for the norml omponent of the veloity t the wll is now

u y | y=δ = v sink .
@IHFPQA e vlidtion of the numeril implementtion of v sink is provided in eppendix qF pigure IHFII shows the omputtionl domins of oth p nd yp setups with the respetive oundry onditionsF he initil rdius of the drop is R 0 = 0.5 with n initil ontt ngle θ 0 = 90 °F he equilirium ngle is set to θ e = 70 °F he visosities nd densities rtio etween oth phses re xote tht the hoie of elet numer grees with the shrp interfe limit pproxE imtion of ue nd peng ue 8 peng PHII

μ = µ g µ l = 10 2 , ρ = ρ g ρ l = 10 2 . etting fg 3 2 √ 2 Cn ∇C • n = cos(θ e )g ′ (C) C = 1 C = 1
Pe = 1 < 16 √ μ Ca Cn = 3.392.
he geometril ftor ounting for the visous )ow is set to a = 3F his vlue ws found y (tting the rdius of urvture nd the right hnd side of iqution IHFIWF he p simultions re rried out using preepiwExEgrD sript ompiltion for preepiw reht PHIP to solve the oupled system of xvierEtokes ghnErillird equtions developedD y F v© is nd F fgheriF he grid size used in the p simultion is ∆ PF = 0.0625 ε nd the one in the yp simultion is ∆ VOF = 0.577 δF he results otined re summrized in the following (guresF sn pigure IHFISD we present tle of the di'erent models with the foremenE tioned tegoriztionsF he olor ode is s followsX sn greenD the models implemented nd tested in the yp frmework tht re presented in ghpter W nd ghpter IHF he toy model is speil type of slipE dynmi ngle omintion were the the omputtionl oundry is lifted o' the solid sustrte nd neessittes speil tretment of the norml veloity omponent in the oundry ondition @eppendix qAF sn lueD the models the were implemented nd tested ut not presented in this disserttionF sn prtiulrD the qxfg pth method pulln PHIW ws the (rst itertion of the qxfg where we modi(ed the veloities round the ontt line without hnging the xvier oundry onditionF sn tht seD onstnt ontt ngle ws usedF sn redD the models tht re yet to e testedF sn greyD the omintions of oundry onditions nd ontt ngles tht hve no mening9 in our numeril frmeworkF por exmpleD in the se of qxfg with dynmi ontt ngleD the ft thtD on one hnd the ontt line speed Ca CL is funtion of the dynmi ngle θ d D ut on the other hnd θ d is itself funtion of the ontt line speed Ca CL lerly leds to prdoxF he free dynmi ontt ngle in the se of slip or superEslip model leds to n illEposed prolemD where the ondition on the ontt ngle is lking nd does not enter the oundry ondition of the momentum equtionF

θ d = f (θe, Ca, L, l) ux = f (λ 1 , λ 2 , U ) θ d = f (θe, Ca, L, l) ux = f (λ 1 , λ 2 , U ) etion IHFIFI θ = θe lidtion eppendix i preding droplet etion IHFIFP ux = f (ε, θ d , θe, U ) θ = θe th method pulln PHIW ux = f (ε, θ d , θe, U ) etion IHFPFI θ d = f (θapp, κ, ∆, h)

Conclusion

sn the (rst prt of this disserttionD we proposed level set method oupled with novel gut gell pproh to tkle twoEphse tefn prolems s well s n djointE sed optimiztion proedureF he key fetures of the level set relted lgorithms reD @iA n impliitEexpliit sheme to solve the level set dvetion llowing us to relx the usul gpv onditionD @iiA the highEorder tohnsenEgolell method used to ompute the norml grdient ross the interfeD nd @iiiA the suEell resolution reinitiliztion proedure to retin the signed distne funtion property s the interfe movesF woreoverD we use the gut gell method oupled with grnkE xiolson time integrtor tht llows us to solve the twoEphse prolem for ny given geometryF he method ws then extended to solve the xvierEtokes equtions in the )uid phseD nd ompres fvorly with previous omputtions on the onset of yleighEfénrd instilitiesF his gut gell method is urrently eing extended to higher dimensionsF sn the ontext of lssil twoEphse tefn prolemsD the djointEsed optiE miztion proedure is shown to e roust lgorithm to ontrol the shpe of melting or solidi(tion frontD even in the presene of dendriti instilities nd nisotropi e'etsF he omprison with derivtiveEfree methods yields fvorle resultsD s the grdientEsed method onverge fsterF xeverthelessD the (nl vlue of the ost funtionl isD s expetedD lower in the derivtiveEfree sesF es future workD one ould onsider to employ hyrid optimiztion methodsD suh s the qE hgys @qrdientEenhned hnmi gyordinte serh using esponse urfe modelsAF he derivtion of the ontinuous djoint in the presene of )ow in the liquid phse ould lso e diretion to followF sn tht seD the use of n inomplete djoint ! were the veloity (eld is known9 quntity in the djoint prolem ! might e requiredF sn this prtiulr seD to void ny sustntil simpli(tionsD one ould lso onsider the disretizeEthenEdi'erentite9 pproh were disreteE djointEsed optimiztion lgorithm ! s opposed to ontinuous one ! would e uiltF sn the seond prt of this work we presented severl ontt line models in olumeEyfEpluid frmeworkF hese models were vlidted nd tested on physil setupsF hesriing the physis of the dynmi ontt line with shrpEinterfe method is still hllengeF he models yield di'erent results depending on the physil prolems onsidered nd onsensus on the est9 pproh to dopt is still detedF yur work is n ttempt to enrih the (eld with vrious new numeril methods for ontt linesF sn prtiulrD the qxfg sed on the yp methodD were the ngle is funE tion of the )owD presents ttrtive gridEindependent propertiesD provided tht the smllest length sle is wellEresolvedF yn the other hndD the toy model model is n interesting wy of modeling evportion ! mehnism driving of the motion of the Conclusion ontt line ! through simple oundry onditions in the yp methodF his model still needs to e tested on wider rnge of prolemsF es (nl remrkD one ould implement these models in our level set methodD thus modeling phseEhnge phenomen in the presene of ontt linesF en interesting dvntge would e to use the lredy existing djointEsed ! or derivtiveEfree ! optimiztion proedures nd test them on se where the ontrol vrile isD for exmpleD the ontt ngleF Part III Appendices his proess is repeted for the omponents of the o'Edigonl elements of the strinErte tensorD de(ned in the ontinuous se s

s αβ = 1 2 ∂u α ∂x β + ∂u β ∂x α , α ̸ = β,
nd in the disrete se s the nodeEentered (eld

∀ (α, β) ∈ {x, y} 2 , α ̸ = β, strain αβ (U, D) = 1 2W αβ   δB αβ U α δξ β + δ A αβ β -B αβ D α δξ β - δA αβ δξ β D α β   + 1 2W βα δB βα U β δξ α + δ A βα α -B βα D β δξ α - δA βα δξ α D β α
. @fFRA st should (nlly e noted tht the ltter formul @iqution fFRA is lso vlid in the digonl se @α = βAD in whih se it simply redues to iqution fFQF Viscous transport term rior to proeeding with the disretiztion of the visE ous trnsport termD it should (rst e noted thtD in the se where the seond rgument @NA of the divergene opertor @iqution QFIQ summed over αA mthes the (rst rgument @QAD iqution QFIQ my e simpli(ed using the identities preE sented y worinishi PHIH s

div (Q, Q) = β B β δQ β δξ β . @
(u • ∇) u = ∇ • (u ⊗ u) -∇ • u = ∇ • (u ⊗ u) , whih in disrete form n e written s conv α U, U † , D, D † = β      δA β U β α U † α β δξ β +    δ B β β -A β D β α δξ β - δB β δξ β D β β α    U † α + D † α 2     
. @fFTA his multiliner opertor is typilly evluted t U † = U nd D † = D ut the distintion might er signi(neD in the ontext of irt lineriztion for exmE ple where distintion pplies etween U whih is typilly frozen wheres U † is updtedF his disretiztion n e onsidered s the generliztion of the enE tered sheme to the ut ell methodD whih n e demonstrted s followsF sn the ontinuous seD

∀ (α, β} ∈ {x, y} 2 , u α ∂u β u α ∂x β = ∂u β u 2 α /2 ∂x β + u 2 α 2 ∂u β ∂x β , @fFUA
whihD upon summtion over αD yields similr eqution for the spei( kineti energy k ≡ ∥u∥ 2 /2D ultimtely onserved in the invisid limitF he proposed disE retiztion of the onvetive trnsport term @iqution fFTA preserves this property t the disrete levelF sing the identities presented y worinishi PHIHD it n e e shown tht

∀ (α, β) ∈ {x, y} 2 , U † α δA β U β α U † α β δξ β = δA β U β α U † α U † α β /2 δξ β + U †2 α 2 δA β U β α δξ β , @fFVA
where • denotes the permnent produt ϕψ x i+ 1 /2,j = ϕ i+1,j ψ ij + ψ i+1,j ϕ ij 2 , @fFWA lso introdued y worinishi PHIH nd esily extended to other dimensions nd rrngements s previously done for di'erentition nd interpoltionF iqution fFVD together with the ontinuity opertor @iqution fFIAD n e used to show tht ∀α ∈ {x, y}

U † α conv α U, U † , D, D † = β      δA β U β α U † α U † α β /2 δξ β +    δ B β β -A β D β α δξ β - δB β δξ β D β β α    U † α D † α 2      + U †2 α 2 cont (U, D) α . @fFIHA
his identity n e interpolted in eh diretion αD nd summed over αD to ultiE mtely stte the proposed disretiztion @iqution fFTA onserves kineti energyD in Flower.jl package description sn this setionD we detil the tuli ode Flower.jl used in the (rst prt of the disserttionF he ode rhiteture n e summrized s followsX he numeril nd physil prmeters ! de(ned y the user when initilizing prolem ! the meshesD disrete opertors nd (elds re ontined in di'erent mutle dt strutures tht enter the min solving funtionsF wo min funtionsD with options determining the type of prolem to onsiderD re de(nedF he (rst one ! run_forward(...) ! solves the forwrd tefn prolem @elgorithm RA nd the seond one ! run_backward(...) ! solves the djoint tefn prolem @elgorithm TAF e now present typil tuli (le tht uses Flower.jl to solveD for exmpleD the yleighEfénrd onvetion se @etion SFRAF e strt y de(ning the numeril nd physil vriles ontined in the num strutureX § ¤ he simultions re run until (nl time t/t µ = 100 where t µ = D 2 /µ is the visous time sleF e look t the evolution in time of the mximl dimensionless veloity sled with the pillry numer

Ca max = µ|u| ∞ σ .
sn pigures hFP to hFRD the results will e ompred two y two @noEontt or onttA the show the e'et of the oundry ondition on the spurious urrentsF pigure hFP shows the spurious urrents in the se where no shift is ppliedF por ll the omintion of prmetersD the mximum veloity eventully onverges to zero within mhine preisionF he (rst ler oservtion is tht the dmping time in whih the spurious urrents re killed is strongly dependent on the vple numerF es expetedD s La inreses @lines from lk to purpleAD the time to reh mhine preision deresesF sn pigure hFP @noEonttAD for v a PHHD we oserve ler onvergene s the grid size is deresedF sn pigure hFP @onttAD the time to dmping is gretly mpli(ed in the D/∆ = 16 sesF his might ome from numeril errors in the height funtions due to n insu0ient resolutionF pigure hFQ shows the spurious urrents in the se where ∆/4 shift is ppliedF he previously grid ligned volume frtionsD where c = 1D re now displedF sn the noEontt seD we see similr results thn in the previous (gureF roweverD the grid onvergene is 'eted s the highest resolution ses re not the ones with smller dmping timesF sn pigure hFQD we see spurious urrents tht do not vnish s time proeedsF por v a PHHD t the highest resolution D/∆ = 64D the mximl veloity does not onverge towrds mhine uryF his is n exmple of prolemti spurious urrentsF sn prtil simultionsD the )ow round the ontt line ould e onsiderly 'eted leding to n umultion of numeril errorsF sn pigure hFR @∆/2 shiftAD howeverD we notie the opposite situtionF he ontt with the sustrte seems to hve regulrizing e'et on the solution s ll the veloities tend to mhine preisionF sn the noEontt seD ll the simultions with grid D/∆ = 32 exhiit nonEvnishing spurious urrents regrdless of the vple numerF his ounterEintuitive se might e n exmple of speil se9 where the omintion of these spei( shiftD grid size nd initil position of the drop leds to numeril errors in the omputtion of the urvture nd therefore in the surfe tension termF where a is positive onstntF sn this studyD we only onsider the dvetion eqution of the olor funtion

∂ t c + ∇ • (cu) = 0.
he presried inompressile veloity (eld will indue osilltions of the drop in oth vertil nd horizontl diretionsF he ngle formed t the ontt line will e 'eted y this motion nd vry in timeF e will ompre the oserved nuE meril ontt ngle with the nlytil one θ exact D given y the following relE tion prike PHPH

θ exact = π 2 + tan 1 -1 tan θ 0 e 2aS , with S = sin(πt) π .
e rry out two sets of simultions in order to ssess the methodF sn the (rst oneD the ontt ngle is onstnt nd set to 90 °wheres in the seond oneD we impose the extrpolted ngleF he (nl time is set to t f = 2 nd the grid sizes onsidered rnge from IT to IPV points per dimeterF pigure pFI summrizes the results otinedF he onvergene of the ontt ngle with respet to the nlytil solution is only otined in the θ d = θ ext se thus vlidting the extrpoltion methodF

  PHIQD poures et al. PHIRD in et al. PHIRF edjointEsed methods hve lso een employed for the purE pose of sensitivity nlysis or ontrol in )ows in the presene of lrge grdiE ents yu 8 tmeson PHIID frmn et al. PHISD vemke et al. PHIWD pikl et al. PHPH ! for exmple in )mes or interfes ! showing gret promiseD nd therefore re dopted here to rry out the optimiztion proedureF

Figure 1 . 1 :

 11 Figure 1.1: Schematic of the liquid and solid subdomains considered in the Stefan problem.

  Figure 1.2: Photography of a crystal growing in an under-cooled liquid bath exhibiting a pattern formation with primary and secondary branches [Fujioka 1978].

Figure 1 . 3 :

 13 Figure 1.3: View of the level set signed distance function ϕ. The interface corresponds to the interface of ϕ and the 0-level plane.

Figure 1 . 4 :

 14 Figure 1.4: Schematic of the two-phase Stefan problem in the level set framework.

Figure 1 . 5 :

 15 Figure 1.5: Photography of an atmospheric surface layer [Pandey et al. 2018].

  Figure 1.6: Schematic of the melting boundary problem.

with onstnts β 1

 1 Figure 2.1: Schematic of the computation of the second term of the cost functional. The values of ϕ f and ϕ d are computed in the mixed cells (black dots) determined by the 0-level set of the nal level set function. The dotted lines represent the -2,-1,1,2 level sets of both functions.

  Figure 3.1: Example of an exact and approximate representation of a Cartesian cell cut by the interface.

Figure 3 . 4 :

 34 Figure 3.4: Intersection of the uid domain with Cartesian elements.

Figure 3

 3 Figure 3.8: Convergence study of the Cut Cell method coupled with a Crank-Nicolson scheme when solving the heat equation inside a stationary circle with a Dirichlet boundary condition T D = 1 imposed at the interface. The top gures show the position of the interface in red and the temperature eld at nal time t f = 0.03125 for N = 16, 32, 64, 128. The middle gure show the normalized error in temperature eld with respect to the reference solution taken for N = 256. The bottom gure shows the convergence rate of the method in mixed cells, full cells and in all cells.

Figure 3

 3 Figure 3.9: Convergence study of the Cut Cell method coupled with a Crank-Nicolson scheme when solving the heat equation outside a stationary circle with a Dirichlet boundary condition T D = 1 imposed at the interface and insulated boundary conditions. The top gures show the position of the interface in red and the temperature eld at nal time t f = 0.03125 for N = 16, 32, 64, 128. The middle gure show the normalized error in temperature eld with respect to the reference solution taken for N = 256. The bottom gure shows the convergence rate of the method in mixed cells, full cells and in all cells.

Figure 3 .

 3 Figure 3.10: Convergence study of the Cut Cell method coupled with a Crank-Nicolson scheme when solving the heat equation inside a stationary square with a Dirichlet boundary condition T D = 1 imposed at the interface. The top gures show the position of the interface in red and the temperature eld at nal time t f = 0.03125 for N = 16, 32, 64, 128. The middle gure show the normalized error in temperature eld with respect to the reference solution taken for N = 256. The bottom gure shows the convergence rate of the method in mixed cells, full cells and in all cells.

Figure 3 .

 3 Figure 3.11: Convergence study of the Cut Cell method coupled with a Crank-Nicolson scheme when solving the heat equation inside a stationary crystal with a Dirichlet boundary condition T D = ε κ κ imposed at the interface. The top gures show the position of the interface in red and the temperature eld at nal time t f = 0.0078125 for N = 32, 64, 128, 256. The middle gure show the normalized error in temperature eld with respect to the reference solution taken for N = 512. The bottom gure shows the convergence rate of the method in mixed cells, full cells and in all cells.

Figure 4 . 1 :

 41 Figure 4.1: Schematic of gradient calculation. The Dirichlet value T d is imposed at the interface centroid in the partial cell and the temperatures T A and T B are determined via a quadratic interpolation from the neighboring 3 points in the vertical direction (dotted lines).

Figure 4

 4 Figure 4.3: Velocity elds of a crystal-shape geometry. The red curves represent the interface location and the velocity elds correspond, from left to right, to narrow bands widths of 0, 6 and 12.

Figure 4

 4 Figure 4.4: Schematic of the diamond-cell strategy for the calculation of the Hamiltonian.The red arrow corresponds to the forward coecient and the blue arrow to the backward coecient.

Figure 4

 4 Figure 4.5: The top gures represent the normalized error eld for dierent CFL = 1, 4, 16 and xed N = 64. The red (white, respectively) curve represents the initial (nal, respectively) 0-level set. In the bottom gure the error in L 2 norm is plotted for dierent resolutions and CFL numbers. The vertical blue lines correspond to a xed number of points for varying CFL numbers.

  Figure 4.6: The -2 to 1 level sets of the function ϕ are shown with a 0.1 step. The red curve represents the 0-level set. After 90 iterations the level set function is a true signed distance function.

Figure 4 . 7 :

 47 Figure 4.7: Example of empty cells becoming partial cells from the point of view of the solid phase. The interface at time t n moves in the normal direction with speed V . At time t n+1 , the newly initialized value T new located at the partial cell centroid is extrapolated from T A and T B .

  .jl an in-house Julia package . . . . . . . . . . . . . 67 5.2 Analytical solutions of two-dimensional problems . . . . 68 5.2.1 Planar interface . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2.2 Growing Frank's spheres . . . . . . . . . . . . . . . . . . . 69 5.3 Crystal growth in an under-cooled liquid bath . . . . . . 72 5.3.1 Surface tension eects . . . . . . . . . . . . . . . . . . . . 72 5.3.2 Grid renement eects . . . . . . . . . . . . . . . . . . . . 76 5.3.3 Eect of the molecular kinetic coecient . . . . . . . . . . 77 5.3.4 Anisotropy eects . . . . . . . . . . . . . . . . . . . . . . 78 5.4 Rayleigh-Bénard convection with a melting boundary . 81 5.1 Flower.jl an in-house Julia package sn this ghpterD we vlidte our numeril sheme with respet to nlytil soE lutions of the twoEphse tefn prolemF woreoverD we present numeril results of growing rystl in n underEooled liquid th nd yleighEfénrd onvetion with melting oundryF he ode is written in tuli feznson et al. PHIU high level sienti( progrmE ming lngugeF he key fetures of the tuli lnguge reX Just-In-Time compilationX the ode is ompiled during exeution of progrmD t run timeF Multiple dispatchX dispth proess to hoose whih of funtion9s methods to ll sed on the numer of rguments givenD nd on the types of ll of the funtion9s rgumentsF

5. 2 .

 2 Figure 5.1: Convergence of the temperature eld in the planar interface case. L 2 and L ∞ norms of the error in temperature in all cells as a function of the number of points per dimension.

Figure 5 . 2 :

 52 Figure 5.2: Convergence of the radius towards the analytical solution as a function of time for the growing Frank sphere.

Figure 5

 5 Figure 5.3: Positions of the interface at times t = 1, 1.5, 2 for the dierent grids. The dotted lines correspond to the initial radius R(t = 1) = 1.56 and the analytical radius at nal time R(t = 2) = 2.206. The heat map corresponds to the normalized error in temperature eld for the highest resolution N = 128.

Figure 5

 5 Figure 5.4: Convergence of the temperature eld in the growing Frank sphere case. L 2 norm of the error in temperature for mixed, full and all cells as a function of the number of points per dimension.

Figure 5 . 5 :

 55 Figure 5.5: Example of numerical instabilities in crystal growth due to the absence of surface tension eects. The initial condition is a solid crystal at temperature T = 0 surrounded by an under-cooled temperature eld with T ∞ = -0.5. The color map represents the nal temperature eld.

Figure 5 . 9 :

 59 Figure 5.9: Eect of the molecular kinematic coecient on the dendrite growing tips. The initial solid crystal at a temperature T = 0 is surrounded by an under-cooled temperature eld with T ∞ = -0.5. The color map represents the nal temperature eld. The nal time is t f = 0.4 and the interface is plotted with a time step of 0.02.

Figure 5 .

 5 Figure5.10: Anisotropy eects on the crystal growth for a prescribed angle of symmetry α 0 = π/2. The initial condition is a six-fold solid at temperature T = 0 surrounded by an under-cooled temperature eld with T ∞ = -0.8. The color map represent the nal temperature eld. The nal time is t f = 0.09 and the interface is plotted with a time step of 0.0045.

Figure 5 .

 5 Figure 5.11: Anisotropy eects on the crystal growth for a prescribed angle of symmetry α 0 = π/4. The initial condition is a six-fold solid at temperature T = 0 surrounded by an under-cooled temperature eld with T ∞ = -0.8. The color map represent the nal temperature eld. The nal time is t f = 0.09 and the interface is plotted with a time step of 0.0045.

Figure 5 .

 5 Figure 5.12: Eective Rayleigh number in logarithmic scale as a function of time for the indicated global Rayleigh numbers. The dotted black line indicates the critical Rayleigh number.

Figure 6 . 1 :

 61 Figure 6.1: Example of convex cost functional as a function of two parameters a 1 and b 1 .

Figure 6 . 3 :

 63 Figure 6.3: Normalized cost functional as a function of the iterations of the L-BFGS optimization procedure for the melting circle case.

  Figure 6.4: Iterations 0, 2, 3 and 9 of the optimization procedure for the Mullins-Sekerka test case. The blue curve represents the desired shape and the red one the nal position of the interface at a given iteration. The nal time is t f = 0.5 and the interface is plotted with a time step of 0.05. The color map corresponds to the error in temperature eld. The inset shows the actuator w at a given iteration (red) and the desired one (blue).

Figure 6 . 7 :

 67 Figure 6.7: Normalized cost functional as a function of the iterations of the L-BFGS optimization procedure for the growing crystals case.

  es future workD one ould solve optimiztion ses where the )ow in the liquid phse is onsideredF he ontrol of the onset of the yleighEfénrd instility through the glol yleigh numer ould e one the golsF he derivtion of the ontinuous djoint in tht se might required some simpli(tionsF he )ow veloity u would only e solved in the forwrd prolemD leding to n -phase ows in the one-uid formulation . . . . . . . 95 7.2 Phase-Field extensions . . . . . . . . . . . . . . . . . . . . 97

Figure 7 . 1 :

 71 Figure 7.1: Example of the discrete representation of the density and viscosity in the one-uid formulation.

Contents 8 . 1 Figure 8

 818 Figure8.1: Schematic of the static contact line. The contact line formed at the intersection of a gas-liquid and a solid boundary is element of many natural and several technological processes. In the static case, the equilibrium angle is determined by Young's law.

Figure 9 . 1 :

 91 Figure 9.1: Construction of the 2D height-functions near the contact line [Afkhami & Bussmann 2008].

  Figure 9.2: Example of the quadtree discretization and the corresponding tree representation.

Figure 9 . 3 :

 93 Figure 9.3: Example of geometrical ux estimation.

  PHPPF he xfg orresponding to slip model in Ph on(gurtion with the sustrte t y = 0 n e written s follows u x | y=0 -λ ∂u x ∂y | y=0 = U, u y | y=0 = 0, @WFUA with u x | y=0 nd u y | y=0 the x nd y omponent of the veloity t the solid oundryD λ the slip lengthD nd U the presried veloity of the moving sustrteF x y

Figure 9 . 4 :

 94 Figure 9.4: Left side : velocities proles at the solid interface for the no-slip and slip boundary conditions. Right side : 3x3 stencils with the ghost boundary layer used to impose the boundary condition.

  110

  Figure 9.6: Time series of a curtain coating simulation for Re = 12, Ca = 1.33. The dimensionless time t ⋆ is scaled with the viscous time scale t µ = µ l /ρ l h 2 c .

Figure 9

 9 Figure 9.7: Example of resolution of the contact line for λ = 10 µm with 64 grid spacing per slip length. The color map corresponds to the x-component of the velocity eld. Top: foot of the curtain. Bottom: from left to right, successive zooms on the contact line.

  Figure 9.8: Steady state solutions for the three dierent ow congurations. From left to right: bead pulling (Re = 12, Ca = 1.33), beneath the liquid curtain (Re = 30, Ca = 2.6), heel formation (Re = 35, Ca = 1.36).

Figure 9

 9 Figure 9.9: Stability analysis of the reduced curtain coating system and comparison with computations of Liu [Liu et al. 2016].

  Error (%) inection point distance.

Figure 9 .

 9 Figure 9.10: Convergence study of the reduced curtain coating system for Re = 30, Ca = 2.6 and λ = 10 µm. The reference solution is taken for 64 grid points per slip length.

Figure 9 .

 9 Figure 9.12: Stability analysis and comparison with the experiments of Marston and computations of Liu. The model parameters are: µ l = 117 mPa.s, µ g = 0.018 mPa.s, σ = 67 mN.m -1 , h c = 2.6 cm, θ m = 67°and λ = 450 nm.

  IF he superEslip model with further regulriztion of the veloity y introduE ing seond order slip lengthF PF he generlized xvier oundry ondition where the unompensted oung9s stress is tken into ountF QF he toy model where we relte the mss di'usion existing in the phseE(eld to simple dynmi ngle reltionF 10.4 Summary of models and applications . . . . . . . . . . . 143 10.1 The super-slip model 10.1.1 The porous substrate analogy

Figure 10 . 1 :

 101 Figure 10.1: Dimensionless radius as a function of the dimensionless time in logarithmic scale for λ 1 = 0, 0.1 and λ 2 = 0, 0.01.

Figure 10

 10 Figure 10.2: Local capillary number depending on the contact line speed as a function of the dimensionless time in logarithmic scale for λ 1 = 0, 0.1 and λ 2 = 0, 0.01.

Figure 10 . 3 :

 103 Figure 10.3: Contact angle as a function of the dimensionless time in logarithmic scale for λ 1 = 0, 0.1 and λ 2 = 0, 0.01.

Figure 10

 10 Figure 10.4: Shear stress at the contact line as a function of the dimensionless time in logarithmic scale for λ 1 = 0, 0.1 and λ 2 = 0, 0.01.

Chapter 10 .

 10 Advances in models for contact lines 10.2 The generalized Navier boundary condition 10.2.1 The uncompensated Young stress he min ide of the generlized xvier oundry ondition @qxfgA is to dd the unompensted oung stress to the xfg in order to llow devition etween the dynmi ontt ngle nd the equilirium oneF his stress omes from the devition of the )uidE)uid interfe from the stti on(gurtion in et al. PHHQD in et al. PHHTF elthough the qxfg ws formulted in di'use interfe modelD prike prike PHPH showed tht the model is wellEposed in the shrp interfe limitF he oundry ondition for the tngentil veloity is

Figure 10 . 6 :

 106 Figure 10.6: Example of steady state meniscus for Ca = 0.03. The solid substrate on the left is moving up at a speed U . The color map shows the volume fractions with the red phase corresponding to the viscous liquid being withdrawn.

Figure 10

 10 Figure 10.7: Normalized contact line position as a function of the dimensionless time for Ca = 0.1 and θ e = 80 °, 90 °.

Figure 10

 10 Figure 10.8: Contact angle as a function of the dimensionless time in logarithmic scale for Ca = 0.1 and θ e = 80 °, 90 °. The blue lines correspond to the equilibrium angles.

Figure 10

 10 Figure 10.9: Shear stress at nal time as a function of x for Ca = 0.1 and θ e = 80 °, 90 °.

  (a) No-slip boundary condition with a constant contact angle θe. (b) Generalized Navier boundary condition.

Figure 10 .

 10 Figure 10.10: Comparison between the no-slip model and the GNBC in the withdrawing plate setup. The capillary number is Ca = 0.03 and the equilibrium angles θ e = 60 °, 90 °. The coarse, medium and ne grids correspond to grid spacings ∆/l c = 0.01953125, 0.009765625, 0.004882813 respectively.

Figure 10 .

 10 Figure 10.11: Schematic of the full PF and the toy model VOF computational domains with their respective boundary conditions.

Figure 10

 10 Figure 10.12: Dimensionless diameter as a function of time for the PF and VOF simulations. The inset shows the spurious initial oscillation of the displacement in the VOF simulation.

  f (θapp, κ, ∆, h) etion IHFPFI engle trnsport eppendix p ux = f (λ, U ) etion WFPFI θ = θe gurtin oting etion WFQFI pulln et al. PHPH ux = f (λ, U ) θ d = f (θe, Ca, L, l) wh v© is et al. PHPH ux = f (ε, λ, U )

  Figure10.15: Summary of the dierent models for the dynamic contact line.

  fFSA hereforeD the disretiztion of the visous trnsport termD ∇ • (2µs)D is performed similrly to tht of the strinErte opertorD y trnslting the de(nition of the pities to yield ∀α ∈ {x, y} , visc α (S) = β B αβ δS αβ δξ β , where S = (S αβ ) is de(ned s funtion of U nd D y iqutions fFQ nd fFRFConvective transport term he onvetive term in the momentum trnsport eqution long α ∈ {x, y} is rewritten in onservtive form using the divergeneEfree onditionD

  points in the x-direction ny = 64 # Number of points in the y-direction x = LinRange(0.0, R * L0, nx+1) # x and y segements with the origin at (0,0) y = LinRange(0.0, L0, ny+1) # Define the numerical and physical parameters num = Numerical(CFL = 0.5, # CFL condition, defines the time step x

  Figure D.1: Initialization of a 2D drop of diameter D = 1 in a 2 × 2 domain without contact (left gures) or with contact with the bottom boundary (right gures). From top to bottom, the shift in the x-direction is 0, ∆/4 and ∆/2. The color map corresponds to the volume fraction c.

  Figure E.2: Validation of the super-slip implementation on a Poiseuille ow with λ 1 = 0, ∆ and λ 2 = 0.

  Figure G.2: Validation of the sink velocity imposed at the computational boundary placed at +δ with δ = 0.02, 0.05, 0.2 denoted by the dotted black lines in the gures. The equilibrium angle is θ e = 70 °for the three cases considered.

  The level set method . . . . . . . . . . . . . . . . . . . . . 9 1.3 Thermal convection and uid ow extension . . . . . . . 12

	Contents
	1.1 Classical two-phase Stefan problem formulation . . . . . 7
	1.2

1.1 Classical two-phase Stefan problem formulation

  Minimization problemhis hpter presents the optimiztion proedure for the twoEphse tefn prolem desried in etion IFI together with the respetive ost funtionl nd the orreE sponding djoint equtionsF he ontrol vrile is xeumnn oundry ondition tht ts on the domin oundryF he djoint het equtions nd djoint level set equtions re derived in the sme fshion s fernuer 8 rerzog PHIIF he resultE ing grdientEsed optimiztion lgorithmD s well s the derivtiveEfree methodsD used throughout the vrious ses onsidered in this study re desriedF fy onsidE ering the twoEphse tefn prolem with the qisEhomson reltionD the generi

	Chapter 2
	Adjoint-based optimization
	algorithm

Contents 2.1 Minimization problem . . . . . . . . . . . . . . . . . . . . 15 2.2 Continuous adjoint derivation . . . . . . . . . . . . . . . . 18 2.2.1 Lagrange functional . . . . . . . . . . . . . . . . . . . . . 18 2.2.2 Adjoint temperature . . . . . . . . . . . . . . . . . . . . . 19 2.2.3 Adjoint level set . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.4 The gradient equation . . . . . . . . . . . . . . . . . . . . 25 2.3 Gradient-based optimization procedure . . . . . . . . . . 26 2.1 forwrd twoEphse tefn prolem @pA n e rest in the level set frmework given elowD whih is lter used in the djoint derivtionF Chapter 2. Adjoint-based optimization algorithm pind funtion T

  4.3. The advection equation 57 in wikul et al. PHIRF gonsider p to e (nite volume of ell nd e pq the edge etween p nd qD q ∈ N (p) where N (p) is the set of neighoring (nite volumesF he length of e pq is normlized to 1F vet n pq e the outer norml vetor to e pq with respet to pF pinllyD let us denote φp the onstnt reonstrution of ϕ p in the (nite volume p nd φpq the onstnt reonstrution of ϕ on e pq F sntegrting volume pD iqution RFW yields

	ˆp ∂ϕ ∂t	dx + ˆp ∇ • F ϕ	∇ϕ |∇ϕ|	dx -ˆp ϕ∇ • F	∇ϕ |∇ϕ|	dx = 0.	@RFIHA
	epplying tokes9 theorem nd using the onstnt reonstrutions of ϕD we otin
		ˆp ∂ϕ ∂t -	dx + q∈N (p) q∈N (p) φp ˆepq φpq F ϕ ˆepq 1 F |∇ϕ|	1 |∇ϕ| ∇ϕ • n pq ds = 0. ∇ϕ • n pq ds

vet us denote |∇ϕ pq | the reonstruted rmiltonin |∇ϕ| on the edge e pq nd (ϕ q -ϕ p ) /1 the norml derivtive ∇ϕ • n pq on the sme edgeF e otin the following expression ˆp ∂ϕ ∂t dx + q∈N (p)

  Actuation function parametrization . . . . . . . . . . . . 85 6.2 Adjoint-based optimization cases . . . . . . . . . . . . . . 86 6.2.1 Asymmetric melting of a solid circle . . . . . . . . . . . . 86 6.2.2 Controlling the Mullins-Sekerka instability . . . . . . . . . 88 6.2.3 Growing crystals with topology changes . . . . . . . . . . 90 6.3 Comparison with derivative-free methods . . . . . . . . . 93

	Chapter 6
	Shape optimization results
	Contents
	6.1

  Table 6.1: Optimization parameters and nal results of the dierent considered cases. The columns J and ∇J correspond to the number of FP and AP calls respectively.

	6.3. Comparison with derivative-free methods		93
	6.3 Comparison with derivative-free methods
	irle	1	10	0	10 3	66	16	1.81 × 10 3
	wullinsEekerk 1	10 0.1	10 4	54	10	1.89 × 10 3
	qrowing rystls 1	1	0	10 2	25	20	3.21 × 10 3

  iggers PHHRD ilson et al. PHHTD pulln et al. PHPHD v© is et al. PHPHD v© is et al. PHPP will e the fous our study nd will e detiled in ghpter VF sn the ontext of numeril methods for prtil di'erentil equtionsD the resoE lution of the the inompressile twoEphse xvierEtokes equtions remins hlE lengeF he inompressiility onstrint must e tking into ount t eh instnt nd the high density nd visosity rtios in presene of distorted interfes leds to illEonditioned liner systemsF woreoverD in order to void ny imlne etween disrete surfe tension nd pressure grdient termsD n urte representtion of the interfil fore is required rdovelli 8 leski IWWWD ryggvson et al. PHIID opinet PHIVF 96 Chapter 7. Two-phase Navier-Stokes equations and Phase-Field extensions he inompressileD vrileEdensityD xvierEtokes equtions with surfe tenE sion n e written s the momentum eqution ugmented with the pillry e'et

  PHHTD viu et al. PHITD viu et al. PHIVD pulln et al. PHPHD v© is et al. PHPHF sn etion WFPFID we will present our results on urtin oting setup using the xvier oundry onditionF xonethelessD the slip model still su'ers from logrithmi divergene of the pillry pressure s shown in hevuhelle et al. PHHUF e further regulriztion of the solution is possile y onsidering seondEorder slip @in nlogy with the porous medi theoryA Paradoxes and singularities in the dynamic case 103sn tht se the ontt line speed is lso proportionl to the urvture of the veloity pro(le in the vertil diretionF sn our workD this model will e referred to s the superEslip model @etion IHFIAF enother prmeter to tke into ount is the ontt ngleF hen the ontt line enters in motionD the ngle formed y the )uidE)uid interfe nd the solid sustrte vriesD eoming dynmi ontt ngleD denoted θ d F he modeling of this ngle is very ontroversil mong expertsD there is no onsensul theory on whether it is n preditle quntity ! funtion of the ontt line speed nd the ulk prmeters ! or n output of the system gox IWVTD hikhmurzev IWWRD flke et al. IWWWD flke 8 hikhmurzev PHHPD iggers PHHRD iggers 8 ivns PHHSF sn the goxEoinov theoryD the dynmi ontt ngle is funtion of the normlized ontt line Ca CL g (θ d ) = g (θ e ) + Ca CL log PHIW the uthors demonstrted tht there exists seond prdox relted to tht dynmi ontt ngleF prom kinemtis reltionsD they showed tht the sher veloity t the ontt line is proportionl to the vrition of the dynmi ontt ngle θ d CL ̸ = 0F he wy to relx this singulrity is to llow devition etween the dynmi ontt ngle nd the equilirium one @θ e AF sn the generlized xvier oundry ondition @qxfgAD derived from thermodynmis onsidertions in et al. PHHQD in et al. PHHTD 104 Chapter 8. Singularities in ows with contact lines the unompensted oung9s stress is introdued s n dditionl term in the xfg

	8.2. L l	,	@VFSA
		∂u ∂y	| y=0 ∼	dθ d dt	.	@VFTA
	gonsidering now stedy stte motion where the ontt line speed is not null with
	slip model @iqution VFPAD we hve
		U CL ∼ λ	dθ d dt	.	@VFUA
	sn stedy stteD the ngle is onstntD implying tht
				dθ d dt	= 0,	@VFVA
	whih ontrdits previous ssumption on the ontt line speed U
	U CL ∼ λ	∂u ∂y	| y=0 + λ 2 ∂ 2 u ∂y 2 | y=0 .	@VFRA

where g is funtion desriing the urved interfe @de(ned lterAF he other pE rmeters re L nd l orresponding to two distint length slesF e omprison etween di'erent options for L nd l n e found in vegendre 8 wglio PHISF hey onluded tht the models inorporting dynmi ontt ngles with nuE meril uto' length l etter represent some experimentsF sn etion WFQFPD we will use dynmi ngle model implemented in our yp frmework nd ompre it to spreding drop experimentF woreoverD in prike et al.

  The Volume-Of-Fluid method . . . . . . . . . . . . . . . . 105 9.2 Classical models for contact lines . . . . . . . . . . . . . . 109 9.2.1 The Navier boundary condition . . . . . . . . . . . . . . . 109 9.2.2 Cox-Voinov relation for dynamic contact angles . . . . . . 110 9.3 Application to physical systems . . . . . . . . . . . . . . . 111 PHIIF roweverD this methodD suh s othersD is sujet to spurious urrentsD with or without ontt lines @eppendix hAF hese urrents re lower when onsidering high pillry numersD whih is usully the se in the ontext of moving ontt linesF sn our study on the moving ontt lineD we use the free softE wre Basilisk D pltform for the solution of prtil di'erentil equtions on dptive grtesin meshesD developed y téphne opinet nd oE uthors efkhmi 8 fussmnn PHHVD efkhmi 8 fussmnn PHHWD opinet PHHWD opinet PHISD efkhmi et al. PHIUD opinet PHIVF es mentioned in ghpter UD we onsider the inompressile xvierEtokes equtions with vrile density nd surfe tensionX

	Chapter 9
	Numerical methods for contact
	lines
	Contents
	9.1

9.

3.1 Onset of wetting failure in a curtain coating setup . . . . 111 9.3.2 Squalane drop spreading with dynamic contact angle . . . 119 9.1 The Volume-Of-Fluid method he yp method for representing )uid interfes oupled with )ow solver is wellEknown to e suited for solving interfil )ows rdovelli 8 leski IWWWD opinet 8 leski IWWWD ryggvson et al.

  dynmi ontt line introdues prdox t the triple pointD where the noEslip oundry ondition or hirihlet oundry ondition t the solid interfe indues nonEintegrility of the solutionF xeverthelessD in the yp method suh ontrdiE tion does not tke ple s the volume frtion is dveted using the veloity hlfEE ell wy from the wllF sn their pperD the uthors of efkhmi et al. PHIUD showed tht there exists numeril slip9 tht is mesh dependentF sn order to ontrol suh slipD it is useful to introdue xvier oundry ondition @xfgA to expliitly de(ne physil slip length λD tht will e used s (tting prmeter in our simulE tions vegendre PHIQD ui et al. PHIRD vegendre 8 wglio PHISD pulln et al. PHPHD v© is et al. PHPHD v© is et al.

	9.2. Classical models for contact lines	109
	the rightEhnd neighor is delimited with dshed lineF he frtion of this volume
	oupied y the (rst phse is indited y the drk grey tringleF
	9.2 Classical models for contact lines
	9.2.1 The Navier boundary condition

he

  is the dynmi numeril ngle imposed t the ontt lineD θ e is stti equilirium ontt ngleD property of the sustrte nd the pillry numer is expressed s funtion of the ontt line speed U CL

	9.3. Application to physical systems			111
	rtio L/l etween the two length sles where L numeril ftor to e determinedF
	prom generl nlysisD gox otined
	g (θ d ) = g (θ e ) + Ca CL log	L l	,	@WFIHA
	where θ d
	9.2.2 Cox-Voinov relation for dynamic contact angles
	he theory of the moving ontt line t smll pillry numers ws founded y
	oinov oinov IWUT nd generlized to ritrry visosity rtios y gox gox IWVTF
	he ltter lri(ed the struture of low pillry numer prolems in terms of the

  Cn ≪ 1D the solution to this eqution hs oundry lyer hrterD nd we introdue new rdil oordinte p suh tht 10.3. The toy model: mimic Phase-Field at lower cost 137 r = R + Cn •p nd expnd iqution IHFIH t zero nd (rst order

	1 Re	1 δ 2 ∼	1 Ca Cn Re	ϕ δ	,
	nd	1 δ	∼		1 P e	ϕ δ 2 .
	olving for δ yields
		δ =	Ca Cn Pe	,	@IHFWA
	whih grees with the sling of ueEhouEpeng ue et al. PHHUF e now derive
	expression for hemil potentil on urved interfeF
	he (rst step is to lulte the vlue of the hemil potentil on urved
	interfeF por dropletD C pprohes +1 inside the dropletD nd -1 outsideF
	essuming plnr geometry nd irulr interfe of rdius RD the hemil
	potentil s de(ned in iqution IHFV n e written s
	ϕ = -Cn 2 ∂ 2 C ∂r 2 +	1 r	∂C ∂r	+ Ψ ′ (C) .	@IHFIHA
	ine the dimensionless interfe width is

  Bibliography euEelEud et al. PHIV wF yF euEelEudD F opinet nd rF eF helepiF A conservative and well-balanced surface tension modelF tournl of gomputE tionl hysisD volF QUID pges VWT!WIQD ytoer PHIVF @gited on pge ITQFA efkhmi 8 fussmnn PHHV F efkhmi nd wF fussmnnF Height functions for applying contact angles to 2D VOF simulationsF sntF tF xumerF wethF plFD volF SUD noF RD pges RSQ!RUPD PHHVF @gited on pges IHSD IHT nd IIIFA efkhmi 8 fussmnn PHHW F efkhmi nd wF fussmnnF Height functions for applying contact angles to 3D VOF simulationsF sntF tF xumerF wethF plFD volF TID noF VD pges VPU!VRUD PHHWF @gited on pges IHSD IHT nd IIIFA efkhmi et al. PHIU F efkhmiD tF fuongiornoD eF quionD F opinetD F rE dovelli nd F leskiF Transition in a numerical model of contact line dynamics and forced dewettingF tF gomputF hysFD volF QURD PHIUF @gited on pges IHSD IHWD IPUD IQH nd IRSFA elmgren IWWQ F elmgrenF Variational algorithms and pattern formation in dendritic solidicationF t gompF hysFD volF IHTD noF PD pges QQU!QSRD IWWQF @gited on pge UHFA et al. PHHS wF fenziD qF rF qolu nd tF viesenF Numerical solution of saddle point problemsF et xumeriD volF IRD S PHHSF @gited on pge ISUFA fernuer 8 rerzog PHII wF uF fernuer nd F rerzogF Optimal Control of the Classical Two-Phase Stefan Problem in Level Set FormulationF sew tF iF gomputFD volF QQD noF ID pges QRP!QTQD PHIIF @gited on pges QD IHD ISD PP nd ISIFA feyer 8 veveque PHHT F F feyer nd F tF vevequeF Analysis of a One-Dimensional Model for the Immersed Boundary MethodF volF PWD pges QQP! QTRD U PHHTF @gited on pge QHFA feznson et al. PHIU tF feznsonD eF idelmnD F urpinski nd F fF hhF Julia: A Fresh Approach to Numerical ComputingF sew eviewD volF SWD noF ID pges TS!WVD PHIUF @gited on pge TUFA flke 8 hikhmurzev PHHP F hF flke nd F hF hikhmurzevF Dynamic Wetting by Liquids of Dierent ViscosityF tournl of olloid nd interfe sieneD volF PSQD pges IWT!PHPD PHHPF @gited on pges QD IHQD IIT nd IIUFA flke et al. IWWW F hF flkeD wF frke nd F hF hikhmurzevF Experimental evidence of nonlocal hydrodynamic inuence on the dynamic contact angleF hysis of pluidsD volF IID noF VD pges IWWS!PHHUD IWWWF @gited on pges RD WSD IHQD IIQD IIT nd IIUFA fonn PHHW hF fonnF Wetting and spreadingF eviews of wodern hysisD volF VID noF PD pges UQW!VHSD PHHWF @gited on pge IQSFA frmn et al. PHIS uF frmnD F eF yliver nd F mnF Adjoint-based sensitivity analysis of amesF gomF heoF wodellingD volF IWD noF ID pges PW!STD PHISF @gited on pge PFA frttkus 8 weiron IWWP uF frttkus nd hF sF weironF Numerical simulations of unsteady crystal growthF sew tF epplF wthFD volF SPD pge IQHQD IWWPF @gited on pge IFA frown et al. PHHI hF vF frownD F gortez nd wF vF winionF Accurate Projection Methods for the Incompressible NavierStokes EquationsF tournl of gomputtionl hysisD volF ITVD R PHHIF @gited on pge ISUFA ghn 8 rillird IWSV tF F ghn nd tF iF rillirdF Free Energy of a Nonuniform System. I. Interfacial Free EnergyF he tournl of ghemil hysisD volF PVD noF PD pges PSV!PTUD IWSVF @gited on pge WVFA glhoun 8 veveque PHHH hF glhoun nd F tF vevequeF A Cartesian Grid Finite-Volume Method for the Advection-Diusion Equation in Irregular Geome-triesF tournl of gomputtionl hysisD volF ISUD I PHHHF @gited on pge QHFA grlson PHIP eF grlsonF Capillarity and dynamic wettingF hh thesisD ur oyl snstitute of ehnologyD PHIPF @gited on pges R nd WUFA ghen et al. IWWU F ghenD fF werrimnD F ysher nd F merekF A Simple Level Set Method for Solving Stefan ProblemsF t gompF hysFD volF IQSD noF ID pges V!PWD IWWUF @gited on pges VD IHD USD UT nd VVFA gheny 8 fotell PHIH F gheny nd yF fotellF The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous ows in complex moving geometries with good conservation proper-tiesF tournl of gomputtionl hysisD volF PPWD pges IHRQ!IHUTD PHIHF @gited on pges QH nd ISRFA ghorin IWTV eF tF ghorinF Numerical Solution of the Navier-Stokes EquationsF emerin wthemtil oietyD volF PPD IWTVF @gited on pge ISUFA gox IWVT F qF goxF The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous owF tF pluid wehFD volF ITVD pges ITW!IWRD IWVTF @gited on pges RD IHQ nd IIHFA grnk IWVU tF grnkF pree nd woving foundry rolemsF yxford siene puE litionsF glrendon ressD IWVUF @gited on pge TVFA hevuhelle et al. PHHU yF hevuhelleD gF tossernd nd F leskiF Forced dewetting on porous mediaF tournl of pluid wehnisD volF SURD pge QRQ!QTRD PHHUF @gited on pges RD IHP nd IPIFA hupont 8 vegendre PHIH tFEfF hupont nd hF vegendreF Numerical simulation of static and sliding drop with contact angle hysteresisF tournl of gomputE tionl hysisD volF PPWD noF UD pges PRSQ!PRUVD PHIHF @gited on pge IIWFA hurismy 8 elonso PHIP uF hurismy nd tF elonsoF Adjoint-based techniques for uncertainty quantication in turbulent ows with combustionF RPnd esee pluid hynmis gonferene nd ixhiitD pges PS!PVD PHIPF @gited on pge PFA iggers 8 ivns PHHS tF iggers nd F ivnsF Comment on "Dynamic Wetting by liquids of dierent viscosity", by T. D. Blake and Y. D. ShikhmurzaevF tournl of olloid nd interfe sieneD volF PVHD pges SQU!VY disussion SQWD PHHSF @gited on pges RD IHQ nd IIUFA iggers PHHR tF iggersF Toward a description of contact line motion at higher capillary numbersF hysis of pluidsD volF ITD noF WD pges QRWI!QRWRD PHHRF @gited on pges RD WS nd IHQFA pvier et al. PHIW fF pvierD tF urseed nd vF huheminF RayleighBénard convection with a melting boundaryF tF pluid wehFD volF VSVD pges RQU!RUQD tnury PHIWF @gited on pges QD IP nd VIFA Bibliography pikl et al. PHPH eF piklD F ve ghende nd F ydiF Control and Optimization of Interfacial Flows Using Adjoint-Based TechniquesF pluidsD volF SD noF QD PHPHF @gited on pges P nd IHFA poures et al. PHIR hF pouresD gF gul(eld nd F hmidF Optimal mixing in twodimensional plane Poiseuille ow at nite Péclet numberF t pluid wehFD volF URVD pges PRI!PUUD PHIRF @gited on pge PFA prnk IWSH pF gF prnkF Radially symmetric phase growth controlled by diusionF roF F oF vondD IWSHF @gited on pge TWFA prike et al. PHIV wF prikeD wF uöhne nd hF fotheF On the Kinematics of Contact Line MotionF roF epplF wthF wehFD volF IVD noF ID PHIVF @gited on pge IPTFA prike et al. PHIW wF prikeD wF uöhne nd hF fotheF A kinematic evolution equation for the dynamic contact angle and some consequencesF hysi hX xonlF henFD volF QWRD pges PT!RQD PHIWF @gited on pges IHQ nd IPTFA prike PHPH wF prikeF Mathematical modeling and Volume-of-Fluid based simulation of dynamic wettingF hh thesisD PHPHF @gited on pges RD IPTD IUQ nd IURFA pujiok IWUV F pujiokF hh thesisD grnegieEwefon niversityD IWUVF @gited on pge WFA pulln et al. PHPH F pullnD F leski nd F opinetF Dynamic wetting failure in curtain coating by the Volume-of-Fluid method -Volume-of-Fluid simulations on quadtree meshesF iurF hysF tF peil opisD volF PPWD noF IHD pges IWPQ!IWQRD PHPHF @gited on pges RD WSD IHPD IHWD III nd IRSFA pulln et al. PHPQ F pullnD F ve ghende nd F ydiF Adjoint-based optimization of two-dimensional Stefan problemsF tournl of gomputtionl hysisD volF RUSD pge IIIVUSD PHPQF @gited on pge QFA pulln PHIW F pullnF tudy of the dynmi ontt line y volume of )uid method on hierrhil dptive meshF wster9s thesisD oronne niversitéD PHIWF @gited on pges IRR nd IRSFA qiles 8 iere PHHH wF fF qiles nd xF eF iereF An Introduction to the Adjoint Approach to DesignF plow urF gomusFD volF TSD noF QD pges QWQ!RISD PHHHF @gited on pge PFA qoldstein et al. IWWQ hF qoldsteinD F rndler nd vF irovihF Modeling a No-Slip Flow Boundary with an External Force FieldF tournl of gomputtionl hysisD volF IHSD R IWWQF @gited on pge QHFA rrlow 8 elh IWTS pF rF rrlow nd tF iF elhF Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free SurfaceF hysis of pluidsD volF VD noF IPD pges PIVP!PIVWD IWTSF @gited on pge QQFA rssn et al. PHPI eF rssnD F ydiD F ve ghendeD rF itsh nd eF ettiliF Adjoint-Based Sensitivity Analysis of Steady Char Burn outF gomF heoF wodelingD volF PSD noF ID pges WT!IPHD PHPIF @gited on pges I nd QFA reht PHIP pF rehtF New development in FreeFem++F tF xumerF wthFD volF PHD noF QERD pges PSI!PTSD PHIPF @gited on pge IRIFA rinze 8 iegenlg PHHU wF rinze nd F iegenlgF Optimal control of the free boundary in a two-phase Stefan problemF t gomputF hysFD volF PPQD noF PD pges TSU!TVRD PHHUF @gited on pge QFA ro'mnn 8 prekels IWVP uF rF ro'mnn nd tF prekelsF Real-time control of the free boundary in a two-phase Stefan problemF xumeril puntionl enlysis nd yptimiztionD volF SD noF ID pges RU!UTD IWVPF @gited on pge QFA ruh 8 riven IWUI gF ruh nd vF iF rivenF Hydrodynamic model of steady movement of a solid/liquid/uid contact lineF tF golloid nd snterfe iFD volF QSD pges VS!IHID IWUIF @gited on pges QD IHP nd IQSFA tqmin IWWW hF tqminF Calculation of two-phase NavierStokes ows using phase-eld modelingF tF gomputF hysFD volF ISSD noF ID pges WT!IPUD IWWWF @gited on pge WVFA tqmin PHHH hF tqminF Contact-line dynamics of a diuse uid interfaceF tF pluid wehFD volF RHPD pges SU!VVD PHHHF @gited on pges RD WU nd WVFA tmeson et al. IWWV eF tmesonD vF wrtinelli nd xF eF iereF Optimum Aerodynamic Design Using the Navier-Stokes EquationsF heorF gompF pluid hynFD volF IHD noF ID pges PIQ!PQUD IWWVF @gited on pge PFA tmeson IWVV eF tmesonF Aerodynamic design via control theoryF t ienti( gomputingD volF QD noF QD pges PQQ!PTHD IWVVF @gited on pge PFA tvierre et al. PHHT iF tvierreD gF uikD pF tF ermolen nd F vn der wgF A comparison of numerical models for one-dimensional Stefan problemsF t gompF epplied wthemtisD volF IWPD noF PD pges RRS!RSWD PHHTF @gited on pge IFA tohnsen 8 golell IWWV rF tohnsen nd F golellF A Cartesian Grid Embedded Boundary Method for Poisson's Equation on Irregular DomainsF t gompF hysFD volF IRUD noF ID pges TH!VSD IWWVF @gited on pge SQFA tuniper PHIH wF tuniperF Triggering in the horizontal Rijke tube: non-normality, transient growth and bypass transitionF tF pluid wehFD volF TTUD pges PUP! QHVD PHIHF @gited on pge PFA turi 8 ryggvson IWWT hF turi nd qF ryggvsonF A Front-Tracking Method for Dendritic SolidicationF t gompF hysFD volF IPQD noF ID pges IPU!IRVD IWWTF @gited on pges ID W nd UTFA ung 8 rs IWWS F ung nd xF rsF Control of the freezing interface motion in two-dimensional solidication processes using the adjoint methodF snterntionl t xumeril wethods in ingineeringD volF QVD pges TQ!VHD IWWSF @gited on pge QFA unner IWVS F unnerF Control of Stefan problems by means of linear-quadratic defect minimizationF xumerishe wthemtikD volF RTD noF QD pges RPW!RRPD IWVSF @gited on pge QFA v© is et al. PHPH F v© isD F tohnssonD F pullnD fF ressD qF emergD F fgheri nd F leskiF Steady moving contact line of water over a noslip substrate: Challenges in benchmarking phase-eld and volume-of-uid methods against molecular dynamics simulationsF iurF hysF tF peF opFD volF PPWD noF IHD pges IVWU!IWPID eptemer PHPHF @gited on pges RD WSD IHPD IHWD IIH nd IRSFA v© is et al. PHPP F v© isD wF ellegrinoD tF undinD qF emergD F leskiD fF ress nd F fgheriF Nanoscale sheared droplet: volume-of-uid, phase-eld and no-slip molecular dynamicsF tournl of pluid wehnisD volF WRHD pge eIHD PHPPF @gited on pges RD WS nd IHWFA vi 8 eskin PHHH wF gF vi nd gF F eskinF An Immersed Boundary Method with Formal Second-Order Accuracy and Reduced Numerical ViscosityF tourE nl of gomputtionl hysisD volF ITHD S PHHHF @gited on pge QHFA vnger IWVH tF F vngerF Instabilities and pattern formation in crystal growthF eviews of wodern hysisD IWVHF @gited on pges ID V nd UPFA vvi 8 wrmur PHHR fF vvi nd eF wrmurF The exponential power law: partial wetting kinetics and dynamic contact anglesF golloids nd urfes eX hysiohemil nd ingineering espetsD volF PSHD noF ID pges RHW!RIRD PHHRF @gited on pges R nd IIWFA vegendre 8 wglio PHIS hF vegendre nd wF wglioF Comparison between numerical models for the simulation of moving contact linesF gomputers 8 pluidsD volF IIQD pges P!IQD wy PHISF @gited on pges RD IHQD IHWD IIW nd IRSFA vegendre PHIQ hF vegendreF Numerical simulation of spreading dropsF golloids urfD volF RQPD pge PW!QUD HW PHIQF @gited on pge IHWFA vemke et al. PHIQ wF vemkeD tF eiss nd tF esterhennF Adjoint-based analysis of thermoacoustic couplingF sgxeewD pges PITQ!PITTD PHIQF @gited on pge PFA vemke et al. PHIW wF vemkeD vF giD tF eissD rF itsh nd tF esterhehnF Adjoint-based sensitivity analysis of quantities of interest of complex combustion modelsF gomF heoF wodellingD volF PQD noF ID pges IVH!IWTD PHIWF @gited on pge PFA vimre et al. PHPQ eF vimreD F opinetD gF tosserndD F ue nd eF qhigoF A hybrid level-set / embedded boundary method applied to solidication-melt problemsF tournl of gomputtionl hysisD volF RURD pge IIIVPWD PHPQF @gited on pges ID QD IHD TPD TVD UV nd VIFA viuF 8 xoedl IWVW hF g viuF nd tF xoedlF On the Limited Memory Method for Large Scale OptimizationF wthemtil rogrmmingD volF QD noF RSD pge SHQ!SPVD IWVWF @gited on pge PTFA viu et al. PHIT gFEF viuD iF ndreD wF F grvlho nd F uumrF Dynamic wetting failure and hydrodynamic assist in curtain coatingF tournl of pluid wehnisD volF VHVD pge PWH!QISD PHITF @gited on pges IHPD IIID IIQD IIS nd IITFA viu et al. PHIV gFEF viuD wF F grvlho nd F uumrF Dynamic wetting failure in curtain coating: Comparison of model predictions and experimental obser-vationsF ghemil ingineering ieneD volF IWSD PHIVF @gited on pges IHP nd IIUFA wrsden et al. PHHV eF vF wrsdenD tF eF peinstein nd gF eF ylorF A computational framework for derivative-free optimization of cardiovascular geome-triesF gomputer wethods in epplied wehnis nd ingineeringD volF IWUD noF PIEPRD pges IVWH!IWHSD PHHVF @gited on pge PFA wrston et al. PHHW tF yF wrstonD F rwkinsD F F heent nd wF tF rF imE monsF Inuence of surfactant upon air entrainment hysteresis in curtain coatingF ixperiments in pluidsD volF RTD noF QD pges SRW!SSVD wr PHHWF @gited on pges QD IIQD IIT nd IIUFA wvriplis IWWS hF tF wvriplisF Unstructured mesh generation and adaptivity. ehnil reportD IWWSF @gited on pge QHFA wikul 8 yhlerger PHIH uF wikul nd wF yhlergerF A New Level Set Method for Motion in Normal Direction Based on a Semi-Implicit Forward-Backward Diusion ApproachF sew tF iF gomputFD volF QPD noF QD pges ISPU!ISRRD PHIHF @gited on pges ST nd SVFA 188 Bibliography wikul et al. PHIR uF wikulD wF yhlerger nd tF ránF Inowimplicit/outow-explicit nite volume methods for solving advection equa-tionsF epplied xumeril wthemtisD volF VSD pges IT!QUD PHIRF @gited on pges ST nd SUFA win PHIH gF winF On reinitializing level set functionsF t gompF hysFD volF PPWD noF VD pges PUTR!PUUPD epril PHIHF @gited on pges TP nd TQFA wogensen 8 iseth PHIV F uF wogensen nd eF xF isethF Optim: A mathematical optimization package for JuliaF tournl of ypen oure oftwreD volF QD noF PRD pge TISD PHIVF @gited on pge TVFA worinishi et al. IWWV F worinishiD F F vundD yF F silyev nd F woinF Fully Conservative Higher Order Finite Dierence Schemes for Incompressible FlowF tournl of gomputtionl hysisD volF IRQD noF ID pges WH!IPRD IWWVF @gited on pges QH nd QIFA worinishi PHIH F worinishiF Skew-symmetric form of convective terms and fully conservative nite dierence schemes for variable density low-Mach number owsF tournl of gomputtionl hysisD volF PPWD noF PD pges PUT!QHHD PHIHF @gited on pges QHD ISS nd ISTFA wullins 8 ekerk IWTR F F wullins nd F pF ekerkF Stability of a Planar Interface During Solidication of a Dilute Binary AlloyF tournl of epplied hysisD volF QSD noF PD pges RRR!RSID IWTRF ulisherX emerin snstitute of hysisF @gited on pges ID UH nd VVFA xoedl 8 right PHHT tF xoedl nd F tF rightF xumeril optimiztionF pringer series in opertion reserh nd (nnil engineeringF pringerD Pnd éditionD PHHTF @gited on pge ISUFA ysher 8 ethin IWVV F ysher nd tF eF ethinF Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formula-tionsF t gomputF hysFD volF UWD noF ID pges IP!RWD IWVVF @gited on pge IFA yu 8 tmeson PHII uF yu nd eF tmesonF Unsteady Adjoint Method for the Optimal Control of Advection and Burger's Equations Using High-Order Spectral Dierence MethodF emerin snstitute of eeronutis nd estronutisD pges I!IVD PHIIF @gited on pge PFA ndey et al. PHIV eF ndeyD tF hF heel nd tF humherF Turbulent superstructures in Rayleigh-Bénard convectionF xt gommunD volF WD noF ID pge PIIVD wy PHIVF xumerX I ulisherX xture ulishing qroupF @gited on pge IPFA eng et al. IWWW hF engD fF werrimnD F ysherD rF ho nd wF ungF A PDE-Based Fast Local Level Set MethodF t gompF hysFD volF ISSD noF PD pges RIH!RQVD IWWWF @gited on pge SRFA erot PHII tF fF erotF Discrete Conservation Properties of Unstructured Mesh SchemesF ennul eview of pluid wehnisD volF RQD pges PWW!QIVD jn PHIIF @gited on pge ISRFA eskin IWUP gF F eskinF Flow patterns around heart valves: A numerical methodF tournl of gomputtionl hysisD volF IHD IH IWUPF @gited on pge QHFA eskin IWUU gF F eskinF Numerical analysis of blood ow in the heartF tournl of gomputtionl hysisD II IWUUF @gited on pge QHFA eskin PHHQ gF F eskinF The immersed boundary methodF et xumeriD volF IID U PHHQF @gited on pge QHFA ierret et al. PHHU F ierretD F pF goelho nd rF utoF Multidisciplinary and multiple operating points shape optimization of three-dimensional compressor bladesF truturl nd wultidisiplinry yptimiztionD volF QQD noF ID pges TI!UHD PHHUF @gited on pge PFA ironneu IWUR yF ironneuF On optimum design in uid mechanicsF t pluid wehFD volF TRD noF ID pges WU!IIHD IWURF @gited on pge PFA opinet 8 leski IWWW F opinet nd F leskiF A front-tracking algorithm for accurate representation of surface tensionF sntF tF xumerF wethods pluidsD volF QHD noF TD pges UUS!UWQD IWWWF @gited on pges IHS nd IHUFA opinet PHHQ F opinetF Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometriesF tournl of gomputtionl hysisD volF IWHD noF PD pges SUP!THHD eptemer PHHQF @gited on pge IHVFA opinet PHHW F opinetF An Accurate Adaptive Solver for Surface-tension-driven Interfacial FlowsF tF gomputF hysFD volF PPVD noF ITD pges SVQV!SVTTD PHHWF @gited on pges TRD IHS nd IHUFA opinet PHIS F opinetF A quadtree-adaptive multigrid solver for the Serre GreenNaghdi equationsF tF gomputF hysFD volF QHPD pges QQT!QSVD PHISF @gited on pges IHS nd IHUFA opinet PHIV F opinetF Numerical Models of Surface TensionF ennuF evF pluid wehFD volF SHD noF ID pges RW!USD PHIVF @gited on pges WSD IHSD IHU nd ITQFA in et al. PHHQ FEF inD FEF ng nd F hengF Generalized Navier Boundary Condition for the Moving Contact LineF gommunF wthF iFD volF ID noF PD pges QQQ!QRID HT PHHQF @gited on pge IIUFA in et al. PHHQ FEF inD FEF ng nd F hengF Molecular scale contact line hydrodynamics of immiscible owsF hysF evF iD volF TVD noF ID pge HITQHTD PHHQF @gited on pges RD WVD IHQ nd IPTFA 190 Bibliography in et al. PHHT FEF inD FEF ng nd F hengF A variational approach to moving contact line hydrodynamicsF tF pluid wehFD volF STRD pges QQQ!PVD PHHTF @gited on pges RD IHQ nd IPTFA uirós odríguez et al. PHPP eF uirós odríguezD F pullnD F ve ghende nd F ydiF A Conservative Cartesian Cut Cell Method for the Solution of the Incompressible Navier-Stokes Equations on Staggered MeshesF in pressD PHPPF @gited on pges Q nd ISQFA in et al. PHIR F inD gF gul(eld nd F uerswellF Designing a more nonlinearly stable laminar ow via boundary manipulationF t pluid wehFD volF UQVD pges IXI!IPD PHIRF @gited on pge PFA ose IWWQ wF iF oseF An enthalpy scheme for Stefan problems in several dimen-sionsF eppF xumeril wthFD volF IPD pge PPWD IWWQF @gited on pge IFA rdovelli 8 leski IWWW F rdovelli nd F leskiF Direct numerical simulation of free-surface and interfacial owF ennuF evF pluid wehFD volF QID noF ID pges STU!THQD IWWWF @gited on pges WSD IHS nd IHUFA hmidt et al. PHIQ F hmidtD gF sliD F hulz nd xF F qugerF Threedimensional large-scale aerodynamic shape optimization based on shape cal-culusF esee tD volF SID noF IID pges PTIS!PTPUD PHIQF @gited on pge PFA egl et al. IWWV qF eglD gF uik nd pF ermolenF A conserving discretization for the free boundary in a two-dimensional Stefan problemF t gomputF hysFD volF IRID pges I!PID IWWVF @gited on pge IFA herer 8 ott PHHS F iF herer nd tF xF ottF High-order compact nitedierence methods on general overset gridsF tournl of gomputtionl hysisD volF PIHD pges RSW!RWTD IP PHHSF @gited on pge QHFA hikhmurzev IWWR F hF hikhmurzevF Mathematical modeling of wetting hy-drodynamicsF pluid hynF esFD volF IQD noF ID pges RS!TRD jn IWWRF @gited on pges RD WSD IHQ nd IIUFA noeijer 8 endreotti PHIQ tF rF noeijer nd fF endreottiF Moving Contact Lines: Scales, Regimes, and Dynamical TransitionsF ennul eview of pluid weE hnisD volF RSD noF ID pges PTW!PWPD PHIQF @gited on pge IQSFA tefn IVWI tF tefnF Über die Theorie der Eisbildung, insbesondere über die Eisbil-dung im PolarmeereF ennF hysik ghemieD volF RPD pges PTW!PVTD IVWIF @gited on pge IFA ui et al. PHIR F uiD rF hing nd F hF wF peltF Numerical Simulations of Flows with Moving Contact LinesF ennuF evF pluid wehFD volF RTD noF ID pges WU!IIWD PHIRF @gited on pge IHWFA n 8 rs PHHT vF n nd xF rsF A level set simulation of dendritic solidication with combined features of front-tracking and xed-domain meth-odsF tournl of gomputtionl hysisD volF PIID noF ID pges QT!TQD PHHTF @gited on pge UVFA ryggvson et al. PHII qF ryggvsonD F rdovelli nd F leskiF hiret nuE meril simultions of gs!liquid multiphse )owsF gmridge niversity ressD PHIIF @gited on pges WS nd IHSFA oinov IWUT yF oinovF Hydrodynamics of wettingF pluid hynFD volF IID pge UIRD IWUTF @gited on pges R nd IIHFA ilson et al. PHHT wF ilsonD tF ummersD F hF hikhmurzevD eF glrke nd F hF flkeF Nonlocal hydrodynamic inuence on the dynamic contact angle: Slip models versus experimentF hysil reviewF iD ttistilD nonlinerD nd soft mtter physisD volF UQD pge RITHTD PHHTF @gited on pges QD WSD IHP nd IIUFA oods IWWP eF F oodsF Melting and dissolvingF tF pluid wehFD volF PQWD pge RPWD IWWPF @gited on pge IFA right 8 xoedl PHHT F tF right nd tF xoedlF xumeril optimiztionF pringerD PHHTF @gited on pge PTFA ng IWWU F ngF The Adjoint Method for the Inverse Design of Solidication Processes with ConvectionF hh thesisD gornell niversityD IWWUF @gited on pge QFA ue 8 peng PHII F ue nd tF tF pengF Wall energy relaxation in the CahnHilliard model for moving contact linesF hysis of pluidsD volF PQD noF ID PHIIF @gited on pges IQS nd IRIFA ue et al. PHHU F ueD gF hou nd tF tF pengF Spontaneous shrinkage of drops and mass conservation in phase-eld simulationsF tournl of gomputtionl hysisD volF PPQD noF ID pges I!WD epril PHHUF @gited on pge IQTFA hn et al. PHHW FErF hnD tF hngD F vi nd rF FErF ghungF Adaptive Particle Swarm OptimizationF siii rnstions on ystemsD wnD nd gyerE netisD rt f @gyernetisAD volF QWD noF TD pges IQTP!IQVID PHHWF @gited on pge WQFA rler IWWS fF rlerF Stefan's work on solid-liquid phase changesF ingineering enlysis with foundry ilementsD volF ITD noF PD pges VQ!WPD IWWSF @gited on pge IFA

	182	Bibliography
	emerg PHHQ qF emergF Semisharp Phase Field Method for Quantitative Phase
	Change SimulationsF hysF evF vettFD volF WID he PHHQF @gited on pges R
	nd WUFA
	erkw 8 vm IWUU eF erkw nd F F vmF Computational Design of
	the Basic Dynamical Processes of the UCLA General Circulation ModelF sn
	wethods in gomputtionl hysisX edvnes in eserh nd epplitionsD
	volume IU of Methods in Computational Physics: Advances in Research and
	ApplicationsD pges IUQ!PTSF ilsevierD IWUUF @gited on pge QQFA
	ernold IWTT F ernoldF Sur la géométrie diérentielle des groupes de Lie de
	dimension innie et ses applications à l'hydrodynamique des uides par-
	faitsF ennles de l9institut pourierD volF ITD pges QIW!QTID IWTTF @gited
	on pge ISRFA
	fly et al. IWWU F flyD F hF qroppD vF gF wsnnes nd fF pF mithF Ecient
	Management of Parallelism in Object Oriented Numerical Software LibrariesF
	sn iF ergeD eF wF fruset nd rF F vngtngenD editorsD wodern oftwre
	ools in ienti( gomputingD pges ITQ!PHPF firkhäuser ressD IWWUF @gited
	on pge TVFA
	fell et al. IWVW tF fF fellD F golell nd rF wF qlzF A second-order projection
	method for the incompressible navier-stokes equationsF tournl of gomputE
	tionl hysisD volF VSD noF PD pges PSU!PVQD IWVWF @gited on pge IHVFA

fenzi

The two-phase Stefan problem 7

Chapter 5. Numerical simulations of two-phase Stefan problems

(a) ∆/4 shift without contact. (b) ∆/4 shift with contact.Figure D.3: Evolution of the maximal dimensionless velocity as a function of time for the Laplace numbers and grid resolutions indicated for a static drop with a ∆/4 initial shift, and with or without contact.

Remerciements

Chapter 9. Numerical methods for contact lines

Transport and shape calculus theorems

sn this setionD we present the trnsport nd shpe lulus theoremsD similrly to fernuer 8 rerzog PHIID tht re used in the ontinuous djoint derivtion @ghpter PAF Theorem A (Reynolds transport theorem) The derivative of the quantity

where V is the velocity eld in which the control volume Ω(t) moves.

Corollary A (Integration by Parts in Time in Moving Domains) For g = g(x, t) and h = h(x, t), we have ˆtf 0 ˆΩ(t)

Theorem B (Derivative of boundary integral) Let J(Ω) = ´Γ f ds be a boundary integral, the derivative is given by

where V is the velocity eld and κ is the mean curvature of Γ.

Appendix A. where w is the normal velocity of the moving surface S t and ḟ is the parameter-time

Corollary C (Integration by Parts in Time on a Moving Surface)

where g and h are restrictions of ĝ and ĥ to S t .

Cut Cell method for Navier-Stokes equations his setion presents the proposed disretiztion of the inompressile xvierE tokes equtions for n isotropi xewtonin )uid in the ontext of the gut gell methodF hetils of this disretiztion ! for sttionry geometries ! n e found in uirós odríguez et al. PHPPF he equtions red

where u nd p respetively denote the )uid9s veloity nd pressure (eldsD ρ its onstnt density nd g the grvittionl elertionF edditionllyD µ denotes the )uid9s onstnt dynmi visosity nd

the strinErte tensorF P = (P ij ) represents the @ellEenteredA pressure (eldD nd

the @stggeredA grtesin omponents of the veloity (eldF pinllyD D = (D x , D y ) denotes the @stggeredA oundry onditions to e pplied on the veloity (eldF Velocity divergence and pressure gradient vet Ω f ij = Ω ij ∩Ω f denote the suE set of Ω ij wet y the )uidD u the ontinuous )uid veloity (eld nd d the oundry onditionF henD tokes9 divergene theorem

sttes tht the volume integrl of the veloity divergene mthes the net volume )uxesD summed over the surfes immersed in the )uid itself nd djent to the oundryF he former termD referred to s homogeneousD qunti(es the exhnge of volume with the neighoring )uid elementsD nd the ltterD referred to s heterogeE neousD qunti(es this exhnge with the exterior domin through the oundryF his deomposition is re)eted t the disrete level y disretizing the volumeE integrted veloity divergene s

. @fFQA the sense tht the rte of hnge of the disrete kineti energy

is result of n exhnge with the neighoring )uid elements @(rst term in the rightEhnd side of iqution fFIHA nd ross the oundry @seond termAF Semi-discrete system he feEentered mss mtries ppering in front of the rte of hnge nd ody fores re digonl with oe0ients V = (V α ) @the volume of the stggered ontrol volumesD de(ned in etion QFIA nd re denoted s

qthering ll the termsD the proposed semiEdisrete momentum equtions then red @α ∈ {x, y}A

with divergeneEfree ondition cont (U, D) = 0. @fFIPA he system is losed with the disrete strinErte tensor SD de(ned s funtion of U nd D s followsD

where the opertors strain αβ re de(ned y iqutions fFQ nd fFRF ell of the opertors ppering in iqutions fFIID fFIP nd fFIQ re liner in ll dependent vriles @P D U nd SA nd oundry ondition D with the exeption of the onvetive trnsport opertors @(conv α ) de(ned in iqution fFTA whih is qudrti when evluted t U † = U nd D † = DF Projection method he disretiztion of the forementioned inompressE ile xvierEtokes equtions results in sddle point system of equE tions fenzi et al. PHHSD sometimes lso lled urushEuuhnEuker @uuA sysE tem xoedl 8 right PHHT in optimiztionF e wide rnge of lgorithms hve een devised to e0iently solve sddle point systems @or pproximtion thereofAF sn the (eld of )uid mehnisD ommon pproh is the frtionl step method ghorin IWTVF sn the present workD the method referred to s projetion method ss @mssA y frown et al. PHHID whih ensures seond order disretiztion of the equtionsD is employedF sn this projetion methodD the onvetive term is disretized using the expliit seondEorder edmsEfshforth sheme nd the visous term is disretized using the Appendix B. Cut Cell method for Navier-Stokes equations impliit grnkExiolson shemeF he (rst step of the method onsists of otining n intermedite veloity (eld U ⋆ y solving

where τ denotes the time step nd the supersript n the itertion numerF he oundry onditions pplile to U ⋆ @the predited veloity (eldA nd used in S ⋆ re those of the veloity (eld t the next time step

sn the projetion stepD the veloity (eld is updted y projeting U ⋆ using the intermedite pressure (eld Φ n+1 D whih is otined y solving the following oisson eqution τ cont pres Φ n+1 , 0 = cont U ⋆ , D n+1 , @fFISA with homogeneous xeumnn oundry onditions eing used for the intermedite pressure @0AF he veloity (eld is ultimtely orreted s

he pressure is (nlly updted s

where the lst term ensures the seond order ury of the pressure (eldF hus frD only hirihlet oundry onditions for the veloity (eld hve een onsideredD whih re pired with homogeneous oundry onditions for the pressure in the projetion stepF gses will e onsidered in the following setion where xeumnn oundry ondition re required long the out)ow oundriesF elong their viinityD hirihlet oundry ondition for the pressure is employed in order to uphold the omptiility eqution fFPF pinllyD the use of periodi ndGor xeumnn oundry onditions gives rise to rnk de(ieny in the vplin opertorF his results in the pressure (eld eing known up to onstntF his knowledge is exploited in the itertive solution of the oisson eqution y projeting the updtes in the spe of zeroEmen solutionsF Appendix C. Flower.jl package description xextD the meshes @gpD guD gv AD gut gell opertors in oth phses @opSD opLAD (elds in oth phses @phSD phLAD struture to sve the dt s the forwrd9 simultion progresses @fwd A nd the level set funtion @gp.phiA re initilizedX § ¤ ¦ ¥

e n now run the simultion using the min funtion run_forward(...)F sn this seD we solveX he het eqution @iqution IFIA in the solid phseF he tefn ondition @iqution IFPA t the interfeF he di'usionEonvetion eqution nd the xvierEtokes equtions with the foussinesq pproximtion @iqution IFIIA in the liquid phseF Appendix D Spurious currents study for static drops he presene of spurious urrents when studying )ows t low pillry numer re known issue in the gpEyp method opinet PHIVD s well s for other methods @level setD frontEtrkingAF xeverthelessD some improvements in the lne of surfe fores nd disrete pressure grdients re still possileF sn euEelEud et al. PHIVD the uthors presented new numeril sheme to model surfe tension for n interfe represented y levelEset funtion tht onserves )uid momentum nd reovers vple9s equilirium extlyF sn this setionD we look t the spurious urrents for stti drop in losed dominF he system is hrterized y the vple numerD representing the rtio of surfe tension to the visous dissiption

where D is hrteristi length sle of the systemD the dimeter of the drop in our seF imilrly to euEelEud et al. PHIVD we will look t the e'et of the vple numer nd the grid resolution on the spurious urrents y mesuring the mximl veloity s simultion proeeds in timeF sn ddition to thtD we will study the e'et of n initil shift in the xEdiretion ! therefore modifying the volume frtion in the interfil ells ! nd the e'et of 90 °ontt ngle imposed t the ottom wllF e will onsider T @2 × 3A ses @presented in pigure hFIA e drop initilized in the middle of the domin @noEonttA or in ontt with the sustrteF e initil shift in the xEdiretion @HD ∆/4 or ∆/2A with ∆ the grid sizeF sn eh seD we will vry he vple numer X v a PHHD RHHD VHHD IPHH he grid resolution X D/∆ a ITD QPD TR sn our simultionsD the drop is initilized with dimeter D = 1 in 2 × 2 dominF he densities of oth phses nd the surfe tension re set to IF he visosity @equl in oth phsesA will vry s

we n de(ne the disrete equivlent s

with u x [ghost] the tngentil veloity t the ghost ellD u x [ ] the tngentil veloity of the ell inside the dominD u x [1] the tngentil veloity of the ell one grid point wy from the wll nd ∆ the grid spingF e vlidte the numeril implementtion y onsidering oiseuille )ow in hnnel ounded y two wlls y ∈ [-0.5, 0.5] where we pply onstnt pressure grdient in the xEdiretionD suh tht

fy pplying the superEslip oundry ondition on oth wlls

we otin the following nlytil solution for the veloity (eld s funtion of y sn this etionD we present vlidtion of the of the implementtion of the free ontt ngle s de(ned in iqution IHFSF he ide is to otin y extrpoltion numeril ontt ngle θ ext t the wll using n pprent ngle θ app loted one ell ove the ontt lineF he dynmi ontt ngle θ d t the wll is then set to this extrpolted vlueF his method llows θ d to e solution of the )ow @insted of eing imposedAF elgorithm U summrizes the steps needed to ompute θ ext t given instntF Algorithm 7: engle extrpoltion pseudoEode for each cell do vote the ell one grid point ove the ontt line gompute the pprent ngle θ app using the unit norml n gompute the (rst order derivtive of the height funtion h x gompute the interfe urvture κ gompute the extrpolted ngle

end epply the extrpolted ngle t the wll through height funtions

prom kinemti onsidertionsD prike prike PHPH derived n nlytil solution for the trnsport of the ontt ngle in the se of n inompressile )owF sing this reltionD we vlidte our extrpoltion method in the yp frmework y rrying simultions of n osillting drop for di'erent grid sizesF e onsider drop of dimeter D = 1 in 2 × 2 domin tht is initilized over stti sustrte with ontt ngle θ 0 = 90 °F he veloity (eld in the whole dominD is set to Appendix G Relaxation of the impermeability condition in the toy model sn the toy model @etion IHFQAD omputtionl oundry is introdued t disE tne +δ from the wll where we pply xvier oundry onditionF hen onE sidering spreding dropD the impermeility ondition @u y | y=δ = 0A needs to e relxed in order to tke into ount the mss )ux through tht omputtionl oundryF o tht endD we introdue sink veloity v sink D suh tht

o lulte the mount of mss lost during one time stepD we ompute the re of the polygon formed y the two ontt line pointsD t times t n-1 nd t n D nd their liner projetion onto the tul wll @t y = 0AF sn prtieD the ngle used for the projetion is n pprent ngle loted one ell ove the ontt line @similrly to the extrpoltion ngle method desried in eppendix pAF he omputed re is then trnslted to the sink veloity y the simple formul

where τ is the time stepD L the length of the wll nd A the singed re ! de(ned lterF elgorithm V summrizes the required steps t given instntF e vlidte the method y onsidering drop spreding over sustrte for di'erent vlues of δF he drop is initilized with dimeter D = 1 in 2 × 2 domin nd we pply onstnt ontt ngle θ e = 70 °F he visositiesD densities nd the surfe tension re set to IF he simultions re run until the equilirium position is rehedF sn pigure qFID we show the totl mss s funtion of time when pplying the sink veloity for δ = 0.2, 0.05, 0.02F he initil mss inside the drop dereses s spreding oursF sn pigure qFPD we ompre the (nl interfes for three sesX @iA toy model with sink veloity @iiA toy model without sink veloity @iiiA full dominF prom the resultsD we n oserve the sink veloity llows us to mth to the (nl true equilirium shpe @full dominAF he error in (nl interfe is funtion of δF es δ is inresedD the liner extrpoltion proedure ! tht does not tke into ount the urved interfe t the ontt line ! indues higher di'erene in (nl shpesF sn prtieD when using the toy model on rel setupD δ is smll ompred to the domin sizeF Appendix G. Relaxation of the impermeability condition in the toy model