CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Coeff 11.f90 - calculate the coefficients

Coeff 11.f90 - calculate the coefficients

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
Line 63: Line 63:
Diff_e = Area_e * Gam_e / Del_e
Diff_e = Area_e * Gam_e / Del_e
 +
  Diff_w = Area_w * Gam_w / Del_w  
  Diff_w = Area_w * Gam_w / Del_w  
 +
Diff_s = Area_s * Gam_s / Del_s
Diff_s = Area_s * Gam_s / Del_s
 +
Diff_n = Area_n * Gam_n / Del_n  
Diff_n = Area_n * Gam_n / Del_n  
 +
                 Flux_e = Area_e * Con_e * Ro_e
                 Flux_e = Area_e * Con_e * Ro_e
 +
Flux_w = Area_w * Con_w * Ro_w
Flux_w = Area_w * Con_w * Ro_w
 +
Flux_s = Area_s * Con_s * Ro_s
Flux_s = Area_s * Con_s * Ro_s
 +
 +
Flux_n = Area_n * Con_n * Ro_n
Flux_n = Area_n * Con_n * Ro_n
 +
Aw(i,j) = Diff_w + max(    Con_w,0.)
Aw(i,j) = Diff_w + max(    Con_w,0.)
 +
Ae(i,j) = Diff_e + max(-1.* Con_e,0.)
Ae(i,j) = Diff_e + max(-1.* Con_e,0.)
 +
As(i,j) = Diff_s + max(    Con_s,0.)
As(i,j) = Diff_s + max(    Con_s,0.)
 +
An(i,j) = Diff_n + max(-1.* Con_n,0.)
An(i,j) = Diff_n + max(-1.* Con_n,0.)
 +
         CheckFlux(i,j) =  Flux_e - Flux_w + Flux_s - Flux_n
         CheckFlux(i,j) =  Flux_e - Flux_w + Flux_s - Flux_n
 +
Ap(i,j)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) + CheckFlux(i,j)
Ap(i,j)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) + CheckFlux(i,j)
Line 121: Line 135:
!--------------------------------QUICK SCHEME-------------------------
!--------------------------------QUICK SCHEME-------------------------
-
!------------------------------------------------------------------------
 
-
 
-
!------------------------------------------------------------------------
 
-
 
-
!------------------------------------------------------------------------
 
!-------------------------------- HLPA SCHEME----------------------------
!-------------------------------- HLPA SCHEME----------------------------

Revision as of 01:57, 19 September 2005

Subroutine Coef_1(nf)

include 'icomm_1.f90'

Dimension F_out(nx,ny),CheckFlux(nx,ny)

Character Filename*10


! calculation of fluxes

! all geometry has rectangular 2D notation


Do 100 I= 2,NXmax

Do 100 J= 2,NYmax


Gam_e = ( Gam(i+1,j ) + Gam(i ,j ) ) * 0.5

Gam_w = ( Gam(i-1,j ) + Gam(i ,j ) ) * 0.5

Gam_s = ( Gam(i ,j-1) + Gam(i ,j ) ) * 0.5

Gam_n = ( Gam(i ,j+1) + Gam(i ,j ) ) * 0.5


Ro_e = ( Ro(i+1,j ) + Ro(i ,j ) ) * 0.5

Ro_w = ( Ro(i-1,j ) + Ro(i ,j ) ) * 0.5

Ro_s = ( Ro(i ,j-1) + Ro(i ,j ) ) * 0.5

Ro_n = ( Ro(i ,j+1) + Ro(i ,j ) ) * 0.5

Area_w = Y_xi(i-1,j-1)

Area_e = Y_xi(i ,j-1)

Area_s = X_et(i-1,j-1)

Area_n = X_et(i-1,j )

Del_e = Del_X_et(i ,j )

Del_w = Del_X_et(i-1,j )

Del_n = Del_Y_xi(i ,j )

Del_s = Del_Y_xi(i ,j-1)

! upwind differencing (all other will be included into the source term)

Con_e = Area_e * ( F(i,j,1) + F(i+1,j ,1) ) * 0.5 Con_w = Area_w * ( F(i,j,1) + F(i-1,j ,1) ) * 0.5

Con_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5

Con_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5

Diff_e = Area_e * Gam_e / Del_e


		Diff_w = Area_w * Gam_w / Del_w 


Diff_s = Area_s * Gam_s / Del_s


Diff_n = Area_n * Gam_n / Del_n


               Flux_e = Area_e * Con_e * Ro_e


Flux_w = Area_w * Con_w * Ro_w


Flux_s = Area_s * Con_s * Ro_s


Flux_n = Area_n * Con_n * Ro_n


Aw(i,j) = Diff_w + max( Con_w,0.)


Ae(i,j) = Diff_e + max(-1.* Con_e,0.)


