<div dir="ltr">Dear Paul,<br><br>Thank you for the reply. I tried to implement the code in the usrchk subroutine but found that the values of dgtq and dg are zeros. I tried to identify the wall boundary conditions and found the corresponding the elements , faces and used the code snippet. I am sure I am missing something here.. Any thoughts ?<br>
<br>However the way you said with the postx works. It did give me a profile.out and I am currently working on a way to post-process it. I have close to 25k elements and probably a few hundred elements near the blade and I am not quite sure how to grab the coordinates to enter into the <a href="http://prof.in">prof.in</a> file. Ideally I would like to get the yplus/z plus for the first couple of elements or so. I tried the same way of grabbing the elements and faces near the area of interest, but I am clueless as to how I access the co-ordinate data of those elements. <br>
<br>
I had another question regarding span-wise averaging of shear stress,
Once I calculate the shear stress using the above routine, could I use
the z_average() routine in navier5.f to do the span-wise averaging ? I was looking at the turbchannel.usr where the wall shear stress and averaging process is done but did not make a huge progress..<br><br><code><br>
subroutine usrchk<br> include 'SIZE'<br> include 'TOTAL'<br> include 'ZPER'<br><br> integer e,f,pf<br> real x0(3)<br> character*3 cb_surf<br> real n1,n2,n3<br>
real dg(3,2)<br> real dgtq(3,4)<br> real xm0 (lx1,ly1,lz1,lelt)<br> real ym0 (lx1,ly1,lz1,lelt)<br> real zm0 (lx1,ly1,lz1,lelt)<br> real sij (lx1,ly1,lz1,3*ldim-3,lelv)<br> real pm1 (lx1,ly1,lz1,lelv)<br>
real visc(lx1,ly1,lz1,lelv)<br><br> n=nx1*ny1*nz1*nelv<br> <br> call rzero(dgtq,12)<br> <br> nface = 2 * ndim<br> cb_surf = 'W '<br> do e = 1,nelv<br> flag = 1<br>
do f = 1,nface<br> if(cbc(f,e,1).eq.cb_surf.and.x.gt.-25E-02 .and. x.lt.125E-02)then <br> flag = 0 ! Wall Element <br> pf = eface1(f) ! convert from preproc. notation<br>
<br> js1 = skpdat(1,pf)<br> jf1 = skpdat(2,pf)<br> jskip1 = skpdat(3,pf)<br> js2 = skpdat(4,pf)<br> jf2 = skpdat(5,pf)<br> jskip2 = skpdat(6,pf)<br> <br> if (if3d.or.ifaxis) then<br>
i = 0<br> a = 0 <br> do j2=js2,jf2,jskip2 <br> do j1=js1,jf1,jskip1<br> i = i+1<br> n1 = unx(i,1,f,e)*area(i,1,f,e)<br> n2 = uny(i,1,f,e)*area(i,1,f,e)<br> n3 = unz(i,1,f,e)*area(i,1,f,e)<br>
a = a + area(i,1,f,e)<br>c<br> v = visc(j1,j2,1,e)<br>c<br> s11 = sij(j1,j2,1,1,e)<br> s21 = sij(j1,j2,1,4,e)<br> s31 = sij(j1,j2,1,6,e)<br>c<br> s12 = sij(j1,j2,1,4,e)<br>
s22 = sij(j1,j2,1,2,e)<br> s32 = sij(j1,j2,1,5,e)<br>c<br> s13 = sij(j1,j2,1,6,e)<br> s23 = sij(j1,j2,1,5,e)<br> s33 = sij(j1,j2,1,3,e)<br>c<br> dg(1,1) = pm1(j1,j2,1,e)*n1 ! pressure drag<br>
dg(2,1) = pm1(j1,j2,1,e)*n2<br> dg(3,1) = pm1(j1,j2,1,e)*n3<br>c<br> dg(1,2) = -v*(s11*n1 + s12*n2 + s13*n3) ! viscous drag<br> dg(2,2) = -v*(s21*n1 + s22*n2 + s23*n3)<br> dg(3,2) = -v*(s31*n1 + s32*n2 + s33*n3)<br>
c<br> r1 = xm0(j1,j2,1,e)<br> r2 = ym0(j1,j2,1,e)<br> r3 = zm0(j1,j2,1,e)<br>c<br> do l=1,2<br> do k=1,3<br> dgtq(k,l) = dgtq(k,l) + dg(k,l)<br> enddo<br> enddo<br>
c<br> dgtq(1,3) = dgtq(1,3) + (r2*dg(3,1)-r3*dg(2,1)) ! pressure<br> dgtq(2,3) = dgtq(2,3) + (r3*dg(1,1)-r1*dg(3,1)) ! torque<br> dgtq(3,3) = dgtq(3,3) + (r1*dg(2,1)-r2*dg(1,1))<br>c<br> dgtq(1,4) = dgtq(1,4) + (r2*dg(3,2)-r3*dg(2,2)) ! viscous<br>
dgtq(2,4) = dgtq(2,4) + (r3*dg(1,2)-r1*dg(3,2)) ! torque<br> dgtq(3,4) = dgtq(3,4) + (r1*dg(2,2)-r2*dg(1,2))<br> enddo<br> enddo<br> endif<br> endif<br> enddo<br> enddo<br>
<br> return<br> end<br><br></code><br><br><br>Thanks and Regards<br>Shriram<br><br></div>