TALK=T;RUN(1,1)
************************************
PQ1 description
************************************
Title: Tubeflow test
Author: Elena Pankova
Start Date: 22.02.2013
Editor: Brian Spalding
for the purposes of:
(1) removing split of tube into five sections
(2) increasing the number of fluids
(3) reducing the number of turbulence models
(4) removing out-of-date comments and deactivations
(5) introducing explanatory comments
(6) removal of ydimrat
(7) providing print options for inforout only
(8) utilising new inforin capability
(9) removing the parabolic optio entirel, in the grounds that
* it is not usable only for incompressile fluids,
* even for then it provides only an economy of computer
memory which users will not notice; and
* its nature would have to be explained if it were retained.
(10) place flow-configuration option in 'general' rather them
'other numerical'
(11) toend and howmany introduced
Last editing date: 23.03.13
Purpose: To test the 10-item structure
Brief description:
Flow is through a straight pipe of circular cross-section with thick
wall allowing heat exchange with an external environment.
Calculation modes are:
1. fully-developed, with uniform longitudinal temperature gradient
2. fully-developed, with fixed wall temperature
3. elliptic, i.e whole tube, but not parabolic
************************************
PQ1 Part 1 Declarations and settings
************************************
vredit=f ! command preventing satellite from double reading
! and re-writing of Q1. It renders unnecessary the
! introduction of saveXbegin and saveXend around In-Form
! statements when VR objects are present.
Utilities
char(formula)
real(dummy)
xml-group general
insert declaration and settings below
CHAR(floform) ! flow formulation
floform=developing_flow! full_dev_uniform-Tgrad; full_dev_uniform-T$
wall
xml-group geometry
insert declaration and settings below
real(xang) ! circumferential extent, radians
xang=3.14159
real(diam) ! tube inner diameter, m
diam=0.05
! command renew
! dbs 13.03.13 An explanatory note on this is needed
real(thck) ! wall thickness, m
thck=0.002
! command renew
real(lng1) ! tube length, m
lng1=5.0
! dbs 13.03.12 one length segment only
xml-group variables solved
insert declaration and settings below
solution
boolean(solp1) ! solve p1
solp1=T ! F
boolean(solu1) ! solve u1
solu1=T ! F
boolean(solv1) ! solve v1
solv1=T ! F
boolean(solw1) ! solve w1
solw1=T ! F
boolean(solt1) ! solve tem1
solt1=T ! F
xml-group material properties
insert declaration and settings below
char(flname) ! name of fluid in case 089
flname=use_specified_properties ! saturated_water;$
SAE_5W-30_engine_oil;ethylene_glycol;glycerin;air;$
ammonia;carbon_dioxide;carbon_monoxide;hydrogen;$
nitrogen;oxygen;superheated_water_vapour
real(prno) ! Prandtl number
prno=7.0
fluid values of water at 20 degC as in PROPS file
real(dnsf) ! fluid density, kg/m**3
dnsf=998.23
dnsf=1.0 !****
!**** indicates values sometimes introduced via manual
! editing during testing
real(sphf) ! fluid specific heat, J/(kg*degC)
sphf=4181.8
sphf=1.0 !****
real(nulf) ! fluid kinematic viscosity, m**2/s
nulf=1.006E-6
nulf=1.e-5 !****
real(cndf) ! fluid thermal conductivity, W/(m*degC)
cndf=sphf*nulf*dnsf/prno
real(prno) ! Prandtl number sphf*dnsf*nulf/cndf
prno=7.0
cndf=0.597
cndf=0.025 !****
real(texf) ! fluid thermal expansient coeff. degC**-1
texf=1.18e-4
dnsf=1
sphf=1
cn
! tube-wall values for copper at 27 degC
real(dnst) ! tube wall density, kg/m**3
dnst=8954.0
real(spht) ! tube wall specific heat, J/(kg*degC)
spht=383.0
real(cndt) ! tube wall thermal conductivity, W/(m*degC)
cndt=381.0
xml-group models
insert declaration and settings below
char(turb) ! turbulence model
! dbs 13,02.13 de-activate kl and ko models
turb=lvel!none;kemodl
turb=none
turb=kemodl
boolean(gravty) ! gravity
gravty=F ! T
gravty=t !****
real(gravax) ! xcart-dir gravitational acceleration, m/s**2
gravax=9.81
real(gravay) ! ycart-dir gravitational acceleration, m/s**2
gravay=0.0
real(gravaz) ! zcart-dir gravitational acceleration, m/s**2
gravaz=0.0
xml-group initial conditions
insert declaration and settings below
xml-group boundary conditions
insert declaration and settings below
char(winset) ! inlet velocity setting
winset=from_Re!from_set_value
real(reno) ! Reynolds number
reno=5000
char(wfin) ! inlet velocity, m/s
wfin=1.0
real(tfin) ! inlet temperature, degC
tfin=80.0
tfin=100.0 !***
tfin=21.0 !***
real(pout) ! outlet gauge pressure, bar
pout=1.e5
pout=0.0
real(tenv1) ! external temperature section 1, degC
tenv1=20.0
! dbs 13.03.13 de-activate 5-length-segment provision
real(envco1) ! external heat-tr, coeff. sect.1, W/(m**2*degC)
envco1=1.e+4
! dbs 13.03.13 de-activate 5-length-segment provision
real(foures) ! fouling resistance m**2 degC/W
foures=0.0
envco1=1.e+7 !*****
! note that the next sequence has the format which enables pq1ed
! to create the appropriate list lines in tubeflow.xml
! The ! following tgrad therefore does NOT signify that a comment
! follows. Instead alternative settings follow, separated by
! semi-colons
xml-group output
insert declaration and settings below
Print to inforout
integer(howmany) ! how many output sets?
howmany=1
! dbs 22.03.13 print-out features relating to RESULT file
! deactivated because these are unlikely to interest SimScene
! users
boolean(prgen) ! print general
prgen=T ! F
boolean(prgrge) ! grid and geometry
prgrge=T ! F
boolean(prrepr) ! ref props
prrepr=T ! F
boolean(prdino) ! dimensionless nos.
prdino=T ! F
boolean(prboco) ! boundary conds.
prboco=T ! F
boolean(prtumo) ! turbulence models
prtumo=T ! F
boolean(prnum) ! numerical
prnum=T ! F
boolean(prntall) ! print all?
prntall=F ! T
xml-group computational grid
insert declaration and settings below
boolean(polar) ! indicator of polar grid
polar=T ! F
polar=f !***
integer(nxx) ! circumferential interval no.
nxx=36
!!! reduced to 1 when absence of gravity render flow 2D
nxx=1 !******
integer(nrgy) ! radial regions
nrgy=2
integer(nrgy1) ! interval no. in radial region 1
nrgy1=20
integer(nrgy2) ! interval no. in radial region 2
nrgy2=2
! dbs 13.03.13 This option now withdrawn for reason below
real(ydimrat) ! y-cell diminution ratio
ydimrat=1.0
! Note that 0.9 was at first made the default; but when
! nrgy1 was made large, this led to errors resulting
! from excessively small grid intervals. It seems to be
! a dangerous grid option for that reason.
integer(nrgz) ! number of longitudinal regions
! dbs 13.03.13 now one region only
nrgz=1
integer(nrgz1) ! longitudinal interval no 1
nrgz1=50
xml-group numerical
insert declaration and settings below
integer(numits) ! number of iterations
numits=100
real(relxw1) ! relaxation factor for longitudinaal velocity
relxw1=0.5
real(relxt1) ! relaxation factor for temperature
relxt1=0.5
boolean(toend) ! continue run to end
toend=T!F
xml-end
************************************
PQ1 Part 2 Include Frommenu
************************************
INCL(frommenu.htm)
mesg(after frommenu
floform
************************************
PQ1 Part 3 Consequential settings
************************************
GROUP 1. Run title and other preliminaries
Text(tubeflow)
title
floform
BOOLEAN(fdf,tgrad,twall)
if(:floform:.eq.full_dev_uniform_tgrad) then !-----------
fdf=t !
tgrad=t !
numits=50*numits !
nrgy1=nrgy1*4 ! this increased grid fineness seems to be !
! necessary for accurate Nusselt no, !
mesg(numits and nrgy1 increased by satellite !
endif !-----------
if(:floform:.eq.full_dev_uniform_twall) then !-----------
fdf=t !
tgrad=t !
numits=50*numits !
nrgy1=nrgy1*4 ! this increased grid fineness seems to be !
! necessary for accurate Nusselt no, !
mesg(numits and nrgy1 increased by satellite !
endif !-----------
if(:floform:.eq.developing_flow) then
fdf=f
endif
fdfsol=fdf
TEXT(fluid is :flname: flow is :floform: grav=:gravty:
gravty
if(.not.gravty) then
nxx=1
xang=0.1
numits=100 !*****
solu1=f
if(fdf) then
solp1=f ! no need to solve for pressure when nx=1
endif
endif
mesg(fdf=:fdf:
if(foures.gt.0.0) then ! add fouling resistance
do ii=1,5
envco:ii:=1/(foures+1.0/envco:ii:)
enddo
endif
************************************************************
IRUNN = 1
GROUP 2. Transience; time-step specification
STEADY = T
************************************************************
GROUP 3. X-direction grid specification
CARTES = .not.polar
cartes=t !****
NX = 1
XULAST =1.
XFRAC(1)=1.
grdpwr(x,nxx,xang,1.0)
nx
xulast
************************************************************
GROUP 4. Y-direction grid specification
NY = 1
YVLAST =1.
YFRAC(1)=1.
nregy=nrgy
iregy=1;grdpwr(y,nrgy1,-0.5*diam,ydimrat)
! dbs 13.03.13 ydimrat replaced by 1.0
iregy=1;grdpwr(y,nrgy1,-0.5*diam,ydimrat)
iregy=1;grdpwr(y,nrgy1,-0.5*diam,0.95) !*****
iregy=1;grdpwr(y,nrgy1,-0.5*diam,1.0) !*****
yvlast
iregy=2;grdpwr(y,nrgy2,thck,1.0)
ny
yvlast
************************************************************
GROUP 5. Z-direction grid specification
NZ = 1
ZWLAST =1.
ZFRAC(1)=1.
if(fdf) then ! fully-developed flow
nrgz1=1;lng1=0.1*diam
grdpwr(z,nrgz1,lng1,1.0)
ydimrat=1.0
nz
else ! elliptic flow
nregz=nrgz
iregz=1;grdpwr(z,nrgz1,lng1,1.0)
! dbs 13.03.13 5 z regions reduced to 1
endif ! end of first if(fdf)
nz
zwlast
************************************************************
GROUP 6. Body-fitted coordinates or grid distortion
************************************************************
GROUP 7. Variables stored, solved & named
ONEPHS = T
* Y in SOLUTN argument list denotes:
* 1-stored 2-solved 3-whole-field
* 4-point-by-point 5-explicit 6-harmonic averaging
STORE(A) ! Variable added to show that cn149 changes to cn148
if(solp1) then
solve(p1)
if(fdf) then
solutn(p1,y,n,n,n,n,n) ! store but do not solve
endif
endif
if(solu1) then
solve(u1)
endif
if(solv1) then
solve(v1)
if(fdf) then
solutn(v1,y,n,n,n,n,n) ! store but do not solve
endif
endif
if(solw1) then
solve(w1)
if(fdf) then
terms(w1,p,n,p,p,p,p)
endif
endif
! Important note: the order in shich the solve and ator staements
! appear is very important because thos other than pressures,
! velocities, enhalpies , ke and eps are given variable indices
! of which start large (usually 150) and then diminish; but the
! lower-index variables are calsolated before higher-index ones
! in earth.
! it is therefore important that all property variables which
! depend depend upon temeperature should be stored before tem1
! in order that temerature is calculated before they are; since,
! for gaces the enu_expression contains rho1, enul should be
! stored before rho1, so that rho1 is calculated before it.
! Otherwise initial values will be incorrect.
store(enul,cp1,rho1)
store(cond)
if(solt1) then
solve(tem1)
terms(tem1,n,p,p,p,p,p) ! de-activate kinetic heating
endif
store(shrz) ! needed for prrad etc, and for nearwall profil
if(prskin) then ! various storage settings for print-out purposes
store(skin)
endif
store(skin) ! needed until e1start.for of 12.06.12 reached !
! earth delivery !
if(prstan) then
store(stan)
endif
if(prhtco) then
store(htco)
endif
if(:turb:.ne.none) then
turmod(:turb:)
store(enut,len1)
if(:turb:.eq.lvel) then
solve(ltls)
store(wdis,wgap)
endif
if(:turb:.eq.kemodl) then
if(fdf) then
terms(ke,p,n,p,p,p,p)
terms(ep,p,n,p,p,p,p)
endif
endif
endif
! dbs 13.03.13 kl and ko models removed
store(prps)
************************************************************
mesg(group 8
GROUP 8. Terms (in differential equations) & devices
* Y in TERMS argument list denotes:
* 1-built-in source 2-convection 3-diffusion 4-transient
* 5-first phase variable 6-interphase transport
DIFCUT =0.5 ;ZDIFAC =1.
GALA = F ;ADDDIF = F
ISOLX = -1 ;ISOLY = -1 ;ISOLZ = -1
************************************************************
GROUP 9. Properties of the medium (or media)
press0=1.e5 ! additive to p1 to create absolute pressure
temp0=273 ! additive to p1 to create absolute pressure
! note that the following declarations appear in
! always called-at-the-start library case 014.htm
CHAR(ABT,ABP) ! absolute temperature and pressure for 089
ABT=(TEM1+TEMP0); ABP=(P1+PRESS0) ! default setting
mesg(flname is :flname:
condtn=if(iy.le.:nrgy1:)
condtn
(property enul is 0.0!if(prps.gt.0))
if(:flname:.eq.use_specified_properties) then
mesg(use specified properties
! the following are provided to enable the boundary-condition
! settings to be the same whether or not constant properties
! are being used
! declarations
char(rho_expression)
char(cp_expression)
char(enu_expression)
char(cond_expression)
! settings
rho_expression=dnsf
cp_expression=sphf
enu_expression=nulf
cond_expression=cndf
(property rho1 is :dnsf:!:condtn:)
(property enul is :nulf:!:condtn:)
(property cp1 is :sphf:!:condtn:)
(property prndtl(tem1) is -:cndf:!:condtn:)
(property prndtl(tem1) is -:cndt:!if(iy.gt.:nrgy1:))
(stored var cond is :cndf:!:condtn:) ! These 2 lines affect
(stored var cond is :cndt:!if(iy.gt.:nrgy1:)) ! print-out only
dvo1dt=texf
else ! use properties according to fluid name
fluid_name=:flname: ! set name for use in next-loaded library case
#$089
! where the following character strings are
! declared and set
! rho_expression
! cp_expression
! enu_expression
! emu_expression
! cond_expression
if(:flname:.eq.saturated_water) then
dvo1dt=1.18e-4
endif
if(:flname:.eq.SAE_5W-30_engine_oil) then
dvo1dt=7.0e-4
endif
if(:flname:.eq.Ethylene_Glycol) then
dvo1dt=6.5e-4
endif
if(:flname:.eq.Glycerin) then
mesg(Glycerin
endif
endif ! end of if(.not.use_specified properties)
************************************************************
Group 9.1 make reference properties
(make refdns is 0.0)
(store1 refdns is rho1[,1,]!IF(iz.eq.1))
(make refkvis is 0.0)
(store1 refkvis is enul[,1,]!IF(iz.eq.1))
(make refsph is 0.0)
(store1 refsph is cp1[,1,]!IF(iz.eq.1))
!
(make refcnd is 0.0)
(store1 refcnd is cond[,1,]!IF(iz.eq.1))
************************************************************
GROUP 11. Initialization of variable or porosity fields
fiinit(tem1)=tfin
fiinit(p1)=pout
patch(inisoli,inival,1,nx,nrgy1+1,ny,1,nz,1,1)
coval(inisoli,prps,0.0,103) ! copper
coval(inisoli,rho1,0.0,dnst) ! wall density
coval(inisoli,tem1,0.0,tfin) ! initial wall temperature
coval(inisoli,cond,0.0,cndt) ! wall conductivity
mesg(rho_expression =:rho_expression:
mesg(enu_expression =:enu_expression:
mesg(cond_expression =:cond_expression:
if(:winset:.eq.from_Re) then
wfin=:reno*nulf/diam:
endif
mesg(wfin =:wfin:
(initial of w1 is :wfin:)
(initial of rho1 is :rho_expression:!:condtn:)
(initial of enul is :enu_expression:!:condtn:)
(initial of cond is :cond_expression:!:condtn:)
************************************************************
GROUP 13. Boundary conditions and special sources
************************************************************
Group 13.1 all flow configurations wall-related settings. Needed?
EGWF = T
WALLCO = GRND2
WALL(innNer,north,1,nx,nrgy1,nrgy1,1,nz,1,lstep) !****
coval(inNner,w1,grnd2,0.0) ! which produce inival
coval(inNner,tem1,grnd2,0.0) ! patch
Note: The above lines were introduce for early, not-longer-used
formulation. They have been left as a reminder that their
unexpected production of an inival patch should be
investigated when time permits.
************************************************************
Group 13.2 Gravity momentum source
mesg(group 13.2
gravty
if(gravty) then ! Gravity source
patch(intube,phasem,1,nx,1,nrgy1,1,nz,1,1)
dummy=nx
dummy=xang/dummy
if(gravax.gt.0.0) then
formula=:gravax*dvo1dt:*(tem1-:tenv1:)*sin((2*ix+1)*:dummy/2:)
(source of u1 at intube is :formula:)
formula=:gravax*dvo1dt:*(tem1-:tenv1:)*cos(ix*:dummy:)
(source of v1 at intube is :formula:)
endif
if(gravay.gt.0.0) then
endif
if(gravaz.gt.0.0) then
endif
endif ! End of gravity source
endif ! end of if(gravty)
mesg(after gravity
fdf
(print wfin is wfin)
wfin
************************************************************
if(.not.fdf) then
Group 13.3 Elliptic ! Inlet boundary
mesg(at inlet density=:dnsf: velocity=:wfin:
patch(inlet,low,1,nx,1,nrgy1,1,1,1,lstep)
! the next two give almost identical results
(source of w1 at inlet is coval(0,:wfin:)!onlyms) ! res 2
(source of w1 at inlet is coval(onlyms,:wfin:)) ! res 3
! either of the next two will serve, results being almost identical
(source of p1 at inlet is :wfin:*rho1) ! res1
(source of p1 at inlet is coval(fixflu,:wfin:*rho1)) ! res2
coval(inlet,tem1,onlyms,tfin) !
************************************************************
Group 13.4 Elliptic flow ! Heat loss boundary
integer(izf,izl) !
izf=1;izl=nrgz1 !
patch(tubewal1,north,1,nx,ny,ny,izf,izl,1,lstep) !
coval(tubewal1,tem1,envco1,tenv1) !
************************************************************
Group 13.5 Elliptic flow ! Outlet boundary
patch(outlet,high,1,nx,1,ny,nz,nz,1,lstep) !!
coval(outlet,p1,1.e10,pout) !!
coval(outlet,tem1,onlyms,same) !!
coval(outlet,w1,onlyms,0.0) !!
mesg(end of 13.5 fdf=:fdf !
else ! end of if(.not.fdf) ----------------------------
************************************************************
Group 13.6 fully-developed-flow momentum-source settings
mesg(13.6
spmsfl=:wfin:*rho1*:0.5*(diam*0.5)**2*xulast:
mesg(spmsfl is :spmsfl: !
spmsfl
dnsf
(print spmsfl is spmsfl)
(make camsfl) ! calculated tube-fluid mass flow rate !
(store1 camsfl is sum(w1*ahigh*rho1)) !
(print camsfl is camsfl)
(make spmsfl) ! specified fluid mass flow rate !
(store1 spmsfl is camsfl!if(isweep.eq.1)) !
(print spmsfl is spmsfl)
!
(make florat is 1.0) ! calculated flow/specified flow !
(store1 florat is camsfl/spmsfl) !
(print florat is florat)
eop
(make pgrmul is 1.0) ! prgrad multiplier to adjust flow rate !
(store1 pgrmul is pgrmul*(0.01+1/florat)/1.01) !
(store1 pgrmul is pgrmul*(1+1/florat)/2) !
(store1 pgrmul is pgrmul*(90+1/florat)/91) !
(print presgrad_mult is pgrmul) !
!
(make pgrad is 0.0) !
real(volum) !
volum=0.5*0.25*diam**2*xulast*zwlast !
(store1 pgrad is shrz[,:nrgy1:,]*4.0/:diam:) ! first estimate !
(print press_grad is pgrad) !
!
dbs
(make pgrad is 2.0)
(make pgrad is 20.0)
real(volum) !
volum=0.5*0.25*diam**2*xulast*zwlast !
(store1 pgrad is 0.9*shrz[,:nrgy1:,]*4.0/:diam:!if(isweep.eq.2))
(store1 pgrad is (pgrad+1.0-florat)!if(isweep.gt.1))
! This is the most sucessful in procuring convergence so far.
! I am suspicious of it because it is not dimensionless
(print pgrad is pgrad)
(make frccoe is 0.0) !
dummy=2.0/(:wfin:**2) !
(store1 frccoe is shrz[,:nrgy1:,]*:dummy:/refdns) !
(print frdcoe is frccoe)
(store1 fcre is frccoe*reyno) !
!
mesg(whole !
patch(whole,volume,1,nx,1,nrgy1,1,nz,1,1) !
(source of w1 at whole is pgrad*pgrmul) !
(source of w1 at whole is pgrad) !
!
!
************************************************************ !
Group 13.7 fully-developed-flow heat-source settings !
There are several ways of computing the heat flow to the external!
fluid, namely !
1. as the net source of the tubewall source !
patch(tubewal1,north,1,nx,ny,ny,1,1,1,lstep) !
coval(tubewal1,tem1,envco1,tenv1) !
deactivated !
(make heafl1 is 0.0) ! heat flux to environment !
(store1 heafl1 is (-nets(tem1,tubewal1))) !
(print heafl1 is heafl1) !
!
2. as the temperature difference times the ap at iy=ny !
deactivatd !
(stored var aptm is apco(tem1)) ! essential !
(make heafl2 is 0.0) ! heat flux to environment !
(store1 heafl2 is (tem1-:tenv1:)*aptm with if(iy.eq.:ny:)) !
(print heafl2 is heafl2) !
!
3. as the heatflux within the metal !
(stored var antm is anco(tem1)) !
(make heafl3 is 0.0) ! heat flux through wall !
(store1 heafl3 is (tem1-tem1[,+1,])*antm with if(iy.eq.:ny-1:)) !
(print heafl3 is heafl3)
4. at tube inner surface !
this one is used !
(make heafl4 is 0.0) ! heat flux to environment !
(store1 heafl4 is (tem1-tem1[,+1,])*antm with if(iy.eq.:ny-2:)) !
!
patch(fluid,high,1,nx,1,nrgy1,1,nz,1,1) ! tube cross-section patch !
!
(make entflo is 0.0) !
(store1 entflo at fluid is sum((tem1-:tenv1:)*cp1*rho1*w1*ahigh)) !
real(entflin) ! enthalpy flowing in !
entflin=spmsfl*sphf*(tfin-tenv1) !
(make entfli is 0.0) !
(store1 entfli at fluid is sum((:tfin:-:tenv1:)*cp1*rho1*w1*ahigh))!
(make heafac is 1.0) ! should become unity after enough sweeps !
(store1 heafac is entfli/entflo) !
(make mult is 1.0) ! multiplies the dtdz !
(store1 mult is mult*(10.0+heafac)/11.0) !
(store1 mult is mult*(1.0+heafac)/2.0) !
!
real(dtdz) ! This if, strictly speaking, MINUS the temperature !
! gradient !
! first estimate the Nusselt number !
real(reyn,pran,nulam,nutur,nusse) !
reyn=diam*:wfin:*dnsf/nulf ! Reynolds using specified properties !
pran=sphf*dnsf*nulf/cndf !
reyn !
pran !
nulam=4.5 !
nutur=0.023*reyn**0.8*pran**0.3333 !
if(nutur.gt.nulam) then !
nusse=nutur !
else !
nusse=nulam !
endif !
!
! Now follows the source calculation for fully developed flow !
dtdz=(tfin-tenv1)*(nusse*cndf)*diam/(:wfin:*dnsf*sphf) !
mesg(const tgrad=:tgrad: twall=:twall: !
! inital guess not affecting result !
if(tgrad) then !
tgrad !
(source of tem1 at fluid is cp1*rho1*w1*(:dtdz:)*mult) !
endif !
The above is the difference between the enthalpy flow in from !
upstream and the ethalpy flow out to downstream.
************************************************************ !
Group 13.4 fully-developed-flow momentum-source settings !
!
(make effect is 0.0) !
(print heat_flux is heafl4)
(print entfli is entfli)
(store1 effect is heafl4/entfli) !
(print effectiveness is effect) !
if(twall) then !
twall !
(source of tem1 at fluid is cp1*rho1*w1*(tem1-:tenv1:)*$ !
mult*:dtdz:/(tem1[1,1,1]-:tenv1:)) !
endif !
!
************************************************************ !
Group 13.5 fully-developed-flow temperature-source settings !
(make maxtem is 0.0) !
(store1 maxtem is tem1[1,1,1]) !
!
(make maxvel is 0.0) !
(store1 maxvel is w1[1,1,1]) !
(print inlet_vel is :wfin:) !
!
(make ohtc is 0.0) !
(store1 ohtc is heafl4/(anorth[,:nrgy1:,]*(:tfin-tenv1:))) !
************************************************************ !
endif
************************************************************
GROUP 15. Termination of sweeps
LSWEEP = numits ;ISWC1 = 1
if(toend) then
isg21=lsweep
else
isg21=1
endif
numits
isg21
LITHYD = 1 ;LITFLX = 1 ;LITC = 1 ;ITHC1 = 1
SELREF = T
RESFAC =1.0E-03
************************************************************
GROUP 16. Termination of iterations
resref(1)=0.0
selref=f ! best not used now that residual calculations have
! been changed without exploration of consequences
************************************************************
GROUP 17. Under-relaxation devices
OVRRLX =0.
EXPERT = F ;NNORSL = F
conwiz=f ! automatic convergence promotion
spedat(rlxfac,refvel,r,0.1*:wfin:)
spedat(rlxfac,reflen,r,zwlast)
relax(w1,linrlx,relxw1)
relax(tem1,linrlx,relxt1)
if(fdf) then
relax(tem1,linrlx,0.2*relxt1)
endif
relax(enul,linrlx,0.01)
relax(enul,linrlx,1.)
if(:turb:.ne.none) then
relax(enut,linrlx,0.5)
endif
if(:turb:.eq.kemodl) then
relax(ke,linrlx,0.1)
relax(ep,linrlx,0.1)
varmax(ep)= 1.e15
endif
************************************************************
mesg(group 18
GROUP 18. Limits on variables or increments to them
how justify these settings?
varmax(tem1)=tenv1+(tfin-tenv1)*10 !****
varmin(tem1)=tenv1
varmin(w1)=0.0
varmax(w1)=2.5*:wfin:
**********************************************************
GROUP 19. Data communicated by satellite to GROUND
PARSOL = F
************************************************************
GROUP 20. Preliminary print-out
! note that the following is rather clumsy, because the inforout-
! printing coding allows only reals to be printed, A better method
! must be devised
************************************************************
GROUP 21. Print-out of variables
SUBWGR = F
nprint=lsweep/howmany
* Y in OUTPUT argument list denotes:
* 1-field 2-correction-eq. monitor 3-selective dumping
* 4-whole-field residual 5-spot-value table 6-residual table
if(.not.prprps) then
output(prps,n,n,y,n,n,n)
endif
if(:turb:.ne.none) then -------------------------------------------
!
if(.not.prenut) then !
output(enut,n,n,y,n,n,n)
endif
if(.not.prlen1) then !
output(len1,n,n,y,n,n,n) !
endif !
!
if(turmo.eq.2.) then !----------------------------------------- !
! !
if(.not.prltls) then ! !
output(ltls,n,n,y,n,n,n) ! !
endif ! !
! !
if(.not.prwgap) then ! ! !
output(wgap,n,n,y,n,n,n) ! ! !
endif ! ! !
endif ! end eq 2 ------------------------------------------- ! !
endif ! end ifturb.ne. none -------------------------------------
************************************************************
GROUP 22. Spot-value print-out
IXMON = 1 ;IYMON = 1 ;IZMON = 1
NPRMON = 100000 ;NPRMNT = 1 ;TSTSWP = -1
UWATCH = T ;USTEER = T
HIGHLO = F
char(moni)
moni=maxabs
moni=maxmin
#:moni:
************************************************************
GROUP 22.1. Inforin print-out
write(>inforin,SimScene is TubeFlow of March 2013)
write(>>inforin,winset is :winset:)
write(>>inforin,fluid is :flname:)
write(>>inforin,gravty is :gravty:)
write(>>inforin,turbulence model is :turb:)
write(>>inforin,flow_form is :floform:)
write(>>inforin,print_general is :prgen:)
write(>>inforin,print_boun_conds is :prboco:)
write(>>inforin,print_dimensionless is :prdino:)
write(>>inforin,priny_grid&geom is :prgrge:)
write(>>inforin,print_ref_props is :prrepr:)
write(>>inforin,print_numerical is :prnum:)
write(>>inforin,print_all is :prntall:)
write(>>inforin, )
#endpause
(print wsource is nets(w1,inisoli))
real(xsect)
xsect=0.125*xang*diam**2
(make force is 0.0)
(store1 force is (p1[1,1,1]-:pout:)*:xsect:)
(print force is force)
************************************************************
GROUP 22.2. Inforout print-out
if(prgrge.or.prntall) then ! grid- and geometry-related
(print used_nrgy1 is :nrgy1:)
(print used_nrgz1 is :nrgz1:)
(print tube_length is :lng1:)
(print tube_inner_diam is :diam:)
if(.not.fdf) then !
(print length/diameter is :zwlast/diam:) !
endif ! end not fdf !
(print tube_wall_thickness is :thck:)
endif
if(prrepr.or.prntall) then
(print reference_density is refdns)
(print reference_kin_visc is refkvis)
(print reference_spc_heat is refsph) !
(print reference_conductivity is refcnd) !
endif
if(prdino.or.prntall) then ! dimensionless nos. in inforout
(make reyno is 0.0)
dummy=:wfin:
(store1 reyno is min(1.e10,:dummy*diam:/refkvis)) ! min needed to
(print Reynolds_No is reyno) ! prevent overflow
(make prnno is 0.0)
(store1 prnno is refsph*refkvis*refdns/refcnd)
(print Prandtl_No is prnno)
(print effectiveness is effect)
(make nusno is 0.0)
if(fdf) then
(store1 nusno is ohtc*:diam:/refcnd)
(print Nusselt_No is nusno)
(make nuMca is 0.0)
(store1 nuMca is 0.023*reyno^0.8*prnno^0.3333)
(print Nu_from_McAdams is nuMca)
(print 4*friction_coeff is 4*frccoe)
(make stanno is 0.0)
(store1 stanno is ohtc/(refdns*refsph*:wfin:))
(print st*pr**2/3 is stanno*prnno^0.667)
endif
endif
if(prboco.or.prntall) then ! boundary conds.. in inforout
(print inlet_velocity is :wfin:)
(print inlet_temperature is :tfin:)
(print outlet_pressure is :pout:)
endif
if(prntall) then !
(print friction_coeff is frccoe) !
if(fdf) then
(print spec_mas_flo is :spmsfl:) !
(print calc_mas_flo is camsfl) !
(print calc_flow/spec_flow is florat) !
(print presgrad_mult is pgrmul) !
(print press_grad is pgrad) !
endif
(print max_temp_ratio is (maxtem-:tenv1:)/:tfin-tenv1:) !
(print heafl4 is heafl4) !
(print enthalpy_outflow is entflo) ! print-out not required !
(print total_heat_out is entflo+heafl4) ! print-out not required !
(print total_heat_out is entflo) ! print-out not required !
(print spec_enth_inflow is entfli) ! print-out not required !
(print heat_inflow/outflo is heafac) ! should become unity !
(print multiplier is mult) ! print-out not required !
(print max_vel_ratio is maxvel/:wfin:) !
(print 4*friction_coeff is frccoe) !
(print frccoe*reyno is fcre)
(print used_lsweep is :lsweep:)
endif !
************************************************************
GROUP 23. Field print-out and plot control
IDISPA=1;IDISPB=2;IDISPC=NZ
NUMCLS = 5
IPLTF = 1 ;IPLTL = -1 ;NPLT = -1
ISWPRF = 1 ;ISWPRL = 100000
ITABL = 3 ;IPROF = 1
ABSIZ =0.5 ;ORSIZ =0.4
NTZPRF = 1 ;NCOLPF = 50
ICHR = 2 ;NCOLCO = 45 ;NROWCO = 20
yzpr=t
iprof=3
patch(crossstr,profil,1,1,1,ny,nz/2,nz/2,1,1)
coval(crossstr,w1,0.0,0.0)
coval(crossstr,tem1,0.0,0.0)
if(:turb:.ne.none) Then
coval(crossstr,enut,0.0,0.0)
endif
if(:turb:.eq.lvel) then
coval(crossstr,wdis,0.0,0.0)
endif
if(:turb:.eq.kemodl) then
coval(crossstr,ke,0.0,0.0)
coval(crossstr,len1,0.0,0.0)
endif
if(nz.gt.1) then
patch(axial,profil,1,1,1,1,1,nz,1,1)
coval(axial,w1,0.0,0.0)
coval(axial,tem1,0.0,0.0)
patch(outersrf,profil,1,1,ny,ny,2,nz,1,1)
coval(outersrf,#hfl,0.0,0.0)
patch(nearwall,profil,1,1,nrgy1,nrgy1,2,nz,1,1)
if(fdf) then
coval(nearwall,shrz,0.0,0.0)
endif
coval(nearwall,htco,0.0,0.0)
endif
inforout
Use properties at outlet temperature in Reynolds etc numbers
because these change least
mesg(inforout
Note for nip. It is desirable to be able print integers, logicals
and character strings to inforout
(print Euler_No is frccoe*:zwlast/diam:) ! do not print as of no
! interest for tubes
if(.not.fdf) then !
real(deltin) ! inlet temperature difference !
deltin=tfin-tenv1 !
(make mflo is 0.0) !
(store1 mflo is nets(r1,inlet)) !
(print mass_flow_rate is mflo) !
!
(make ohtc is 0.0) !
real(area) !
area=xulast*0.5*diam*zwlast !
(store1 ohtc is hflx/(lmtd*:area:)) !
if(prntall) then !
(print overall_ht_tr/co is ohtc) !
endif !
!
(make velrat is 0.0) !
(store1 velrat is w1[1,1,9*NZ/10]/:wfin:) ! 9/10 needed ? !
if(prntall) then !
(print max/inlet_velocity is velrat) !
endif !
(make deltou is 0.0) !
(store1 deltou is tfou-:tenv1:) !
(print deltou is deltou)
(print deltin is :deltin:)
!
(make lmtd is 0.0) !
(store1 lmtd is (:deltin:-deltou)/loge(:deltin:/deltou)) !
if(prntall) then !
(print log_mean_t_diff is lmtd) !
endif !
(make hflx is 0.0)
(store1 hflx is -nets(tem1,tubewal1))
if(prntall) then !
(print total_heat_flux is hflx) !
!
(print inlet_temp. is :tfin:) !
(print external_temp. is :tenv1:) !
(print outlet_temp. is tfou) !
(print effectiveness is (:tfin:-tfou)/:tfin-tenv1:) !
endif !
endif ! end of not fdf ----------------------------------------- !
if(.not.fdf) then !--------------------------------------
(make frccoe is 0.0) !
dummy=diam/(zwlast*:wfin:**2) !
formula=:dummy:*(p1[1,1,1]-p1[1,1,nz])/refdns !
(store1 frccoe is :formula:) !
(store1 nusno is ohtc*:diam:/refcnd)
(print Nusselt_No is nusno)
(make nuMca is 0.0)
(store1 nuMca is 0.023*reyno^0.8*prnno^0.3333)
(print Nu_from_McAdams is nuMca)
(print 4*friction_coeff is 4*frccoe)
(make stanno is 0.0)
(store1 stanno is ohtc/(refdns*refsph*:wfin:))
(print st*pr**2/3 is stanno*prnno^0.667)
if(prntall) then !
(print friction_coeff is frccoe) !
endif !
(print 4*friction_coeff is 4*frccoe) !
!
(make tfou is 0.0) !
(store1 tfou is :tfin:-hflx/(mflo*cp1[,1,])!IF(IZ.EQ.1)) !
endif ! ???? end of not. fdf) ----------------------------
izprf=24
nz=30
nzprin=1
izdb1=26
debug=t
dbout=t
flag=t
fdf
tstswp=1
PHOTON macro
mesg(gravity = :gravty:
real(vfact)
vfact=zwlast/diam
write(>u,p)
write(>>u, )
if(.not.gravty) then ! 2D only, when gravity is absent
if(fdf) then
write(>>u, )
write(>>u, )
write(>>u,AUTOPLOT )
write(>>u, )
write(>>u, )
write(>>u,file; phi 5 )
write(>>u, )
write(>>u,cl; da 1; w1; colf 1 )
write(>>u, msg W-velocity profile)
write(>>u, pause)
write(>>u, dump; velocrve)
write(>>u,cl; da 1; tem1; colf 1)
write(>>u, msg temperature profile)
write(>>u, pause)
write(>>u, dump; tempcrve)
write(>>u,cl; da 1; enul; colf 1)
write(>>u, msg laminar kinematic viscosity profile)
write(>>u, pause)
write(>>u, dump; enulcrve)
if(:turb:.ne.none) then
write(>>u,cl; da 1; enut; colf 1)
write(>>u, msg turbulent kinematic viscosity profile)
write(>>u, pause)
write(>>u, dump; enutcrve)
endif
write(>>u,msg( enter E then RETURN to quit)
write(>>u,enduse)
endif
write(>>u,:vfact: 1)
write(>>u, )
write(>>u,vi -x)
write(>>u,gr ou x 1 col 1)
write(>>u,window 0.0 0.95 0.4 0.80)
write(>>u,gr ou x 1 col 1)
write(>>u,con tem1 x 1 fi;.001)
write(>>u,pause; dump tem1; con off; red)
if(:flname:.ne.use_specified_properties) then
write(>>u,con enul x 1 y 1 :nrgy1: fi;.001)
write(>>u,pause; dump enul; con off; red)
endif
write(>>u,con w1 x 1 fi ;.001)
write(>>u,pause; dump w1; con off; red)
write(>>u,con p1 x 1 fi ;.001)
write(>>u,pause; dump p1; con off; red)
if(:turb:.ne.none) then
write(>>u,con enut x 1 fi ;.001)
write(>>u,pause; dump enut; con off; red)
endif
write(>>u,msg( enter E then RETURN to quit)
write(>>u,enduse)
else ! 3D when gravity is active
mesg(about to write u
write(>>u,;;;)
write(>>u, )
write(>>u,vi z)
write(>>u,gr z 1)
write(>>u,pause)
integer(izcon)
if(fdf) then
izcon=1
else
izcon=nz/2
endif
write(>>u,con tem1 z :izcon: fi;.0001)
write(>>u,pause; dump tem1; con off; red)
write(>>u,gr of;red;gr ou z 1)
write(>>u,con u1 z :izcon: fi;.0001)
write(>>u,pause; dump u1; con off; red)
write(>>u,con v1 z :izcon: fi;.0001)
write(>>u,pause; dump v1; con off; red)
write(>>u,con w1 z :izcon: fi;.0001)
write(>>u,pause; dump w1; con off; red)
dummy=:wfin:
write(>>u,set ref vec :0.5*dummy:)
write(>>u,vec z :izcon: sh)
write(>>u,pause; dump vec; con off; red)
write(>>u,msg( enter E then RETURN to quit)
write(>>,enduse
endif
************************************************************
GROUP 24. Dumps for restarts
SAVE = T ;NOWIPE = F
NSAVE =CHAM
************************************************************
GROUP 25. VR-related settings
************************************************************
GROUP 26. Graphical-display macros
************************************************************
endit(tem1)=0.0
resref(tem1)=0.0
nyprin=1
iyprf=ny-10
STOP