TN-07-1_ExcaliburNeptuneImplementation1DFluidSolverRealisticBoundaryConditions ============================================================================== .. meta:: :description: technical note :keywords: Excalibur-Neptune,report,2047356-TN-07-1,Task,2.2,BOUT++,1D,fluid,solver,with,realistic,boundary,conditions,:,implementation,Ben,Dudson,,Peter,Hill,,Ed,Higgins,,David,Dickinson,,and,Steven,Wright,University,of,York,David,Moxey,KCL,March,29,,2022,Contents,1,Executive,summary,2,Introduction,3,Implementation,of,selected,test,cases,3.1,Single,species,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,3.2,Fluid,neutrals,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,4,Conclusions,5,References,A,SD1D,Manual,1,1,3,3,4,5,7,7,1,Executive,summary,In,this,report,we,summarise,the,implementation,of,system,2-3,through,the,1D,model,SD1D,[1],including,fluid,neutrals,and,realistic,boundary,conditions.,We,also,discuss,selected,test,cases,from,task,83-2.1,implemented,with,SD1D.,This,provides,a,reference,implementation,of,a,1D,fluid,solver,with,realistic,boundary,conditions,for,use,in,other,work,,including,an,exploration,of,non-intrusive,UQ,in,task,83-2.4.,2,Introduction,Here,we,are,interested,in,an,implementation,with,BOUT++,[2],of,system,2-3.,This,system,is,described,in,detail,in,the,Equations,for,ExCALIBUR/NEPTUNE,Proxyapps,document,[3],so,here,we,provide,a,qualitative,summary,of,the,system,,rather,than,reproduce,the,full,system,of,equations.,The,system,under,consider-,ation,is,1D,,with,the,spatial,dimension,representing,the,distance,along,a,single,1,field,line.,This,is,typically,mapped,onto,the,region,from,the,outboard,midplane,(upstream),to,the,target,(downstream),,in,order,to,represent,the,behaviour,of,the,scrape,off,layer,(SOL).,The,equations,describing,this,plasma,evolution,are,derived,from,summing,Braginskii’s,two,fluid,system,alongside,a,small,number,of,assumptions.,The,end,result,is,that,we,represent,the,plasma,as,a,single,fluid,with,evolving,plasma,number,density,,pressure,(i.e.,energy,density),and,momentum.,In,addition,to,the,plasma,fluid,,the,system,also,includes,a,fluid,neutral,model.,The,neutral,system,takes,the,same,form,as,the,plasma,one,and,again,describes,the,evolution,of,number,density,,pressure,and,momentum.,The,neutral,and,plasma,fluids,are,then,coupled,through,sources,and,sinks,represent-,ing,the,various,plasma-neutral,interactions,such,as,recombination,,ionsiation,,charge,exchange,and,elastic,collisions.,The,1D,model,SD1D,[1],implements,such,a,system,,assuming,Te,=,Ti,,and,is,built,on,top,of,BOUT++,[2].,This,gives,SD1D,access,to,a,range,of,numeri-,cal,approaches,to,spatial,differencing,,time,integration,and,preconditioning,as,provided,by,BOUT++.,This,makes,it,a,useful,tool,for,experimenting,with,dif-,ferent,numerical,approaches.,The,SD1D,manual,,provided,alongside,the,code,,provides,a,brief,summary,of,the,physics,details,and,some,of,the,numerical,con-,siderations.,We,do,not,reproduce,this,here,and,for,convenience,,a,copy,of,the,manual,is,included,in,appendix,A.,The,user,has,a,lot,of,run,time,control,over,the,setup,of,the,system;,choosing,which,physics,terms,are,included,,the,value,of,physics,coefficients,,the,numerical,scheme,etc.,through,the,input,file.,We,shall,make,use,of,this,flexibility,in,task,83-2.4,when,uncertainty,quantification,(UQ),is,investigated.,When,applied,to,a,typical,divertor,system,the,appropriate,upstream,bound-,ary,conditions,are,generally,symmetry,,implemented,as,Neumann-zero,on,the,density,,pressure,and,temperature,(the,derivative,of,which,appears,in,the,heat,flux,,q),and,Dirichlet-zero,on,the,parallel,velocity.,At,the,downstream,,sheath,boundary,conditions,are,imposed,with,the,temperature,using,Neumann-zero,and,the,parallel,velocity,adopting,V(cid:107),≥,vs.,The,SD1D,implementation,handles,this,inequality,by,checking,the,value,of,V(cid:107),in,front,of,the,boundary,and,adopting,a,Dirichlet,condition,if,V(cid:107),<,vs,,switching,to,Neumann,if,V(cid:107),≥,vs,already.,2,3,Implementation,of,selected,test,cases,In,the,report,for,83-2.1,[4],a,number,of,generic,test,cases,for,1D,plasma,fluid,models,are,outlined.,Here,we,discuss,the,implementation,of,a,subset,of,these,in,SD1D,and,show,selected,results.,We,note,that,we,use,branch,bout-next,of,SD1D,at,commit,7bd6bc91,and,then,build,this,using,BOUT++,commit,080f3b27.,This,is,not,the,BOUT++,commit,which,will,be,used,by,default,if,SD1D,is,not,supplied,with,an,existing,BOUT++,build.,3.1,Single,species,Here,we,drop,the,neutral,fluid,entirely,and,only,evolve,the,single,plasma,fluid.,In,particular,we,focus,on,the,“Half,source,,sheath,boundary”,of,83-2.1,in,which,there,are,particle,and,energy,sources,distributed,over,the,upstream,half,of,the,domain,and,sheath,boundaries,are,enforced,at,the,target.,This,case,corresponds,to,example,case-03,shipped,with,SD1D.,We,also,note,that,this,is,setup,identi-,cally,to,case-02,of,SD1D,,except,it,includes,heat,conduction.,Whilst,this,does,not,drastically,modify,the,steady,state,obtained,,as,shown,in,figure,1,,this,can,have,a,significant,impact,on,the,numerical,performance.,For,example,case-02,is,found,to,run,in,around,5,seconds1,on,a,test,machine,,whilst,once,the,heat,condution,is,enabled,it,takes,∼,155,seconds.,We,note,that,case-03,actually,runs,in,around,7,seconds.,The,reason,this,is,not,as,slow,as,case-02,with,heat,condution,enabled,,is,that,it,also,enables,SD1D’s,custom,physics,based,precon-,ditioner,,which,attempts,to,solve,the,heat,conduction,problem,in,the,precondi-,tioner,stage.,This,highlights,the,importance,of,using,a,build,of,BOUT++,with,time,integrators,compatible,with,preconditioners,such,as,CVODE,and,PETSc.,One,can,also,disable,the,physics,motivated,preconditioner,and,instead,enable,a,generic,preconditioner,package,such,as,HYPRE.,For,example,,taking,case-03,and,disabling,the,physics,based,preconditioner,one,can,recover,a,run,time,of,7,seconds,by,switching,to,the,PETSc,based,time,integrator,,“beuler”,,and,using,either,PETSc’s,ILU,preconditioner,(serial,runs),or,HYPRE’s,euclid,(parallel,ILU),through,PETSc.,HYPRE’s,BoomerAMG,can,also,be,a,good,choice,when,running,in,parallel.,Tuning,the,preconditioner,parameters,can,lead,to,further,3,gains,in,more,expensive,test,cases.,(a),(b),Figure,1:,Comparison,of,the,plasma,state,variables,as,a,function,of,arc,length,,ρ,,at,the,final,time,for,case-02,and,case-03.,(c),3.2,Fluid,neutrals,Next,we,turn,our,attention,to,a,more,complete,test,case,which,now,includes,coupling,to,neutrals.,In,particular,,we,introduce,the,full,neutral,fluid,model,of,system,2-3.,The,simulation,setup,is,held,the,same,as,case-03,in,all,other,aspects,aside,from,a,reduction,in,the,recycling,fraction,from,0.9,to,0.2.,This,results,in,case-04,of,SD1D.,The,performance,of,this,case,is,again,rather,sensitive,to,the,preconditioning,approach,taken.,The,default,setup,when,BOUT++,is,con-,figured,with,CVODE,support,will,be,to,use,the,physics,based,preconditioner,,which,gives,a,run,time,of,53s,on,one,core.,Removing,this,preconditioner,sees,the,run,time,shoot,up,to,717,seconds.,Switching,to,the,PETSc,based,“beuler”,1A,large,fraction,(∼,40%),of,this,run,time,is,actually,in,I/O,overhead.,4,method,and,retaining,the,physics,based,preconditionter,sees,the,run,time,at,368,seconds.,If,one,instead,employs,ILU,rather,than,the,physics-based,precondi-,tioner,the,run,time,is,reduced,significantly,to,around,25,seconds.,The,result,of,running,this,test,case,is,shown,in,figure,2.,By,comparison,with,the,earlier,test,cases,one,can,see,that,the,plasma,structure,is,broadly,consistent,with,case-03,except,for,the,region,immediately,in,front,of,the,target,,as,one,would,expect.,We,also,see,that,the,default,resolution,,ny,=,200,,is,slightly,underresolved,when,compared,to,ny,=,600.,(a),(b),Figure,2:,Comparison,of,the,plasma,and,neutral,state,variables,as,a,function,of,arc,length,,ρ,,at,the,final,time,for,case-04.,(c),4,Conclusions,This,report,has,provided,a,brief,introduction,to,the,implementation,of,system,2-3,within,SD1D.,Further,details,are,provided,through,the,SD1D,manual,,in-,cluded,as,an,appendix.,In,addition,we,have,discussed,the,implementation,of,a,5,small,number,of,the,tests,outlined,in,83-2.1.,Whilst,the,SD1D,implementation,of,this,system,provides,a,very,convenient,tool,for,rapid,study,of,system,2-3,in,many,cases,it,is,perhaps,useful,to,note,that,a,yet,more,flexible,implemen-,tation,is,provided,in,Hermes-3,[5].,This,allows,for,the,relaxation,of,some,of,the,assumptions,made,here,(e.g.,Te,=,Ti),as,well,as,providing,a,pathway,to,consistently,extend,beyond,1D.,6,5,References,[1],Ben,Dudson,et,al.,SD1D,SOL,and,Divertor,model,in,1D.,https://github.,com/boutproject/SD1D.,[2],BOUT++,contributors.,BOUT++,manual.,https://bout-dev.,readthedocs.io/.,[3],Wayne,Arter,and,Benjamin,Dudson.,Equations,for,Ex-,CALIBUR/NEPTUNE,Proxyapps.,https://github.com/,ExCALIBUR-NEPTUNE/Documents/blob/main/reports/ukaea_reports/,CD-EXCALIBUR-FMS0021-1.20-M1.2.1.pdf.,[4],Benjamin,Dudson,,Peter,Hill,,Ed,Higgins,,David,Dickinson,,Steven,Wright,and,David,Moxey.,1D,fluid,model,tests.,https:,//github.com/ExCALIBUR-NEPTUNE/Documents/blob/main/reports/,2047356/TN-04.pdf.,[5],Ben,Dudson.,Hermes-3.,https://github.com/bendudson/hermes-3.,A,SD1D,Manual,A.1,Getting,started,First,get,a,copy,of,development,branch,of,BOUT++.,You,can,download,a,tarball,from,https://github.com/boutproject/BOUT-dev,,but,it,is,strongly,recommended,you,use,Git:,$,git,clone,https://github.com/boutproject/BOUT-dev.git,Configure,and,make,BOUT-dev,,including,SUNDIALS.,This,is,available,from,http://computation.llnl.gov/projects/sundials,,and,is,needed,for,pre-,conditioning,to,work,correctly.,$,cd,BOUT-dev,$,./configure,--with-sundials,$,make,7,The,user,manual,for,BOUT++,is,in,subdirectory,of,BOUT-dev,called,”manual”,,and,contains,more,detailed,instructions,on,configuring,and,compiling,BOUT++.,This,will,build,the,core,library,code,,which,is,then,used,in,each,model,or,test,case,(see,the,examples/,subdirectory),Next,download,a,copy,of,SD1D,into,the,BOUT-dev/examples,subdirectory.,This,isn’t,strictly,necessary,,but,it,makes,the,”make”,command,simpler,(other-,wise,you,add,an,argument,BOUT,TOP=/path/to/BOUT-dev/,to,make),BOUT-dev/examples/$,git,clone,https://github.com/boutproject/SD1D.git,BOUT-dev/examples/$,cd,SD1D,BOUT-dev/examples/SD1D,$,make,Hopefully,you,should,see,something,like:,Compiling,sd1d.cxx,Compiling,div_ops.cxx,Compiling,loadmetric.cxx,Compiling,radiation.cxx,Linking,sd1d,Here,the,main,code,is,in,”sd1d.cxx”,which,defines,a,class,with,two,methods:,init(),,which,is,run,once,at,the,start,of,the,simulation,to,initialise,everything,,and,rhs(),which,is,called,every,timestep.,The,function,of,rhs(),is,to,calculate,the,time,derivative,of,each,evolving,variable:,In,the,init(),function,the,evolving,variables,are,added,to,the,time,integration,solver,(around,line,192).,This,time,integration,sets,the,variables,to,a,value,,and,then,runs,rhs().,Starting,line,782,of,sd1d.cxx,you’ll,see,the,density,equation,,calculating,ddt(Ne).,Ne,is,the,evolving,variable,,and,ddt(),is,a,function,which,returns,a,reference,to,a,variable,which,holds,the,time-derivative,of,the,given,field.,BOUT++,contains,many,differential,operators,(see,BOUT-dev/include/difops.hxx),,but,work,has,been,done,on,improving,the,flux,conserving,Finite,Volume,imple-,mentations,,and,they’re,not,yet,in,the,public,repository.,These,are,defined,in,div,ops.hxx,and,div,ops.cxx.,The,atomic,rates,are,used,in,sd1d.cxx,starting,around,line,641,,and,are,defined,in,radiation.cxx,and,radiation.hxx.,8,To,run,a,simulation,,enter:,$,./sd1d,-d,case-01,This,will,use,the,”case-01”,subdirectory,for,input,and,output.,All,the,options,for,the,simulation,are,in,case-01/BOUT.inp.,The,output,should,be,a,whole,bunch,of,diagnostics,,printing,all,options,used,(which,also,goes,into,log,file,BOUT.log.0),,followed,by,the,timing,for,each,output,timestep:,Sim,Time,|,RHS,evals,|,Wall,Time,|,Calc,Inv,Comm,I/O,SOLVER,0.000e+00,5.000e+03,1.000e+04,1.500e+04,2.000e+04,2.500e+04,1,525,358,463,561,455,1.97e-02,1.91e-01,1.30e-01,1.68e-01,2.02e-01,1.65e-01,1.9,89.0,88.8,89.2,89.6,89.2,0.0,0.0,0.0,0.0,0.0,0.0,0.2,0.6,0.6,0.6,0.6,0.6,21.6,76.3,1.1,1.4,1.3,1.1,1.2,9.3,9.2,8.9,8.7,9.1,The,simulation,time,(first,column),is,normalised,to,the,ion,cyclotron,frequency,(as,SD1D,started,life,as,part,of,a,turbulence,model),,which,is,stored,in,the,output,as,”Omega,ci”.,So,each,output,step,is,5000,/,Omega,ci,=,104.4,mi-,croseconds.,The,number,of,internal,timesteps,is,determined,by,the,solver,,and,determines,the,number,of,times,the,rhs(),function,was,called,,which,is,given,in,the,second,column.,If,this,number,starts,steadily,increasing,,it’s,often,a,sign,of,numerical,problems.,To,analyse,the,simulation,,the,data,is,stored,in,the,”case-01”,subdirectory,along,with,the,input.,You,can,use,IDL,or,Python,to,look,at,the,”Ne”,,”NVi”,,”P”,variables,etc.,which,have,the,same,names,as,in,the,sd1d.cxx,code.,See,section,A.6,for,details,of,the,output,variables,and,their,normalisation.,The,evolving,variables,should,each,be,4D,,but,all,dimensions,are,of,size,1,except,for,the,time,and,parallel,index,(200).,Please,see,the,BOUT++,user,manual,for,details,of,setting,up,the,Python,and,IDL,reading,(”collect”),routines.,9,A.1.1,Examples,Case,1:,Without,heat,conduction,(Euler’s,equations),Removing,heat,conduction,reduces,the,system,to,fluid,(Euler),equations,in,1D.,Note,that,in,this,case,the,boundary,condition,(equation,4),is,subsonic,,because,the,adiabatic,fluid,sound,speed,is,cs,=,(cid:17)1/2,(cid:16),γp,n,γ,=,5/3,In,this,case,the,sources,of,particles,and,energy,are,uniform,across,the,grid.,Case,2:,Localised,source,region,The,same,equations,are,solved,,but,here,the,sources,are,only,in,the,first,half,of,the,domain,,applied,with,a,Heaviside,function,so,the,sources,abruptly,change.,Case,3:,Heat,conduction,We,now,add,Spitzer,heat,conduction,,the,κ||e,term,in,the,pressure,equation.,This,coefficient,depends,strongly,on,temperature,,and,severely,limits,the,timestep,unless,preconditioning,is,used.,Here,we,use,the,CVODE,solver,with,precondi-,tioning,of,the,electron,heat,flux.,In,addition,to,improving,the,speed,of,conver-,gence,,this,preconditioning,also,improves,the,numerical,stability.,Case,4:,Recycling,,neutral,gas,The,plasma,equations,are,now,coupled,to,a,similar,set,of,equations,for,the,neutral,gas,density,,pressure,,and,parallel,momentum.,A,fixed,particle,and,power,source,is,used,here,,and,a,20%,recycling,fraction.,Exchange,of,particles,,momentum,and,energy,between,neutrals,and,plasma,occurs,through,ionisation,,recombination,and,charge,exchange.,Case,5:,High,recycling,,upstream,density,controller,This,example,uses,a,PI,feedback,controller,to,set,the,upstream,density,to,1,×,1019m−3.,This,adjusts,the,input,particle,source,to,achieve,the,desired,density,,so,generally,needs,some,tuning,to,minimise,transient,oscillations.,This,is,controlled,by,the,inputs,10,density_upstream,=,1e19,density_controller_p,=,1e-2,density_controller_i,=,1e-3,The,input,power,flux,is,fixed,,specified,in,the,input,as,20MW/m2:,[P],#,Plasma,pressure,P,=,2,*,Ne,*,T,powerflux,=,2e7,#,Input,power,flux,in,W/m^2,The,recycling,is,set,to,95%,frecycle,=,0.95,NOTE:,This,example,is,under-resolved;,a,realistic,simulation,would,use,a,higher,resolution,,but,would,take,longer.,To,increase,resolution,adjust,ny:,ny,=,200,#,Resolution,along,field-line,Rather,than,200,,a,more,realistic,value,is,about,600,or,higher.,A.2,Plasma,model,Equations,for,the,plasma,density,n,,pressure,p,and,momentum,minV||i,are,evolved:,∂n,∂t,(cid:19),p,(cid:18),3,∂,∂t,2,(cid:0)minV||,∂,∂t,=,−∇,·,(cid:0)bV||n(cid:1),+,Sn,−,S,=,−∇,·,q,+,V||∂||p,+,Sp,−,E,−,R,(cid:1),−,∂||p,−,F,(cid:1),=,−∇,·,(cid:0)minV||bV||,j||,=,0,Ti,=,Te,=,1,2,p,en,q,=,5,2,pbV||,−,κ||e∂||Te,11,Which,has,a,conserved,energy:,(cid:90),V,(cid:20),1,2,minV,2,||i,+,(cid:21),p,dV,3,2,The,heat,conduction,coefficient,κ||e,is,a,nonlinear,function,of,temperature,Te:,κ||e,=,κ0T,5/2,e,where,κ0,is,a,constant.,See,section,A.8,for,details.,Operators,are:,∂||f,=,b,·,∇f,∇||f,=,∇,·,(bf,),(1),A.3,Boundary,conditions,A.3.1,Upstream:,Symmetry,Symmetry,boundary,conditions,are,applied,at,the,upstream,side,,corresponding,to,zero,flow,through,the,boundary.,∂||n,=,0,∂||p,=,0,∂||Te,=,0,V||,=,0,nV||,=,0,(2),Since,the,boundary,is,half-way,between,grid,points,,this,is,implemented,as,n0,=,n1,p0,=,p1,Te,0,=,Te,1,V||,0,=,−V||,1,nV||,0,=,−nV||,1,A.3.2,Downstream:,Sheath,Boundary,conditions,are,applied,to,the,velocity,and,the,heat,flux:,12,•,At,the,left,boundary,a,no-flow,condition,is,applied:,V||,=,0,∂||Te,=,0,•,At,the,right,boundary,is,a,sheath,boundary:,V||,≥,vs,∂||Te,=,0,where,the,inequality,is,implemented,by,switching,from,a,Dirichlet,to,a,Neumann,boundary,if,V||,>,vs,in,front,of,the,boundary.,The,critical,speed,into,the,sheath,,vs,is,sensitive,to,assumptions,on,the,thermodynamics,of,the,sheath,,taking,the,form:2,vs,=,(cid:18),e,(Te,+,γTi),mi,(cid:19)1/2,(3),where,Te,is,the,electron,temperature,(in,eV),,Ti,is,the,ion,temperature,,γ,is,the,ratio,of,specific,heats.,For,isothermal,flow,γ,=,1,,for,adiabatic,flow,with,isotropic,pressure,γ,=,5/3,,and,for,one-dimensional,adiabatic,flow,γ,=,3.,Here,we,are,assuming,Te,=,Ti,and,∂||Te,so,take,the,isothermal,case.,This,therefore,becomes:,vs,=,(cid:17)1/2,(cid:16),p,n,(4),Note:,If,the,sheath,velocity,is,subsonic,,then,waves,can,propagate,in,from,the,boundary.,Their,domain,of,dependence,is,outside,the,simulation,domain,,so,these,waves,can,cause,numerical,instabilities.,Several,boundary,conditions,are,available,for,the,density,and,pressure,,including,free,boundaries,and,Neumann,(zero,gradient).,These,are,controlled,by,settings,density,sheath,and,pressure,sheath.,Density,can,have,the,following,values:,0.,Free,boundary,,linearly,extrapolating,the,value,from,inside,the,domain,n−1,=,2n−2,−,n−3,2K-U,Riemann,,J.Phys.,D:Appl.,Phys,24,(1991),493-518,(5),13,1.,Neumann,(zero,gradient),n−1,=,n−2,2.,Constant,flux,n−1/2,=,n−2v−2J−2/,(cid:0)vsJ−1/2,(cid:1),(6),(7),where,the,Jacobian,factors,J,account,for,a,changing,flux,tube,cross-section,area.,Pressure,can,have,the,following,values:,0.,Free,boundary,,linearly,extrapolating,the,value,from,inside,the,domain,1.,Neumann,(zero,gradient),p−1,=,2p−2,−,p−3,p−1,=,p−2,(8),(9),2.,Constant,energy,flux,5,2,pv,+,1,2,nv3,5p−1/2,=,(cid:0)5p−2v−2,+,n−2v3,−2,(cid:1),/vs,−,n−1/2v2,s,(10),A.4,Sources,and,transfer,terms,External,sources,are,•,Sn,=,Source,of,plasma,ions,•,Sp,=,Source,of,pressure,,related,to,energy,source,SE,=,3,2,Sp,In,the,simulations,carried,out,so,far,,these,source,functions,are,both,constant,between,midplane,and,X-point,,and,zero,from,X-point,to,target.,A.4.1,Transfer,channels,There,are,several,transfer,channels,and,sinks,for,particles,,energy,and,momen-,tum,due,to,rates,of,recombination,,ionisation,,charge,exchange,,electron-neutral,14,excitation,,and,elastic,collisions,with,units,of,m−3s−1:,Rrc,=,n2,(cid:104)σv(cid:105)rc,Riz,=,nnn,(cid:104)σv(cid:105)iz,Rcx,=,nnn,(cid:104)σv(cid:105)cx,Rel,=,nnn,(cid:104)σv(cid:105)el,(Recombination),(Ionisation),(Charge,exchange),(Elastic,collisions),where,n,is,the,plasma,density;,nn,is,the,neutral,gas,density;,σcx,is,the,cross-,section,for,charge,exchange;,σrc,is,the,cross-section,for,recombination;,and,σiz,is,the,cross-section,for,ionisation.,Each,of,these,processes’,cross-section,depends,on,the,local,density,and,temperatures,,and,so,changes,in,time,and,space,as,the,simulation,evolves.,•,S,=,Net,recombination,i.e,neutral,source,(plasma,particle,sink).,Calcu-,lated,as,Recombination,-,Ionisation:,S,=,Rrc,−,Riz,•,R,=,Cooling,of,the,plasma,due,to,radiation,,and,plasma,heating,due,to,3-body,recombination,at,temperatures,less,than,5.25eV.,R,=,(1.09Te,−,13.6eV),Rrc,(Recombination),+,EizRiz,(Ionisation),+,(1eV),Rex,(Excitation),+,Rz,imp,(Impurity,radiation),The,factor,of,1.09,in,the,recombination,term,,together,with,factor,of,3/2,in,E,below,,is,so,that,recombination,becomes,a,net,heat,source,for,the,plasma,at,13.6/2.59,=,5.25eV.,Eiz,is,the,average,energy,required,to,ionise,an,atom,,including,energy,lost,through,excitation.,If,excitation,is,not,included,(excitation,=,false),then,following,Togo,et,al.,,Eiz,is,chosen,to,be,30eV.,If,excitation,is,included,,then,Eiz,should,be,set,to,13.6eV.,15,•,E,=,Transfer,of,energy,to,neutrals.,E,=,−,+,+,3,2,3,2,3,2,3,2,TeRrc,(Recombination),TnRiz,(Ionisation),(Te,−,Tn),Rcx,(Charge,exchange)**,(Te,−,Tn),Rel,(Elastic,collisions)**,(**),Note,that,if,the,neutral,temperature,is,not,evolved,,then,Tn,=,Te,is,used,to,calculate,the,diffusion,coefficient,Dn.,In,that,case,,Tn,is,set,to,zero,here,,otherwise,it,would,cancel,and,leave,no,CX,energy,loss,term.,•,F,=,Friction,,a,loss,of,momentum,from,the,ions,,due,to,charge,exchange,and,recombination.,The,momentum,of,the,neutrals,is,not,currently,mod-,elled,,so,instead,any,momentum,lost,from,the,ions,is,assumed,to,be,trans-,mitted,to,the,walls,of,the,machine.,F,=,miV||Rrc,(Recombination),−,miV||nRiz,+,mi,+,mi,(cid:0)V||,−,V||n,(cid:0)V||,−,V||n,(Ionisation),(cid:1),Rcx,(cid:1),Rel,(Charge,exchange),(Elastic,collisions),All,transfer,channels,are,integrated,over,the,cell,volume,using,Simpson’s,rule:,S,=,1,6JC,(JLSL,+,4JCSC,+,JRSR),where,J,is,the,Jacobian,of,the,coordinate,system,,corresponding,to,the,cross-,section,area,of,the,flux,tube,,and,subscripts,L,,C,and,R,refer,to,values,at,the,left,,centre,and,right,of,the,cell,respectively.,A.4.2,Recycling,The,flux,of,ions,(and,neutrals),to,the,target,plate,is,recycled,and,re-injected,into,the,simulation.,The,fraction,of,the,flux,which,is,re-injected,is,controlled,by,frecycle:,16,frecycle,=,0.95,#,Recycling,fraction,The,remaining,particle,flux,(5%,in,this,case),is,assumed,to,be,lost,from,the,system.,Note,that,if,there,are,any,external,particle,sources,,then,this,fraction,must,be,less,than,1,,or,the,number,of,particles,in,the,simulation,will,never,reach,steady,state.,Of,the,flux,which,is,recycled,,a,fraction,fredistribute,is,redistributed,along,the,length,of,the,domain,,whilst,the,remainder,is,recycled,at,the,target,plate,fredistribute,=,0.8,#,Fraction,of,recycled,neutrals,redistributed,evenly,along,length,The,weighting,which,determines,how,this,is,redistributed,is,set,using,redist,weight:,redist_weight,=,h(y,-,pi),#,Weighting,for,redistribution,which,is,normalised,in,the,code,so,that,the,integral,is,always,1.,In,these,expres-,sions,y,is,uniform,in,cell,index,,going,from,0,to,2π,between,the,boundaries.,The,above,example,therefore,redistributes,the,neutrals,evenly,(in,cell,index),from,half-way,along,the,domain,to,the,end.,When,neutrals,are,injected,,some,assumptions,are,needed,about,their,energy,and,momentum,•,When,redistributed,,neutrals,are,assumed,to,arrive,with,no,net,parallel,momentum,(so,nothing,is,added,to,N,Vn),,and,they,are,assumed,to,have,the,Franck-Condon,energy,(3.5eV,currently),•,When,recycled,from,the,target,plate,,neutrals,are,assumed,to,have,a,paral-,lel,momentum,away,from,the,target,,with,a,thermal,speed,corresponding,to,the,Franck-Condon,energy,,and,is,also,added,to,the,pressure,equation.,NOTE:,This,maybe,should,be,one,or,the,other,,but,not,both...,A.5,Neutral,model,The,number,of,equations,solved,is,controlled,by,the,following,parameters,in,the,input,file:,17,[NVn],evolve,=,true,#,Evolve,neutral,momentum?,[Pn],evolve,=,true,#,Evolve,neutral,pressure?,Otherwise,Tn,=,Te,model,Neutral,density,is,always,evolved,,so,turning,off,evolution,of,momentum,and,pressure,(setting,both,of,the,above,to,false),reduces,the,neutral,model,to,a,simple,diffusion,model,(next,section).,By,turning,on,the,momentum,equation,A.5.1,Diffusive,model,In,the,simplest,neutral,model,,neutral,gas,is,modelled,as,a,fluid,with,a,density,nn,which,diffuses,with,a,diffusion,coefficient,Dn:,∂nn,∂t,=,∇,·,(Dn∇nn),+,S,−,nn/τn,(11),The,temperature,of,the,neutrals,is,assumed,to,be,the,same,as,the,ions,Tn,=,Ti.Diffusion,of,neutrals,depends,on,the,neutral,gas,temperature,,and,on,the,collision,rate:,Dn,=,v2,th,n/,(νcx,+,νnn),(12),where,vth,n,=,(cid:112)eTn/mi,is,the,thermal,velocity,of,a,neutral,atom;,νcx,=,nσcx,is,the,charge-exchange,frequency,,and,σnn,=,vth,nnna0,is,the,neutral-neutral,collision,frequency,where,a0,(cid:39),π,(cid:0)5.29,×,10−11(cid:1)2,m2,is,the,cross-sectional,area,of,a,neutral,Hydrogen,atom.,In,order,to,prevent,divide-by-zero,problems,at,low,densities,,which,would,cause,D,to,become,extremely,large,,the,mean,free,path,of,the,neutrals,is,limited,to,1m.,An,additional,loss,term,is,required,in,order,to,prevent,the,particle,inventory,of,the,simulations,becoming,unbounded,in,detached,simulations,,where,recycling,no,longer,removes,particles,from,the,system.,This,represents,the,residence,time,for,neutral,particles,in,the,divertor,region,,which,in,[Togo,2013],was,set,to,around,10−4s.,18,A.5.2,Neutral,fluid,model,A,more,sophisticated,neutrals,model,can,be,used,,which,evolves,the,neutral,gas,momentum,and,energy:,∂nn,∂t,(cid:19),pn,(cid:18),3,∂,∂t,2,(cid:0)minV||n,∂,∂t,=,−∇,·,(cid:0)bV||nnn,(cid:1),+,∇,·,(Dn∇nn),+,S,−,nn/τn,=,−V||n∂||pn,+,∇,·,(κn∇Tn),+,∇,·,(DnTn∇nn),+,E,(cid:1),=,−∇,·,(cid:0)minV||nbV||n,(cid:1),−,∂||p,+,F,where,κn,is,the,neutral,gas,heat,conduction,coefficient.,This,is,assumed,to,be,κn,=,nnv2,th,n/,(νcx,+,νnn),i.e.,similar,to,Dn,for,the,diffusive,neutral,model,,but,with,a,factor,of,nn.,Note,that,if,the,diffusion,term,Dn,is,retained,in,the,neutral,density,(nn),equa-,tion,,then,a,corresponding,term,is,needed,in,the,pressure,(pn),equation.,To,re-,move,these,terms,,set,dneut,to,zero,in,the,input,options,,which,will,set,Dn,=,0.,The,density,diffusion,term,should,not,be,included,if,the,momentum,is,evolved,,and,so,is,switched,off,if,this,is,the,case.,The,continuity,equation,for,nn,is,exact,once,the,flow,is,known,,so,the,diffusive,flux,should,be,contained,in,the,flow,velocity,V||n.,To,see,where,this,comes,from,,assume,an,isothermal,neutral,gas:,∂nn,∂t,(cid:0)minV||n,∂,∂t,=,−∇,·,(cid:0)bV||nnn,(cid:1),+,S,−,nn/τn,(cid:1),=,−∇,·,(cid:0)minV||nbV||n,(cid:1),−,eTn∂||nn,+,F,Dropping,the,inertial,terms,reduces,the,momentum,equation,to,eTn∂||nn,=,F,=,νminn,(cid:0)V||i,−,V||n,(cid:1),where,ν,is,a,collision,frequency,of,the,neutrals,with,the,ions,,due,to,charge,exchange,,recombination,and,ionisation,(i.e.,νcx,+,νnn,as,used,in,the,calculation,19,of,diffusion,coefficient,Dn).,This,gives,an,equation,for,the,neutral,flow,velocity:,V||n,=,V||i,−,eTn,minnν,∂||nn,=,1,nn,v2,th,n,ν,∂||nn,where,vth,=,(cid:112)eTn/mi,is,the,neutral,thermal,speed,,as,used,in,the,calculation,of,Dn.,This,gives,a,flux,of,neutrals,nnV||n,=,nnV||i,−,Dn∂||nn,Hence,the,diffusive,flux,is,included,in,the,balance,between,pressure,gradients,and,friction,in,the,momentum,equation.,A.6,Outputs,Output,quantities,are,normalised,,with,the,normalisation,factors,stored,in,the,output,files,Table,1:,Normalisation,quantities,Name,Nnorm,Tnorm,Cs0,Omega,ci,rho,s0,Description,Density,Temperature,Speed,Time,Length,Units,m−3,eV,m/s,1/s,m,The,following,variables,are,stored,in,the,output,file,if,they,are,evolved:,Name,Description,Ne,NVi,P,Nn,NVn,Pn,Plasma,density,Plasma,flux,Plasma,pressure,Neutral,density,Neutral,flux,Neutral,pressure,Normalisation,Nnorm,[m−3],Nnorm×Cs0,[m−2s−1],e×Nnorm×Tnorm,[Pascals],Nnorm,[m−3],Nnorm×Cs0,[m−2s−1],e×Nnorm×Tnorm,[Pascals],The,following,rates,and,coefficients,are,also,stored:,20,Name,Description,Sink,of,plasma,density,Normalisation,Nnorm×Omega,ci,[m−3s−1],Sink,of,plasma,momentum,mi×Nnorm×Cs0×Omega,ci,[Nm−3],e×Nnorm×Tnorm×Omega,ci,[Wm−3],e×Nnorm×Tnorm×Omega,ci,[Wm−3],Radiative,loss,of,energy,Sink,of,plasma,energy,S,F,R,E,kappa,epar,Plasma,thermal,conduction,Neutral,diffusion,coefficient,Dn,flux,ion,Flux,of,ions,to,target,Note,that,the,R,term,is,energy,which,is,lost,from,the,system,,whilst,E,is,energy,which,is,transferred,between,plasma,and,neutrals.,For,all,transfer,terms,(S,,F,,R),a,positive,value,means,a,transfer,from,plasma,to,neutrals.,To,diagnose,atomic,processes,,turn,on,diagnose,=,true,in,the,input,settings,(this,is,the,default).,Additional,outputs,contain,the,contributions,from,each,atomic,process.,They,have,the,same,normalisation,factors,as,the,corresponding,(S,,F,,R),term.,Name,Srec,Siz,Frec,Fiz,Fcx,Fel,Rrec,Riz,Rzrad,Rex,Erec,Eiz,Ecx,Eel,Description,Sink,of,plasma,particles,due,to,recombination,Sink,of,plasma,particles,due,to,ionisation,(negative),Sink,of,plasma,momentum,due,to,recombination,Sink,of,plasma,momentum,due,to,ionisation,Sink,of,plasma,momentum,due,to,charge,exchange,Sink,of,plasma,momentum,due,to,elastic,collisions,Radiation,loss,due,to,recombination,Radiation,loss,due,to,ionisation,(inc.,excitation),Radiation,loss,due,to,impurities,Radiation,loss,due,to,electron-neutral,excitation,Sink,of,plasma,energy,due,to,recombination,Sink,of,plasma,energy,due,to,ionisation,Sink,of,plasma,energy,due,to,charge,exchange,Sink,of,plasma,energy,due,to,elastic,collisions,21,A.7,Atomic,cross,sections,Cross,sections,are,approximated,with,semi-analytic,expressions,,obtained,from,E.Havlickova,but,of,unknown,origin.,For,the,purposes,of,calculating,these,cross-sections,,any,temperatures,below,1eV,are,set,to,1eV.,The,charge,exchange,cross-section,is,approximated,as:,(cid:40),σiz,=,10−14T,1/3,10−14,if,T,≥,1eV,if,T,<,1eV,(13),with,units,of,[m3/s].,Ionisation,is,calculated,as,σcx,=,,,,5.875,×,10−12,·,T,−0.5151,·,10−2.563/,log10,T,10−6,·,T,−3.054,·,10−15.72,exp(−,log10,T,)+1.603,exp(−,log2,7.638,×,10−21,10,T,),if,T,≥,20eV,if,1eV,<,T,<,20eV,if,T,≤,1eV,(14),Recombination,rates,are,calculated,using,a,9,×,9,table,of,coefficients,so,is,not,reproduced,here.,Figure,3:,Cross-sections,[Thanks,to,E.Havlickova,and,H.Willett],Plots,of,these,cross-sections,are,shown,in,figure,3.,There,are,a,few,anomalies,with,this:,charge,exchange,always,has,the,highest,cross-section,of,any,process,,and,ionisation,has,a,jump,at,20eV.,The,ionisation,and,charge,exchange,rates,do,not,depend,on,density,,but,recombination,does,so,a,typical,range,of,values,22,051015202530Electron,temperature,[eV]10-2010-1910-1810-1710-1610-1510-1410-1310-12Rate,<σv>,[m3s−1]IonisationRecombination,(n,=,1018m−3)Recombination,(n,=,1020m−3)Recombination,(n,=,1022m−3)Charge,exchange,is,shown.,A.8,Heat,conduction,Spitzer,heat,conduction,is,used,κ||e,=,3.2,ne2T,τe,me,(cid:39),3.1,×,104,T,5/2,ln,Λ,(15),which,has,units,of,W/m/eV,so,that,in,the,formula,q,=,−κ||e∇T,,,q,has,units,of,Watts,per,m2,and,T,has,units,of,eV,.,This,uses,the,electron,collision,time:,√,6,τe,=,√,2π3/2(cid:15)2,0,ln,Λe2.5n,meT,3/2,e,(cid:39),3.44,×,1011,T,3/2,e,ln,Λn,in,seconds,,where,T,e,is,in,eV,,and,n,is,in,m−3.,Normalising,by,the,quantities,in,table,1,gives,ˆκ||e,=,3.2ˆn,ˆTe,mi,me,τeΩci,where,hats,indicate,normalised,(dimensionless),variables.,A.9,Non-uniform,mesh,(16),(17),An,example,of,using,a,non-uniform,grid,is,in,diffusion,pn.,The,location,l,along,the,field,line,as,a,function,of,normalised,cell,index,y,,which,goes,from,0,at,the,upstream,boundary,to,2π,at,the,target,,is,(cid:20),(2,−,δymin),l,=,L,y,2π,−,(1,−,δymin),(cid:17)2(cid:21),(cid:16),y,2π,(18),where,0,<,δymin,<,1,is,a,parameter,which,sets,the,size,of,the,smallest,grid,cell,,as,a,fraction,of,the,average,grid,cell,size.,The,grid,cell,spacing,δy,therefore,varies,as,δy,=,L,Ny,(cid:104),1,+,(1,−,δymin),(cid:16),1,−,(cid:17)(cid:105),y,π,(19),This,is,set,in,the,BOUT.inp,settings,file,,under,the,mesh,section:,23,dy,=,(length,/,ny),*,(1,+,(1-dymin)*(1-y/pi)),In,order,to,specify,the,size,of,the,source,region,,the,normalised,cell,index,y,at,which,the,location,l,is,a,given,fraction,of,the,domain,length,must,be,calculated.,This,is,done,by,solving,for,y,in,equation,18.,yxpt,=,π,(cid:20),2,−,δymin,−,(cid:113),(2,−,δymin)2,−,4,(1,−,δymin),fsource,(cid:21),/,(1,−,δymin),(20),which,is,calculated,in,the,BOUT.inp,file,as,y_xpt,=,pi,*,(,2,-,dymin,-,sqrt(,(2-dymin)^2,-,4*(1-dymin)*source,),),/,(1,-,dymin),where,source,is,the,fraction,fsource,of,the,length,over,which,the,source,is,spread.,This,is,then,used,to,calculate,sources,,given,a,total,flux.,For,density:,source,=,(flux/(mesh:source*mesh:length))*h(mesh:y_xpt,-,y),which,switches,on,the,source,for,y,<,yxpt,using,a,Heaviside,function,,then,divides,the,flux,by,the,length,of,the,source,region,fsourceL,to,get,the,volumetric,sources.,A.10,Numerical,methods,All,variables,are,defined,at,the,same,location,(collocated).,Several,different,numerical,methods,are,implemented,,to,allow,testing,of,their,accuracy,and,robustness.,A.10.1,Advection,terms,∇,·,(cid:0)bV||f,(cid:1),Flux,splitting,,MinMod,limiter,The,default,method,uses,a,combination,of,HLL-style,flux,splitting,and,MinMod,slope,limiting.,Terms,of,the,form,∇,·,(bf,),are,implemented,as,fluxes,through,24,cell,boundaries:,∇,·,(bV,f,)i,(cid:39),1,J∆y,(cid:2)Fi+1/2,−,Fi−1/2,(cid:3),(21),where,F,is,the,flux.,This,is,calculated,by,linearly,interpolating,the,velocity,to,the,cell,edges,Vi+1/2,=,1,2,(Vi,+,Vi+1),(22),The,field,being,advected,,f,,,is,reconstructed,from,the,cell,centre,values,fi,onto,cell,edges,f,L,i,and,f,R,i,:,f,L,i,=,fi,−,1,2,s,f,R,i,=,fi,+,1,2,s,where,the,slope,s,is,limited,using,the,MinMod,method:,s,=,,,,0,fi+1,−,fi,fi,−,fi−1,if,sign(fi+1,−,fi),(cid:54)=,sign(fi,−,fi−1),if,|fi+1,−,fi|,<,|fi,−,fi−1|,otherwise,(23),(24),In,order,to,handle,waves,travelling,both,left,and,right,,flux,splitting,handles,characteristics,moving,left,differently,from,characteristics,moving,right.,In,gen-,eral,this,is,problem,dependent,and,computationally,expensive,,so,here,we,adopt,a,simple,approximation,similar,to,an,HLL,splitting3.,We,assume,that,the,fastest,waves,in,the,system,travel,with,speed,a,(the,sound,speed),with,respect,to,the,flow,,so,that,there,are,waves,travelling,with,V,+,a,and,V,−,a.,If,the,flow,speed,is,supersonic,then,these,waves,are,only,in,one,direction,,but,for,subsonic,flows,there,is,a,flux,in,both,directions.,The,fluxes,between,cells,are,calculated,using:,Fi+1/2,=,,,,f,R,i,Vi+1/2,f,L,i+1Vi+1/2,f,R,i,1,2,(cid:0)Vi+1/2,+,a(cid:1),+,f,L,i+1,if,Vi+1/2,>,a,if,Vi+1/2,<,−a,otherwise,(25),(cid:0)Vi+1/2,−,a(cid:1),1,2,(cid:1),,(cid:0)f,R,i,−,f,L,Hence,for,subsonic,flows,the,flux,becomes,Vi+1/2,i+1,i,(cid:39),f,L,where,the,second,term,is,a,diffusion.,When,the,solution,is,smooth,,f,R,i+1,,the,numerical,method,becomes,central,differencing,and,the,diffusion,goes,to,zero,as,∆x2.,Oscillatory,solutions,introduce,dissipation,,and,the,method,becomes,i,+,f,L,i+1,(cid:1),+,a,(cid:0)f,R,1,2,2,3A.,Harten,,P.,D.,Lax,,and,B.,van,Leer,”On,Upstream,Differencing,and,Godunov-Type,Schemes,for,Hyperbolic,Conservation,Laws”,,SIAM,Review,,25(1),,pp.,35-61,,1983,25,increasingly,upwind,as,the,flow,becomes,sonic.,Nonlinear,fluxes,When,advecting,quantities,which,are,a,nonlinear,combination,of,variables,,such,as,nV||,,conservation,properties,can,be,slightly,improved,by,using,the,following,interpolation4,5,6:,(f,g)R,=,(cid:0)f,RgC,+,f,CgR(cid:1),(26),1,2,where,superscript,C,denotes,cell,centre,,and,R,right,hand,side.,This,method,is,implemented,,using,MinMod,interpolation,for,each,variable.,Central,differencing,Central,difference,schemes,have,an,advantage,over,upwind,schemes,,in,that,they,do,not,need,to,take,account,of,wave,speeds.,The,simple,central,differencing,scheme,produces,large,unphysical,oscillations,,due,to,the,decoupling,of,odd,and,even,points,in,collocated,schemes,,but,can,(usually),be,stabilised,by,adding,dissipation.,It,is,implemented,here,for,comparison,with,other,schemes.,Skew,symmetric,central,differencing,A,simple,modification,to,the,central,differencing,scheme,improves,numerical,stability,,coupling,nearby,points7,8,The,idea,is,to,split,the,divergence,terms,into,a,“skew-symmetric”,form,∇,·,(cid:0)bV||f,(cid:1),=,1,2,(cid:2)∇,·,(cid:0)bV||f,(cid:1),+,V||b,·,∇f,+,f,∇,·,(cid:0)bV||,(cid:1)(cid:3),(27),Each,of,the,terms,on,the,right,are,then,discretised,with,standard,2nd-order,central,differences.,This,method,can,avoid,the,need,for,additional,dissipation,,or,be,stabilised,with,a,smaller,viscosity,than,the,simple,central,differencing,method.,4F.N.Felten,,T.S.Lund,“Kinetic,energy,conservation,issues,associated,with,the,collocated,mesh,scheme,for,incompressible,flow”,J.Comp.Phys.,215,(2006),465-484,5F.N.Felten,,T.S.Lund,“Critical,comparison,of,the,collocated,and,staggered,grid,arrange-,ments,for,incompressible,turbulent,flow”,Report,ADP013663,6Y.Morinishi,et,al.,“Fully,Conservative,Higher,Order,Finite,Difference,Schemes,for,In-,compressible,Flow”,J.Comp.Phys.,143,(1998),90-124,7S.Pirozzoli,“Stabilized,non-dissipative,approximations,of,Euler,equations,in,generalized,curvilinear,coordinates”,J.Comp.Phys.,230,(2011),2997-3014,8A.E.Honein,,P.Moin,“Higher,entropy,conservation,and,numerical,stability,of,compressible,turbulence,simulations”,J.Comp.Phys.,201,(2004),532-545,26,A.10.2,Artificial,viscosity,Artificial,viscosity,(viscos,input),is,implemented,as,a,diffusion,of,momentum,in,index,space,,so,that,the,diffusion,coefficient,varies,as,∆y2.,∂,∂t,(cid:0)nV||,(cid:1),i,=,.,.,.,+,ν,(cid:2)(Vi+1,−,Vi),Ji+1/2,−,(Vi,−,Vi−1),Ji−1/2,(cid:3),/Ji,(28),where,J,is,the,Jacobian,,subscript,i,indicates,cell,index,,and,Ji+1/2,=,(Ji,+,Ji+1),/2.,27 :pdfembed:`src:_static/TN-07-1_ExcaliburNeptuneImplementation1DFluidSolverRealisticBoundaryConditions.pdf, height:1600, width:1100, align:middle`