[Nek5000-users] 2D Axisymmetric in semi-circular annulus
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Mon Aug 27 22:05:19 CDT 2012
Hi Adrian,
A quick note that BC "A " is used only in a case when domain touches the axis of rotation like in a case of cylinder. But in a case of cylindrical annulus, you should use 'W ' for axial Vx & radial Vr components of velocity (or 'v ' with appropriate values in userbc when ifield=1). In a case of non-trivial azimuthal velocity or swirl (IFAZIV=T), azimuthal velocity component is placed into t(*,*,*,*,1) so one should use 't ' with appropriate values for TEMP in userbc when ifield=2
If in addition to azimuthal velocity (IFAZIV=T), you are also solving for temperature then you temperature array is placed into an array of the first passive scalar, i.e. t(*,*,*,*,2). Then you can set BC & sources using PS1 for it when ifield=3
Let me know if you have questions.
Aleks
----- Original Message -----
From: nek5000-users at lists.mcs.anl.gov
To: nek5000-users at lists.mcs.anl.gov
Sent: Monday, August 27, 2012 5:40:54 PM
Subject: [Nek5000-users] 2D Axisymmetric in semi-circular annulus
Dear Nek users,
I have been attempting to set up a 2D axisymmetric problem with an azimuthal velocity, in a semi-circular annulus container (the 2D axisymmetric version of a spherical shell), but have encountered problems with the terms evolving the azimuthal velocity blowing up near to the axis, and wondered if anybody would be able to help with this?
I have looked at the vortex2 example, and have copied the relevant aspects to my .usr and .rea files. I have set IFAXIS=T to get the correct viscous operators for cylindrical geometry where (x,y)->(z,R), as well as IFHEAT=T, IFAZIV=T, using temperature as the azimuthal velocity. My test mesh is created using prenek, and contains 4 elements on a polar grid from r=0.5 to 1 and theta=0 to 180. Since elements 2 & 4 have one side touching the R=0 (y=0) axis, I have used "A" boundary conditions, which are applied to side 1 for each of these elements, as is required. The remaining BCs are 'SYM' and 'I' for all of the exterior sides. The .rea file containing the mesh with 4 elements is attached to this email. (In the SIZE file I have set lx1=20,lxd=30,lx2=lx1-2).
For a simple test, if I set the initial conditions to be a uniform rotation about the z-axis (setting temp = y), within several time-steps the code blows up, giving NaN's for each quantity (I have tested various initial conditions, which each give the same error). The line of code that appears to cause problems is the line ~240 of hmholtz.f,
TERM1 = BM1(I,J,1,IEL)*U(I,J,1,IEL)/YM1(I,J,1,IEL)**2,
since, when this term is set to 0, the code does not produce NaN's. (Though the viscous operator for the azimuthal velocity is then incorrect...)
When IFHEAT=IFAZIV=F, i.e. when I am not integrating the temperature/azimuthal velocity, the code runs without problem on the same mesh.
I am wondering if there is a problem with how I have set up the mesh that causes this problem? And if so, is there a way to fix this?
Thanks,
Adrian
Modifications to .usr file are:
c-----------------------------------------------------------------------
subroutine userf (ix,iy,iz,eg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
integer e,f,eg
if(y.gt.0) ffy = temp*temp/y
return
end
c-----------------------------------------------------------------------
subroutine userq (ix,iy,iz,eg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
integer e,f,eg
if(y.gt.0) then
visc = param(2)
qvol = -uy*temp/y
endif
return
end
c-----------------------------------------------------------------------
subroutine useric (ix,iy,iz,ieg)
include 'SIZE'
include 'TOTAL'
include 'NEKUSE'
temp = y
return
end
c-----------------------------------------------------------------------
My .rea file:
****** PARAMETERS *****
2.610000 NEKTON VERSION
2 DIMENSIONAL RUN
118 PARAMETERS FOLLOW
1.00000 p001 DENSITY
0.10000E-06 p002 VISCOS
0.00000 p003
1.00000 p004
0.00000 p005
0.00000 p006
1.00000 p007 RHOCP
0.100000E-06 p008 CONDUCT
0.00000 p009
0.00000 p010 FINTIME
0.100000E+08 p011 NSTEPS
0.100000E-03 p012 DT
0.00000 p013 IOCOMM
0.00000 p014 IOTIME
1000.00 p015 IOSTEP
0.00000 p016 PSSOLVER: 0=default
1.00000 p017
0.500000E-01 p018 GRID < 0 --> # cells on screen
-1.00000 p019 INTYPE
5.00000 p020 NORDER
0.100000E-06 p021 DIVERGENCE
0.100000E-08 p022 HELMHOLTZ
0.00000 p023 NPSCAL
0.100000E-01 p024 TOLREL
0.100000E-01 p025 TOLABS
0.450000 p026 COURANT/NTAU
3.00000 p027 TORDER
1.00000 p028 TORDER: mesh velocity (0: p28=p27)
0.00000 p029 = magnetic visc if > 0, = -1/Rm if < 0
0.00000 p030 > 0 ==> properties set in uservp()
0.00000 p031 NPERT: #perturbation modes
0.00000 p032 #BCs in re2 file, if > 0
0.00000 p033
0.00000 p034
0.00000 p035
0.00000 p036
0.00000 p037
0.00000 p038
0.00000 p039
0.00000 p040
0.00000 p041 1-->multiplicative SEMG
0.00000 p042 0=gmres/1=pcg
0.00000 p043 0=semg/1=schwarz
0.00000 p044 0=E-based/1=A-based prec.
0.00000 p045 Relaxation factor for DTFS
0.00000 p046 reserved
0.00000 p047 vnu: mesh matieral prop.
0.00000 p048
0.00000 p049
0.00000 p050
0.00000 p051
0.00000 p052 IOHIS
0.00000 p053
0.00000 p054 fixed flow rate dir: |p54|=1,2,3=x,y,z
0.00000 p055 vol.flow rate (p54>0) or Ubar (p54<0)
0.00000 p056
0.00000 p057
0.00000 p058
0.00000 p059 !=0 --> full Jac. eval. for each el.
0.00000 p060 !=0 --> init. velocity to small nonzero
0.00000 p061
0.00000 p062 >0 --> force byte_swap for output
0.00000 p063 =8 --> force 8-byte output
0.00000 p064 =1 --> perturbation restart
0.00000 p065 #iofiles (eg, 0 or 64); <0 --> sep. dirs
0.00000 p066 output : <0=ascii, else binary
0.00000 p067 restart: <0=ascii, else binary
0.00000 p068 iastep: freq for avg_all (0=iostep)
0.00000 p069
0.00000 p070
0.00000 p071
0.00000 p072
0.00000 p073
0.00000 p074 verbose Helmholtz
0.00000 p075
0.00000 p076
0.00000 p077
0.00000 p078
0.00000 p079
0.00000 p080
0.00000 p081
0.00000 p082
0.00000 p083
0.00000 p084 !=0 --> sets initial timestep if p12>0
0.00000 p085 dt ratio if p84 !=0, for timesteps>0
0.00000 p086 reserved
0.00000 p087
0.00000 p088
0.00000 p089
0.00000 p090
0.00000 p091
0.00000 p092
20.0000 p093 Number of previous pressure solns saved
5.00000 p094 start projecting velocity after p94 step
5.00000 p095 start projecting pressure after p95 step
0.00000 p096
0.00000 p097
0.00000 p098
3.00000 p099 dealiasing: <0--> off/3--> old/4--> new
0.00000 p100
0.00000 p101 Number of additional modes to filter
1.00000 p102 Dump out divergence at each time step
0.100000E-01 p103 weight of stabilizing filter (.01)
0.00000 p104
0.00000 p105
0.00000 p106
0.00000 p107 !=0 --> add to h2 array in hlmhotz eqn
0.00000 p108
0.00000 p109
0.00000 p110
0.00000 p111
0.00000 p112
0.00000 p113
0.00000 p114
0.00000 p115
0.00000 p116 !=0: x elements for fast tensor product
0.00000 p117 !=0: y elements for fast tensor product
0.00000 p118 !=0: z elements for fast tensor product
4 Lines of passive scalar data follows2 CONDUCT; 2RHOCP
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000
12 LOGICAL SWITCHES FOLLOW
T IFFLOW
T IFHEAT
T IFTRAN
T T F F F F F F F F F IFNAV & IFADVC (convection in P.S. fields)
F F T T T T T T T T T T IFTMSH (IF mesh for this field is T mesh)
T IFAXIS
T IFAZIV
T IFSTRS
F IFMGRID
F IFMODEL
F IFKEPS
F IFCHAR
10.0000 10.0000 -4.00000 -5.50000 XFAC,YFAC,XZERO,YZERO
**MESH DATA** 1st line is X of corner 1,2,3,4. 2nd line is Y.
4 2 4 NEL,NDIM,NELV
ELEMENT 1 [ 1C] GROUP 0
-0.4371138E-07 -0.2185569E-07 0.3535533 0.7071066
1.000000 0.5000000 0.3535533 0.7071066
ELEMENT 2 [ 1D] GROUP 0
0.5000000 1.000000 0.7071066 0.3535533
0.000000 0.000000 0.7071066 0.3535533
ELEMENT 3 [ 1D] GROUP 0
0.5962440E-08 0.1192488E-07 -0.7071066 -0.3535533
0.5000000 1.000000 0.7071066 0.3535533
ELEMENT 4 [ 1 ] GROUP 0
-1.000000 -0.5000000 -0.3535533 -0.7071066
0.8742278E-07 0.4371139E-07 0.3535533 0.7071066
***** CURVED SIDE DATA *****
8 Curved sides follow IEDGE,IEL,CURVE(I),I=1,5, CCURVE
2 1 -0.500000 0.00000 0.00000 0.00000 0.00000 C
4 1 1.00000 0.00000 0.00000 0.00000 0.00000 C
2 2 1.00000 0.00000 0.00000 0.00000 0.00000 C
4 2 -0.500000 0.00000 0.00000 0.00000 0.00000 C
2 3 1.00000 0.00000 0.00000 0.00000 0.00000 C
4 3 -0.500000 0.00000 0.00000 0.00000 0.00000 C
2 4 -0.500000 0.00000 0.00000 0.00000 0.00000 C
4 4 1.00000 0.00000 0.00000 0.00000 0.00000 C
***** BOUNDARY CONDITIONS *****
***** FLUID BOUNDARY CONDITIONS *****
E 1 1 3.00000 1.00000 0.00000 0.00000 0.00000
SYM 1 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 1 3 2.00000 3.00000 0.00000 0.00000 0.00000
SYM 1 4 0.00000 0.00000 0.00000 0.00000 0.00000
A 2 1 0.00000 0.00000 0.00000 0.00000 0.00000
SYM 2 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 2 3 1.00000 3.00000 0.00000 0.00000 0.00000
SYM 2 4 0.00000 0.00000 0.00000 0.00000 0.00000
E 3 1 1.00000 1.00000 0.00000 0.00000 0.00000
SYM 3 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 3 3 4.00000 3.00000 0.00000 0.00000 0.00000
SYM 3 4 0.00000 0.00000 0.00000 0.00000 0.00000
A 4 1 0.00000 0.00000 0.00000 0.00000 0.00000
SYM 4 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 4 3 3.00000 3.00000 0.00000 0.00000 0.00000
SYM 4 4 0.00000 0.00000 0.00000 0.00000 0.00000
***** THERMAL BOUNDARY CONDITIONS *****
E 1 1 3.00000 1.00000 0.00000 0.00000 0.00000
I 1 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 1 3 2.00000 3.00000 0.00000 0.00000 0.00000
I 1 4 0.00000 0.00000 0.00000 0.00000 0.00000
A 2 1 0.00000 0.00000 0.00000 0.00000 0.00000
I 2 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 2 3 1.00000 3.00000 0.00000 0.00000 0.00000
I 2 4 0.00000 0.00000 0.00000 0.00000 0.00000
E 3 1 1.00000 1.00000 0.00000 0.00000 0.00000
I 3 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 3 3 4.00000 3.00000 0.00000 0.00000 0.00000
I 3 4 0.00000 0.00000 0.00000 0.00000 0.00000
A 4 1 0.00000 0.00000 0.00000 0.00000 0.00000
I 4 2 0.00000 0.00000 0.00000 0.00000 0.00000
E 4 3 3.00000 3.00000 0.00000 0.00000 0.00000
I 4 4 0.00000 0.00000 0.00000 0.00000 0.00000
0 PRESOLVE/RESTART OPTIONS *****
7 INITIAL CONDITIONS *****
C Default
C Default
C Default
C Default
C Default
C Default
C Default
***** DRIVE FORCE DATA ***** BODY FORCE, FLOW, Q
4 Lines of Drive force data follow
C
C
C
C
***** Variable Property Data ***** Overrrides Parameter data.
1 Lines follow.
0 PACKETS OF DATA FOLLOW
***** HISTORY AND INTEGRAL DATA *****
0 POINTS. Hcode, I,J,H,IEL
***** OUTPUT FIELD SPECIFICATION *****
6 SPECIFICATIONS FOLLOW
T COORDINATES
T VELOCITY
T PRESSURE
T TEMPERATURE
F TEMPERATURE GRADIENT
0 PASSIVE SCALARS
***** OBJECT SPECIFICATION *****
0 Surface Objects
0 Volume Objects
0 Edge Objects
0 Point Objects
_______________________________________________
Nek5000-users mailing list
Nek5000-users at lists.mcs.anl.gov
https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
More information about the Nek5000-users
mailing list