[Nek5000-users] 2D Axisymmetric in semi-circular annulus
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Tue Aug 28 02:00:01 CDT 2012
Hi Adrian,
My suspicion is that the issue is somehow related to
the SYM bc + curve sides + ifaxis. Generally, SYM on
curve sides implies using the stress formulation.
At present, I'm not certain of the status of the
combination:
IFSTRS T
IFAXIS T
We can look into this as well.
Sorry that this has been holding you up.
Paul
On Mon, 27 Aug 2012, nek5000-users at lists.mcs.anl.gov wrote:
> 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