|
[Sponsors] |
October 3, 2010, 04:41 |
Can anyone explain about this function
|
#1 |
New Member
cxcxcx0505
Join Date: May 2010
Posts: 17
Rep Power: 15 |
Do anyone know about this metric function?Where can I find more information about this?And,I have another question,what is the fmag and fangle for in this function?
Thank you. Subroutine Metric Use young Integer :: index(4) Real :: dum1,dum2,dum3,dum4 Real,Dimension(mx,my) :: xxi,xeta,yxi,yeta,xtau,ytau Real :: fudh,fvdh,fanglet Real :: fxtemp1,fytemp1,fxtemp2,fytemp2 Real :: xvel,yvel,fcfind index(1)=ios-1; index(2)=ioe; index(3)=jos-1; index(4)=joe Do j=jos,joe-1 Do i=ios-1,ioe xtau(i,j)=(nxc(i,j)-xc(i,j))/dt ytau(i,j)=(nyc(i,j)-yc(i,j))/dt End Do End Do Call Compact2D_xi(xc,xxi,mx,my,index) Call Compact2D_xi(yc,yxi,mx,my,index) Call Compact2D_eta(xc,xeta,mx,my,index) Call Compact2D_eta(yc,yeta,mx,my,index) Do j=jos-1,joe Do i=ios-1,ioe vol(i,j)=xxi(i,j)*yeta(i,j)-xeta(i,j)*yxi(i,j) xit(i,j)=(-xtau(i,j)*yeta(i,j)+xeta(i,j)*ytau(i,j))/vol(i,j) etat(i,j)=(xtau(i,j)*yxi(i,j)-xxi(i,j)*ytau(i,j))/vol(i,j) xix(i,j)=yeta(i,j)/vol(i,j) xiy(i,j)=-xeta(i,j)/vol(i,j) etax(i,j)=-yxi(i,j)/vol(i,j) etay(i,j)=xxi(i,j)/vol(i,j) q11(i,j)=(xix(i,j)**2+xiy(i,j)**2)*vol(i,j) q22(i,j)=(etax(i,j)**2+etay(i,j)**2)*vol(i,j) End Do End Do Do j=js,je Do i=is,ie x1=xc(i-1,j-1); x2=xc(i,j-1); x3=xc(i,j); x4=xc(i-1,j) y1=yc(i-1,j-1); y2=yc(i,j-1); y3=yc(i,j); y4=yc(i-1,j) xx=x(i,j); yy=y(i,j) dx1=xx-x1; dx2=xx-x2; dx3=xx-x3; dx4=xx-x4 dy1=yy-y1; dy2=yy-y2; dy3=yy-y3; dy4=yy-y4 xi1=xix(i-1,j-1)*dx1+xiy(i-1,j-1)*dy1 eta1=etax(i-1,j-1)*dx1+etay(i-1,j-1)*dy1 xi2=1.+xix(i,j-1)*dx2+xiy(i,j-1)*dy2 eta2=etax(i,j-1)*dx2+etay(i,j-1)*dy2 xi3=1.+xix(i,j)*dx3+xiy(i,j)*dy3 eta3=1.+etax(i,j)*dx3+etay(i,j)*dy3 xi4=xix(i-1,j)*dx4+xiy(i-1,j)*dy4 eta4=1.+etax(i-1,j)*dx4+etay(i-1,j)*dy4 xic=0.25*(xi1+xi2+xi3+xi4) etac=0.25*(eta1+eta2+eta3+eta4) rr(i,j,1)=(1.-etac)*(1.-xic) rr(i,j,2)=(1.-etac)*xic rr(i,j,3)=etac*(1.-xic) rr(i,j,4)=etac*xic End Do End Do ! Compute x- and y- velocities of each point fudh=(fmag2-fmag1)/dt fvdh=0. fanglet=(fangle2-fangle1)/dt Do j=jos-1,joe Do i=ios-1,ioe fxtemp1=xc(i,j)*cos(stroke)+yc(i,j)*sin(stroke)-fmag1 fytemp1=-xc(i,j)*sin(stroke)+yc(i,j)*cos(stroke) fxtemp2=fxtemp1*cos(fangle1)+fytemp1*sin(fangle1) fytemp2=-fxtemp1*sin(fangle1)+fytemp1*cos(fangle1) fcfind=fxtemp2**2/(aa**2)+fytemp2**2/(bb**2) If (fcfind <= 1.) Then xvel=fudh-fxtemp1*fanglet*sin(fangle2-fangle1)-fytemp1*fanglet*cos(fangle2-fangle1) yvel=fvdh+fxtemp1*fanglet*cos(fangle2-fangle1)-fytemp1*fanglet*sin(fangle2-fangle1) uob(i,j)=xvel*cos(stroke)-yvel*sin(stroke) vob(i,j)=xvel*sin(stroke)+yvel*cos(stroke) ck(i,j)=ckzero Else uob(i,j)=0. vob(i,j)=0. ck(i,j)=ckinf End If End Do End Do cmapxi=0. cmapeta=0. Do j=jos,joe-1 Do i=ios,ioe-1 xgrad=abs(ck(i+1,j)-2.*ck(i,j)+ck(i-1,j)) ygrad=abs(ck(i,j+1)-2.*ck(i,j)+ck(i,j-1)) If (xgrad > 1.) Then cmapxi(i+1,j)=1.; cmapxi(i,j)=1.; cmapxi(i-1,j)=1. End If If (ygrad > 1.) Then cmapeta(i,j+1)=1.; cmapeta(i,j)=1.; cmapeta(i,j-1)=1. End If End Do End Do Return End Subroutine Metric |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Compile problem | ivanyao | OpenFOAM Running, Solving & CFD | 1 | October 12, 2012 09:31 |
[blockMesh] BlockMesh FOAM warning | gaottino | OpenFOAM Meshing & Mesh Conversion | 7 | July 19, 2010 14:11 |
latest OpenFOAM-1.6.x from git failed to compile | phsieh2005 | OpenFOAM Bugs | 25 | February 9, 2010 04:37 |
Error with Wmake | skabilan | OpenFOAM Installation | 3 | July 28, 2009 00:35 |
[blockMesh] Axisymmetrical mesh | Rasmus Gjesing (Gjesing) | OpenFOAM Meshing & Mesh Conversion | 10 | April 2, 2007 14:00 |