[gchnhn]

Hello everybody,

I have downloaded a 2d simple euler fvm solver from http://code.google.com/p/cybo/. When I try to read the code, I come across some several variables qs1,qs2,fs,alpha,dfs, and can't figure out their meanings. Is alpha a factor of diffusion? The paper wrote by the author is not so clear, and he don't reply my email. I hope you can explain to me in more detailed way.

Here is part of the program that I got confused. The following attachment is a zip package of the code. The relevant paper by the author is too large and you can download it from http://code.google.com/p/cybo/downlo...o.pdf&can=2&q=

SUBROUTINE get_flux(flux)
USE euler
USE mesh
USE inputs, ONLY: gamma,mach
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(4,numpts), INTENT(INOUT) :: flux
DOUBLE PRECISION, DIMENSION(numpts) :: div
DOUBLE PRECISION, DIMENSION(4) :: fs,dfs,norm
INTEGER :: i,n1,n2,t1,t2
INTEGER :: bc,e
DOUBLE PRECISION :: dx,dy,qs1,qs2,alpha,c1,c2,len,u1,u2,v1,v2
!! Routine to get the flux (residual) for each cell volume

!! T1 Diagram of how the edge flux is used
!! /\ to update the total flux of T1 and T2.
!! / \
!! Tri 1 / \
!! / \
!! N1/__Edge__\N2
!! \ /
!! \ /
!! Tri 2 \ /
!! \ /
!! \/
!! T2

! Loop over the interior edges to get the flux balance of
! corresponding nodes
div = 0.0
!CALL get_div(div)
DO i=1,size(inter) !对内边循环
t1 = edg(1,inter(i)) ! Node 1 of tri 1
n1 = edg(2,inter(i)) ! Node 2 of tri 1/ node 1 of edge
t2 = edg(3,inter(i)) ! Node 1 of tri 2
n2 = edg(4,inter(i)) ! Node 2 of tri 2/ node 2 of edge

dx = x(n1) - x(n2) ! Get dx for edge
dy = y(n1) - y(n2) ! Get dy for edge

qs1 = (rhou(n1)*dy - rhov(n1)*dx)/rho(n1) ! Reused in flux
qs2 = (rhou(n2)*dy - rhov(n2)*dx)/rho(n2) ! Reused in flux

fs(1) = .5d0*(qs1*rho(n1) + qs2*rho(n2))
fs(2) = .5d0*(qs1*rhou(n1) + qs2*rhou(n2)) + .5d0*(p(n1) + p(n2))*dy
fs(3) = .5d0*(qs1*rhov(n1) + qs2*rhov(n2)) - .5d0*(p(n1) + p(n2))*dx
fs(4) = .5d0*(qs1*(rhoE(n1)+p(n1)) + qs2*(rhoE(n2)+p(n2)))

! Add scalar diffusion
c1 = sqrt( p(n1)*gamma/rho(n1)) !棱边两端的声速
c2 = sqrt( p(n2)*gamma/rho(n2)) !棱边两端的声速
dx = x(t2) - x(t1) ! Get dx for edge
dy = y(t2) - y(t1) ! Get dy for edge
len = sqrt(dx**2 + dy**2)

alpha = ( abs(qs1 + qs2)/2.0d0/len + (c1 + c2)/2.0d0 ) * len
dfs = - alpha/2.0d0*(w(1:4,n1)-w(1:4,n2))

! Add edge fluxes up for each T point
flux(:,t1) = flux(:,t1) - fs / area(t1)
flux(:,t2) = flux(:,t2) + fs / area(t2)

! Add diffusive fluxes for each N point
flux(:,n1) = flux(:,n1) - dfs / area(n1)
flux(:,n2) = flux(:,n2) + dfs / area(n2)

END DO
! Some loop over the boundary edges
.....

END SUBROUTINE