Register Blogs Members List Search Today's Posts Mark Forums Read

October 3, 2010, 04:41
#1
New Member

cxcxcx0505
Join Date: May 2010
Posts: 17
Rep Power: 8
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
cmapxi(i+1,j)=1.; cmapxi(i,j)=1.; cmapxi(i-1,j)=1.
End If
cmapeta(i,j+1)=1.; cmapeta(i,j)=1.; cmapeta(i,j-1)=1.
End If
End Do
End Do

Return
End Subroutine Metric
Attached Images
 untitled.jpg (100.9 KB, 6 views)

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post ivanyao OpenFOAM Running, Solving & CFD 1 October 12, 2012 09:31 gaottino OpenFOAM Native Meshers: blockMesh 7 July 19, 2010 14:11 phsieh2005 OpenFOAM Bugs 25 February 9, 2010 05:37 skabilan OpenFOAM Installation 3 July 28, 2009 00:35 Rasmus Gjesing (Gjesing) OpenFOAM Native Meshers: blockMesh 10 April 2, 2007 14:00

All times are GMT -4. The time now is 04:02.