[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 13:40:49 CDT 2012


P.S. Sorry Adrian -- I misse your point about the spherical shell -- disregard my previous email -- I'll look further.

Best.
Aleks

----- Original Message -----
From: nek5000-users at lists.mcs.anl.gov
To: nek5000-users at lists.mcs.anl.gov
Sent: Tuesday, August 28, 2012 1:35:24 PM
Subject: Re: [Nek5000-users] 2D Axisymmetric in semi-circular annulus

Hi Adrian,

For your geometry, you do not need to use axisymmetric formulation of Nek -- use IFAXI=F, IFAZIV=F and standard userf & userq or take nek5_svn/examples/taylor/taylor.* as a starting point.  

Note though that you do not need to transform a Cartesian geometry into polar one as userdat2() routine does in taylor.usr since you have already generated one in prenek.

Let me know if you have any questions.
Best.
Aleks



----- Original Message -----
From: nek5000-users at lists.mcs.anl.gov
To: nek5000-users at lists.mcs.anl.gov
Sent: Tuesday, August 28, 2012 10:56:46 AM
Subject: Re: [Nek5000-users] 2D Axisymmetric in semi-circular annulus

Hi Adrian,

Sorry for missing a subtle point of your email -- it looks to me that 
vortex2 example is not a relevant example for your case.

Judging by your .rea file, your geometry is a constant z slice, polar 
geometry [r,theta] that does not touch the axis of rotation while vortex2 
example has a constant theata slice [r,z]...

Let me look further and get back to you.
Aleks

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

> Hi Aleks & Paul,
>
> Thanks both for your replies, and for looking into this.
>
> In my problem, I've got a semi-circular annulus that represents the (z,R)-plane for R.ge.0. Two of the elements have one side each that lies entirely along the axis of symmetry y=0 (R=0) -- these are elements 2 & 4, which have side 1 along y=0, in the positive & negative x (z) half-plane, respectively. Is an "A  " BC not appropriate for this case? I thought that an "A  " BC was required for each element that touched the axis of symmetry and that it corresponded with zero radial velocity and zero tangential stress. Have I misunderstood the "A  " boundary condition? ("W  " will apply a tangential stress at the boundary, which is unwanted in my case, so perhaps I should use "SYM" if "A  " is inappropriate).
>
> Thanks,
> Adrian
>
> 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<mailto:nek5000-users at lists.mcs.anl.gov>
> To: nek5000-users at lists.mcs.anl.gov<mailto: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
_______________________________________________
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