As(i,j) = Diff_s + max( Con_s,0.)


An(i,j) = Diff_n + max(-1.* Con_n,0.)


       CheckFlux(i,j) =  Flux_e - Flux_w + Flux_s - Flux_n


Ap(i,j)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) + CheckFlux(i,j)

Sp(i,j)= 0.

!--------------------------------QUICK SCHEME-------------------------

  go to 700
if( (i.GT.2).AND.(i.LT.NXmax).and.(j.GT.4).AND.(j.LT.NYmax) ) then 
 if(Con_e.GT.0.)  Sp(i,j) = Sp(i,j) + Con_e * (-1.) * (-0.125 * F(i-1,j  ,nf) - 0.25 * F(i  ,j  ,nf) + 0.375 * F(i-1,j  ,nf) )
 if(Con_w.GT.0.)  Sp(i,j) = Sp(i,j) + Con_w * (-1.) * (-0.125 * F(i-2,j  ,nf) - 0.25 * F(i-1,j  ,nf) + 0.375 * F(i-2,j  ,nf) )
 if(Con_n.GT.0.)  Sp(i,j) = Sp(i,j) + Con_n * (-1.) * (-0.125 * F(i  ,j-1,nf) - 0.25 * F(i  ,j  ,nf) + 0.375 * F(i  ,j-1,nf) )
 if(Con_s.GT.0.)  Sp(i,j) = Sp(i,j) + Con_s * (-1.) * (-0.125 * F(i  ,j-2,nf) - 0.25 * F(i  ,j-1,nf) + 0.375 * F(i  ,j-2,nf) )
 if(Con_e.LT.0.)  Sp(i,j) = Sp(i,j) + Con_e * (-1.) * ( 0.375 * F(i  ,j  ,nf) - 0.25 * F(i+1,j  ,nf) - 0.125 * F(i+2,j  ,nf) )
 if(Con_w.LT.0.)  Sp(i,j) = Sp(i,j) + Con_w * (-1.) * ( 0.375 * F(i-1,j  ,nf) - 0.25 * F(i  ,j  ,nf) - 0.125 * F(i+1,j  ,nf) )
 if(Con_n.LT.0.)  Sp(i,j) = Sp(i,j) + Con_n * (-1.) * ( 0.375 * F(i  ,j  ,nf) - 0.25 * F(i  ,j+1,nf) - 0.125 * F(i  ,j+2,nf) )

! if(Con_s.LT.0.) Sp(i,j) = Sp(i,j) + Con_s * (-1.) * ( 0.375 * F(i ,j-1,nf) - 0.25 * F(i ,j+0,nf) - 0.125 * F(i ,j+1,nf) )


end if	
  700 continue	

!--------------------------------QUICK SCHEME-------------------------


!-------------------------------- HLPA SCHEME---------------------------- ! go to 600

! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f)

if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then 
  !------------------ w face -------------------

Fww = F(i-2,j,nf)

Fw = F(i-1,j,nf)

Fp = F(i ,j,nf)

Fe = F(i+1,j,nf)

call HLPA(Con_w,Fww,Fw,Fp,Fe,Delta_f)

Sp(i,j) = Sp(i,j) + Con_w * Delta_f

 !------------------ e face--------------------

Fww = F(i-1,j,nf)

Fw = F(i ,j,nf)

Fp = F(i+1,j,nf)

Fe = F(i+2,j,nf)

call HLPA(Con_e,Fww,Fw,Fp,Fe,Delta_f)

   Sp(i,j) = Sp(i,j) + Con_e * Delta_f * (-1.)
 !------------------ s face--------------------

Fww = F(i ,j-2,nf)

Fw = F(i ,j-1,nf)

Fp = F(i ,j ,nf)

Fe = F(i ,j+1,nf)

call HLPA(Con_s,Fww,Fw,Fp,Fe,Delta_f)


  Sp(i,j) = Sp(i,j) + Con_s * Delta_f
 !------------------ n face--------------------

Fww = F(i ,j-1,nf)

Fw = F(i ,j ,nf)

Fp = F(i ,j+1,nf)

Fe = F(i ,j+2,nf)

call HLPA(Con_n,Fww,Fw,Fp,Fe,Delta_f)


  Sp(i,j) = Sp(i,j) + Con_n * Delta_f *(-1.)


end if	
  600 continue	

!-------------------------------- HLPA SCHEME----------------------------

!------------------------------------------------------------------------

100 continue

!----------------------------------------------------------------

NImax = NXmaxp

NJmax = NYmaxp

F_out = CheckFlux

	Filename  ='0_Flx.txt' 
   Call  Out_array(F_out,NImax,NJmax,Filename)

!-------------------------------------------------------------------


Return

End

My wiki