wall shear stress
485 views
Skip to first unread message
Usama Niaz
Jul 1, 2021, 3:02:33 AM7/1/21
to Nek5000
Hi everyone.
I am quite new to nek5000, I am working on airfoil for which I require wall shear stress. I have follow the userchk subroutine from turbchannel tutorial by setting the argument to true, highlighted below. But when I visualize it in ParaView it doesn't show any parameter related to wall shear stress. it there anything I am missing? thanks
subroutine
subroutine userchk
include 'SIZE'
include 'TOTAL'
real x0(3)
data x0 /0.0, 0.0, 0.0/
save x0
integer icalld
save icalld
data icalld /0/
real atime,timel
save atime,timel
integer ntdump
save ntdump
real rwk(INTP_NMAX,ldim+1) ! r, s, t, dist2
integer iwk(INTP_NMAX,3) ! code, proc, el
save rwk, iwk
integer nint, intp_h
save nint, intp_h
logical iffpts
save iffpts
real xint(INTP_NMAX),yint(INTP_NMAX),zint(INTP_NMAX)
save xint, yint, zint
save igs_x, igs_z
parameter(nstat=9)
real ravg(lx1*ly1*lz1*lelt,nstat)
real stat(lx1*ly1*lz1*lelt,nstat)
real stat_y(INTP_NMAX*nstat)
save ravg, stat, stat_y
save dragx_avg
logical ifverbose
common /gaaa/ wo1(lx1,ly1,lz1,lelv)
& , wo2(lx1,ly1,lz1,lelv)
& , wo3(lx1,ly1,lz1,lelv)
real tplus
real tmn, tmx
integer bIDs(1)
save iobj_wall
n = nx1*ny1*nz1*nelv
nelx = NUMBER_ELEMENTS_X
nely = NUMBER_ELEMENTS_Y
nelz = NUMBER_ELEMENTS_Z
if (istep.eq.0) then
bIDs(1) = 1
call create_obj(iobj_wall,bIDs,1)
nm = iglsum(nmember(iobj_wall),1)
if(nid.eq.0) write(6,*) 'obj_wall nmem:', nm
call prepost(.true.,' ')
endif
ubar = glsc2(vx,bm1,n)/volvm1
e2 = glsc3(vy,bm1,vy,n)+glsc3(vz,bm1,vz,n)
e2 = e2/volvm1
if (nfield.gt.1) then
tmn = glmin(t,n)
tmx = glmax(t,n)
endif
if(nid.eq.0) write(6,2) time,ubar,e2,tmn,tmx
2 format(1p5e13.4,' monitor')
if (time.lt.tSTATSTART) return
c
c What follows computes some statistics ...
c
if(ifoutfld) then
if (ldimt.ge.2) call lambda2(t(1,1,1,1,2))
if (ldimt.ge.3) call comp_vort3(t(1,1,1,1,3),wo1,wo2,vx,vy,vz)
endif
if(icalld.eq.0) then
if(nid.eq.0) write(6,*) 'Start collecting statistics ...'
nxm = 1 ! mesh is linear
call interp_setup(intp_h,0.0,nxm,nelt)
nint = 0
if (nid.eq.0) then
nint = INTP_NMAX
call cfill(xint,XCINT,size(xint))
do i = 1,INTP_NMAX
yi = (i-1.)/(INTP_NMAX-1)
yint(i) = tanh(BETAM*(2*yi-1))/tanh(BETAM)
enddo
call cfill(zint,ZCINT,size(zint))
endif
iffpts = .true. ! dummy call to find points
call interp_nfld(stat_y,ravg,1,xint,yint,zint,nint,
$ iwk,rwk,INTP_NMAX,iffpts,intp_h)
iffpts = .false.
call gtpp_gs_setup(igs_x,nelx ,nely,nelz,1) ! x-avx
call gtpp_gs_setup(igs_z,nelx*nely,1 ,nelz,3) ! z-avg
call rzero(ravg,size(ravg))
dragx_avg = 0
atime = 0
timel = time
ntdump = int(time/tSTATFREQ)
icalld = 1
endif
dtime = time - timel
atime = atime + dtime
! averaging over time
if (atime.ne.0. .and. dtime.ne.0.) then
beta = dtime / atime
alpha = 1. - beta
ifverbose = .false.
call avg1(ravg(1,1),vx ,alpha,beta,n,'uavg',ifverbose)
call avg2(ravg(1,2),vx ,alpha,beta,n,'urms',ifverbose)
call avg2(ravg(1,3),vy ,alpha,beta,n,'vrms',ifverbose)
call avg2(ravg(1,4),vz ,alpha,beta,n,'wrms',ifverbose)
call avg3(ravg(1,5),vx,vy,alpha,beta,n,'uvmm',ifverbose)
call avg1(ravg(1,6),t ,alpha,beta,n,'tavg',ifverbose)
call avg2(ravg(1,7),t ,alpha,beta,n,'trms',ifverbose)
call avg3(ravg(1,8),vx,t ,alpha,beta,n,'utmm',ifverbose)
call avg3(ravg(1,9),vy,t ,alpha,beta,n,'vtmm',ifverbose)
call torque_calc(1.0,x0,.true.,.true.) ! compute wall shear
dragx_avg = alpha*dragx_avg + beta*dragx(iobj_wall)
endif
timel = time
! write statistics to file
if(istep.gt.0 .and. time.gt.(ntdump+1)*tSTATFREQ) then
! averaging over statistical homogeneous directions (x-z)
do i = 1,nstat
call planar_avg(wo1 ,ravg(1,i),igs_x)
call planar_avg(stat(1,i),wo1 ,igs_z)
enddo
if (nfield.gt.1) then
! evaluate d<T>/dy at the lower wall
call opgrad(wo1,wo2,wo3,stat(1,6))
call dssum(wo2,lx1,ly1,lz1)
call col2(wo2,binvm1,n)
call interp_nfld(stat_y,wo2,1,xint,yint,zint,nint,
$ iwk,rwk,INTP_NMAX,iffpts,intp_h)
dTdy_w = stat_y(1)
else
dTdy_w = 1.
endif
! extract data along wall normal direction (1D profile)
call interp_nfld(stat_y,stat,nstat,xint,yint,zint,nint,
$ iwk,rwk,INTP_NMAX,iffpts,intp_h)
ntdump = ntdump + 1
if (nid.ne.0) goto 998
rho = param(1)
dnu = param(2)
A_w = 2 * XLEN * ZLEN
tw = dragx_avg / A_w
u_tau = sqrt(tw / rho)
qw = -param(8) * dTdy_w
t_tau = 1/u_tau * qw
Re_tau = u_tau / dnu
tplus = time * u_tau**2 / dnu
write(6,*) 'Dumping statistics ...', Re_tau, t_tau
open(unit=56,file='vel_fluc_prof.dat')
write(56,'(A,1pe14.7)') '#time = ', time
write(56,'(A)')
$ '# y y+ uu vv ww uv tt ut -vt'
open(unit=57,file='mean_prof.dat')
write(57,'(A,1pe14.7)') '#time = ', time
write(57,'(A)')
$ '# y y+ Umean Tmean'
do i = 1,nint
yy = 1+yint(i)
write(56,3)
& yy,
& yy*Re_tau,
& (stat_y(1*nint+i)-(stat_y(0*nint+i))**2)/u_tau**2,
& stat_y(2*nint+i)/u_tau**2,
& stat_y(3*nint+i)/u_tau**2,
& stat_y(4*nint+i)/u_tau**2,
& (stat_y(6*nint+i)-(stat_y(5*nint+i))**2)/t_tau**2,
& stat_y(7*nint+i)/qw,
& -stat_y(8*nint+i)/qw
write(57,3)
& yy,
& yy*Re_tau,
& stat_y(0*nint+i)/u_tau,
& stat_y(5*nint+i)/t_tau
3 format(1p15e17.9)
enddo
close(56)
close(57)
998 endif
return
end
Emmanuel Gillyns
Jul 1, 2021, 5:53:49 AM7/1/21
to Nek5000
Hi,
Setting those parameters to True will not output something that can be visualized in paraview.
The parameters of the subroutine torque_calc() are : scaling (use 1.0), the variable to use (x0), then the two variable you set to true are to show verbose mode in the logfile.
If you want to understand how it works in more details, you should open go to the file "Nek5000/core/navier5.f" at lines 1623 to 1821.
Hope this helps you,
Emmanuel
Reply all
Reply to author
Forward
0 new messages