[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:19:00 CDT 2012


OK --- that does not appear to be the issue... will dig a bit
deeper.

Paul


On Tue, 28 Aug 2012, nek5000-users at lists.mcs.anl.gov wrote:

>
> 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
>> 
> _______________________________________________
> 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