c c =================================================================== subroutine fmod1(maxmx,mvar,meqn,maux,mbc,mx,q,ql,qr, & aux,dx,dt,method,cfl,fm,fp,q1,aux1, & amdq,apdq) c =================================================================== c c # Copyright (C) 2002 Ralf Deiterding c # Brandenburgische Universitaet Cottbus c c # Copyright (C) 2003-2007 California Institute of Technology c # Ralf Deiterding, ralf@cacr.caltech.edu c implicit double precision (a-h,o-z) include "ck.i" include "cuser.i" c parameter (lenma=46, maxm2 = 10002) !# assumes at most 10000 grid points with mbc=2 dimension q(mvar, 1-mbc:maxmx+mbc) dimension ql(mvar, 1-mbc:maxmx+mbc) dimension qr(mvar, 1-mbc:maxmx+mbc) dimension fm(mvar, 1-mbc:maxmx+mbc) dimension fp(mvar, 1-mbc:maxmx+mbc) dimension aux(maux, 1-mbc:maxmx+mbc) dimension q1(1-mbc:maxmx+mbc, meqn) dimension aux1(1-mbc:maxmx+mbc, maux) dimension amdq(1-mbc:maxmx+mbc, meqn) dimension apdq(1-mbc:maxmx+mbc, meqn) dimension method(7) dimension a(lenma, -1:maxm2) dimension fdiv(leNsp+2), hR(leNsp), hL(leNsp) c logical tdiff data tdiff /.false./ !# activate thermal diffusion c mu = Nsp+1 mE = Nsp+2 mT = Nsp+3 c marho = 1 mau = marho+1 mapr = mau+1 maW = mapr+1 macp = maW+1 mavis = macp+1 macon = mavis+1 maX1 = macon+1 maXN = maX1+Nsp-1 maD1 = maXN+1 maDN = maD1+Nsp-1 maTD1 = maDN+1 maTDN = maTD1+Nsp-1 mawsx = maTDN+1 c if (lenma.lt.mawsx) then write (6,*) 'Set lenma in fnavs1g to ',maTDN,'!' stop endif c c # Compute derived quantities do 10 i = 1-mbc, mx+mbc c # Calculate total density a(marho,i) = 0.d0 do k=1,Nsp a(marho,i) = a(marho,i)+q(k,i) enddo a(mau,i) = q(mu,i)/a(marho,i) c # Calculate mass fractions do k=1,Nsp a(maX1+k-1,i) = q(k,i)/a(marho,i) enddo call ckmmwy(a(maX1,i),iwork,rwork,a(maW,i)) a(macp,i) = avgtabip(q(mT,i),a(maX1,i),cpk,Nsp) c c # Calculate mole fractions do k=1,Nsp a(maX1+k-1,i) = a(maX1+k-1,i)*a(maW,i)/Wk(k) enddo c a(mapr,i) = a(marho,i)/a(maW,i)*RU*q(mT,i) c # Mixture viscosity call mcavis(q(mT,i),a(maX1,i),rmcwork,a(mavis,i)) c # Mixture diffusion coefficients call mcadif(a(mapr,i),q(mT,i),a(maX1,i),rmcwork,a(maD1,i)) c # Mixture thermal conductivity call mcacon(q(mT,i),a(maX1,i),rmcwork,a(macon,i)) c # Thermal diffusion ratios for light species if (tdiff) & call mcatdr(q(mT,i),a(maX1,i),imcwork,rmcwork,a(maTD1,i)) c 10 continue c fac = 2.d0*dt/(dx**2) do 110 i = 1, mx+1 c # Difference quotients ux = (a(mau,i)-a(mau,i-1))/dx Tx = (q(mT,i)-q(mT,i-1))/dx px = (a(mapr,i)-a(mapr,i-1))/dx c c # Upwinding based on velocity u, if both sides have c # the same velocity, otherwise use average iR = i iL = i-1 c if (a(mau,i).ge.0.d0.and.a(mau,i-1).ge.0.d0) iR = i-1 c if (a(mau,i).lt.0.d0.and.a(mau,i-1).lt.0.d0) iL = i c c # Construct intermediate values and material properties rho = 0.5d0*(a(marho,iR)+a(marho,iL)) u = 0.5d0*(a(mau,iR)+a(mau,iL)) T = 0.5d0*(q(mT,iR)+q(mT,iL)) p = 0.5d0*(a(mapr,iR)+a(mapr,iL)) W = 0.5d0*(a(maW,iR)+a(maW,iL)) vis = 0.5d0*(a(mavis,iR)+a(mavis,iL)) con = 0.5d0*(a(macon,iR)+a(macon,iL)) cp = 0.5d0*(a(macp,iR)+a(macp,iL)) c call tabintp(q(mT,iR),hR,hms,Nsp) call tabintp(q(mT,iL),hL,hms,Nsp) c c # Diffusive momentum flux fdiv(mu) = (4.d0/3.d0)*vis*ux c # Thermal conductivity, energy change due to diffusive momentum fdiv(mE) = con*Tx+u*fdiv(mu) c c # Stability condition from diffusive momentum and energy flux cflvis = max(fac*dabs(vis/rho),fac*dabs(con/(rho*cp))) c do k=1, Nsp c # Species-dependent difference quotient Xx = (a(maX1+k-1,i)-a(maX1+k-1,i-1))/dx c # Construct species-dependent intermediate values Xk = 0.5d0*(a(maX1+k-1,iR)+a(maX1+k-1,iL)) rhok = 0.5d0*(q(k,iR)+q(k,iL)) Dmix = 0.5d0*(a(maD1+k-1,iR)+a(maD1+k-1,iL)) c Dmix = (a(maD1+k-1,iR)*a(maD1+k-1,iL))/ c & (0.5d0*(a(maD1+k-1,iR)+a(maD1+k-1,iL))) TDmix = 0.5d0*(a(maTD1+k-1,iR)+a(maTD1+k-1,iL)) hk = 0.5d0*(hR(k)+hL(k)) c c # Diffusion driving force dk = Xx + (Xk-rhok/rho) * px/p c c # Diffusive multi-species flux fdiv(k) = dk c # add thermal diffusion if (tdiff) fdiv(k) = fdiv(k) + TDmix*Tx/T fdiv(k) = fdiv(k) * rho*Wk(k)*Dmix/W c # incorporate multi-species flux into energy flux fdiv(mE) = fdiv(mE) + hk*fdiv(k) c # Dufour effect if (Xk.gt.0.d0.and.tdiff) & fdiv(mE) = fdiv(mE) + p*Dmix*TDmix*dk/Xk cflvis = max(cflvis, fac*dabs(Dmix)) enddo c cfl = max(cfl,cflvis) c c # Add diffusive flux approximation do k = 1, mE fm(k,i) = fm(k,i) - fdiv(k) fp(k,i) = fp(k,i) - fdiv(k) enddo c 110 continue end