这是我的代码,你看下电流密度为什么会从3900变到3400
原理是什么,什么公式
! For modeling oxygen transport in Catalyst layer, details can be found in the following literature.
! Mesoscopic modeling impacts of liquid water saturation, and platinum distribution on gas transport resistances in a PEMFC catalyst layer
! Mesoscopic modeling of transport resistances in a polymer-electrolyte fuel-cell catalyst layer: Analysis of hydrogen limiting currents
! Please cite this reference: YT Mu, SR Yang, P He, WQ Tao, Electrochimica Acta 388, 138659
! Please cite this reference: YT Mu, AZ Weber, ZL Gu, WQ Tao, Applied Energy 255, 113895-113908
! Developed by Yutong Mu from Xi'an Jiaotong University, Email: yutongmu@mail.xjtu.edu.cn
module start
integer i,j,k,ij,mm,iter,nprocs,istart,iend,ntot,ncase
parameter (nxmax=100,nprocsx=10,nx=nxmax/nprocsx)
parameter (ny=100,nz=100)
parameter (q = 19) ! D3Q19 model
double precision,parameter::XL=2000.E-9,YL=300.E-9,ZL=300.E-9
integer,dimension(0:nx+1,0:ny+1,0:nz+1)::ls
logical,dimension(0:nx+1,0:ny+1,0:nz+1)::wallg,walle,wallPt
double precision,dimension(0:nx+1,0:ny+1,0:nz+1)::Por_o2,Naf_o2,diff
double precision,dimension(0:nx+1,0:ny+1,0:nz+1)::chmcal,chPt
! D3Q19 discrete velocities
integer,dimension(1:q)::fcx=(/1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 1,-1, 1,-1, 0, 0, 0, 0, 0/)
integer,dimension(1:q)::fcy=(/0, 0, 1,-1, 0, 0, 1, 1,-1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0/)
integer,dimension(1:q)::fcz=(/0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 1, 1,-1,-1, 1, 1,-1,-1, 0/)
double precision,dimension(0:q-1,0:nx+1,0:ny+1,0:nz+1)::fp,ffp,fn,ffn
double precision,dimension(0:q-1)::w ! weight coefficients
double precision::tau_p,tau_n ! BGK relaxation times
double precision::chgin,Dg,Delec0,Delec,Dliq
double precision::GNa_pn,Gup1_pn,Gup2_pn,Nup1_pn,Nup2_pn
double precision::DX,DT,scale,sita,fx,EtaD,j0,kk1,kk2,kk3,PtO,Enerst
double precision::farady,gascons,concref,excoeff,currentref,Etain,pressure,tmeaue
double precision::one,two,react,RPt,RPt0,SumRate,TotRate,lamda_eq,scale_o2
double precision::RGDL,act,henry
double precision::CurDensity,Rcl_delta,Rcl_old
double precision::porosity
double precision::Rcl,xiter0
double precision::RF,SumRF,SumRF0,ECSA
double precision::kgdl
double precision::totPot,totPor
double precision::sumPot,sumPor
double precision::K_GDL,K_CL,miu_liq,distance,theta_CL,theta_GDL
double precision::Pc_CL,Pc_GDL,sigma,por_GDL,sat_GDL
CHARACTER*20::name,fname,name1
integer idx,myid
integer mxp,mxm
integer nallgrp,ierr
end module
PROGRAM MAIN
use start
use mpi
integer :: ramp_duration, iter_load
double precision :: Etain_ramp_step, Etain_initial, Etain_final
integer istatus(MPI_STATUS_SIZE)
integer, parameter :: iter_total = 200000
integer, parameter :: output_freq = 200
integer :: transient_output_unit
double precision :: sumNaf, totNaf
! =================== 1. initialize MPI ===================
call MPI_INIT(ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
nallgrp=MPI_COMM_WORLD
idx=myid
mxp=mod(idx+1,nprocsx)
mxm=mod(idx+nprocsx-1,nprocsx)
istart=1; iend=nx
if(idx.eq.0) istart=3
if(idx.eq.nprocsx-1) iend=nx-2
!====================== initialize structure =====================
dx=3.0d-9
fx=0.5d0
pressure=101325.d0+110000.d0
tmeaue=273.15d0+80.d0
farady=96487.d0
gascons=8.314d0
concref=3.39d0
excoeff=0.5d0
currentref=60.d0
Dg =1.47e-6
! *** MODIFICATION 1: Break the circular dependency ***
! First, set a fixed, stable value for tau_p. It MUST be > 0.5.
tau_p = 0.6d0
! Initialize D3Q19 weights
call init_D3Q19
! Second, calculate dt and other parameters based on the fixed tau_p.
j0=1.d0/7.d0
sita=(1.d0-j0)*(tau_p-0.5d0)/3.d0
scale=Dg/sita
dt=dx**2.d0/scale
EtaD=(1.d0-j0)/3.d0
lamda_eq = 8.7d0
Delec0=(1.14698d0/10000000000.d0)*(lamda_eq**0.708d0)
henry=101325.d0/(4.408d0-0.09712d0*lamda_eq)
kp_pn=henry/gascons/tmeaue
scale_o2=1.d0/5.d0 !scale_o2=10.0d0/5.d0
Delec=scale_o2*Delec0
! Third, calculate tau_n based on the fixed tau_p and other parameters.
tau_n=0.5d0+Delec/Dg*(tau_p-0.5d0)
if(myid.eq.0) write(*,*) 'tau_p (fixed), tau_n (calculated), kp_pn:'
if(myid.eq.0) write(*,*) tau_p,tau_n,kp_pn
GNa_pn=-1.d0/kp_pn-1.d0
Gup1_pn=1.d0/kp_pn-1.d0
Gup2_pn=-2.d0
Nup1_pn=-1.d0/kp_pn+1.d0
Nup2_pn=-2.d0/kp_pn
if(myid.eq.0) write(*,*) GNa_pn,Gup1_pn,Gup2_pn,Nup1_pn,Nup2_pn
! Now that dt is a valid non-zero number, these calculations will be correct.
kk1=currentref*(dt/dx)/concref
kk2=excoeff*farady/gascons/tmeaue
kk3=3000.d0/gascons/tmeaue
Rpt=7.3d0; Rpt0=Rpt*dx/dt
! Ramp loading parameters
Etain_initial = 1.0d0
Etain_final = 0.6d0
iter_load = 20000
ramp_duration = 1
Etain_ramp_step = (Etain_final - Etain_initial) / dble(ramp_duration)
Etain = Etain_initial
Enerst=1.2d0-Etain
PtO=1.d0/(1.d0+exp(22.4d0*(0.818d0-Enerst)))
react=(1.d0-PtO)*kk1*exp(kk2*Etain-kk3*PtO)
one=react*(EtaD*Rpt0+1.d0)/4.d0/farady
two=one+EtaD
if(myid.eq.0) then
write(*,*) "Initial Etain: ", Etain
write(*,*) "Initial PtO, one, two: ", PtO, one, two
transient_output_unit = 20
open(transient_output_unit, file='transient_response.dat')
write(transient_output_unit, '(A)') '# Time(s), CurrentDensity(A/m^2), Rcl(s/m), Avg_O2_Pore, Avg_O2_Nafion'
endif
if(myid.eq.0) write(*,*) "Initial PtO = ", PtO, kk2*Etain, kk3*PtO, one, two
RGDL=20.d0
kgdl=dt/dx/RGDL
if(myid.eq.0) write(*,*) kk1,kk2,kk3,RGDL
call MYTSTRUCTURE
chgtot =0.48952d0
chgin =chgtot
do k=1,NZ
do j=1,NY
do i=1,NX
if(.not.wallg(i,j,k)) then
Por_o2(i,j,k)=chgin
Naf_o2(i,j,k)=chgin/kp_pn
call initialize_distributions(fp(:,i,j,k), Por_o2(i,j,k), w)
call initialize_distributions(ffp(:,i,j,k), Por_o2(i,j,k), w)
call initialize_distributions(fn(:,i,j,k), Naf_o2(i,j,k), w)
call initialize_distributions(ffn(:,i,j,k), Naf_o2(i,j,k), w)
else if(.not.walle(i,j,k)) then
Por_o2(i,j,k)=chgin
Naf_o2(i,j,k)=chgin/kp_pn
call initialize_distributions(fp(:,i,j,k), Por_o2(i,j,k), w)
call initialize_distributions(ffp(:,i,j,k), Por_o2(i,j,k), w)
call initialize_distributions(fn(:,i,j,k), Naf_o2(i,j,k), w)
call initialize_distributions(ffn(:,i,j,k), Naf_o2(i,j,k), w)
else
Por_o2(i,j,k)=chgin
Naf_o2(i,j,k)=chgin/kp_pn
call initialize_distributions(fp(:,i,j,k), Por_o2(i,j,k), w)
call initialize_distributions(ffp(:,i,j,k), Por_o2(i,j,k), w)
call initialize_distributions(fn(:,i,j,k), Naf_o2(i,j,k), w)
call initialize_distributions(ffn(:,i,j,k), Naf_o2(i,j,k), w)
endif
enddo
enddo
enddo
Rcl_old=0.0
do ITER=1,iter_total
if (ITER >= iter_load .and. ITER < iter_load + ramp_duration) then
Etain = Etain + Etain_ramp_step
Enerst=1.2d0-Etain
PtO=1.d0/(1.d0+exp(22.4d0*(0.818d0-Enerst)))
react=(1.d0-PtO)*kk1*exp(kk2*Etain-kk3*PtO)
one=react*(EtaD*Rpt0+1.d0)/4.d0/farady
two=one+EtaD
else if (ITER >= iter_load + ramp_duration) then
if (Etain /= Etain_final) then
Etain = Etain_final
Enerst=1.2d0-Etain
PtO=1.d0/(1.d0+exp(22.4d0*(0.818d0-Enerst)))
react=(1.d0-PtO)*kk1*exp(kk2*Etain-kk3*PtO)
one=react*(EtaD*Rpt0+1.d0)/4.d0/farady
two=one+EtaD
endif
endif
CALL COLLISION_BGK
CALL STREAM_D3Q19
CALL LBM_BOUNDARY_D3Q19
CALL MACRO_BOUNDARY
IF(mod(ITER,output_freq).eq.0) then
CurDensity = TotRate*4.d0*farady/real(ny*nz)
Rcl=chgin*4.d0*farady/CurDensity
Rcl_delta=abs(Rcl-Rcl_old)/(Rcl+1.E-20)
Rcl_old=Rcl
call Cal_RoughFactor
if(myid.eq.0) write(*,*) "RoughFactor= (m Pt 2 * m MEA -2)",SumRF
chgin=0.85d0*chgin+0.15d0*(chgtot*Rcl/(Rcl+RGDL))
totPot=0.d0
totPor=0.d0
totNaf=0.d0
do k=1,nz
do j=1,ny
do i=istart,iend
if((.not.wallg(i,j,k))) then
totPot=totPot+Por_o2(i,j,k)
totPor=totPor+1.d0
else if (.not.walle(i,j,k)) then
totNaf=totNaf+Naf_o2(i,j,k)
endif
enddo
enddo
enddo
call MPI_REDUCE(totPot,sumPot,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr)
call MPI_REDUCE(totPor,sumPor,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr)
call MPI_REDUCE(totNaf,sumNaf,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr)
if(myid.eq.0) then
write(*,*) "iter,CurDensity,Rcl,Rcl*sumRF :"
write(*,*) iter,CurDensity,Rcl,Rcl*sumRF
write(*,*) "chgin,Rcl_delta :"
write(*,*) chgin,Rcl_delta
write(*,*) "henry,tau_n,Delec :"
write(*,*) henry,tau_n,Delec
write(*,*) "sumPot/sumPor, lamda_eq "
write(*,*) sumPot/sumPor,lamda_eq
write(transient_output_unit, '(E12.5, 4(A, E12.5))') &
ITER*dt, ',', CurDensity, ',', Rcl, ',', &
sumPot/(sumPor+1.e-20), ',', sumNaf/((nxmax-4)*ny*nz*porosity + 1.e-20)
endif
if(iter.gt.20000 .and. Rcl_delta.le.1.d-5) then
goto 110
endif
endif
enddo
110 continue
if (myid.eq.0) then
close(transient_output_unit)
endif
call output
call MPI_FINALIZE(ierr)
CONTAINS
!************************************************************************
! *** MODIFICATION 2: This subroutine now only initializes the weights. ***
subroutine init_D3Q19
use start
! D3Q19 weight coefficients
w(0) = 1.0d0/3.0d0
w(1:6) = 1.0d0/18.0d0
w(7:18) = 1.0d0/36.0d0
if(myid.eq.0) then
write(*,*) "D3Q19 BGK model weights initialized"
endif
end subroutine
subroutine initialize_distributions(f, rho, weight)
use start
double precision, intent(out) :: f(0:q-1)
double precision, intent(in) :: rho, weight(0:q-1)
integer :: m
do m = 0, q-1
f(m) = rho * weight(m)
enddo
end subroutine
subroutine equilibrium_D3Q19(feq, rho, u, v, w_vel, weight)
use start
double precision, intent(out) :: feq(0:q-1)
double precision, intent(in) :: rho, u, v, w_vel, weight(0:q-1)
integer :: m
double precision :: cu, c_sq, u_sq
u_sq = u*u + v*v + w_vel*w_vel
do m = 0, q-1
cu = fcx(m+1)*u + fcy(m+1)*v + fcz(m+1)*w_vel
feq(m) = rho * weight(m) * (1.d0 + 3.d0*cu + 4.5d0*cu*cu - 1.5d0*u_sq)
enddo
end subroutine
subroutine MYTSTRUCTURE
use start
use mpi
integer istatus(MPI_STATUS_SIZE)
integer lentot,ii
double precision,dimension(0:nxmax+1,0:ny+1,0:nz+1)::difftot
integer,dimension(0:nxmax+1,0:ny+1,0:nz+1)::lstot
if(myid.eq.0) then
lstot=0
porosity=0
open(10, file="LSTOT.dat")
do k=1,nz
do j=1,ny
do i=1,nxmax
read(10,*) lstot(i,j,k)
if(i.ge.3 .and. i.le.nxmax-2 .and. lstot(i,j,k).eq.0) porosity=porosity+1
enddo
enddo
enddo
close(10)
lstot(0:2,:,:)=0
lstot(nxmax-1:nxmax+1,:,:)=0
lstot(:,0,:)=lstot(:,ny,:)
lstot(:,ny+1,:)=lstot(:,1,:)
lstot(:,:,0)=lstot(:,:,nz)
lstot(:,:,nz+1)=lstot(:,:,1)
porosity=real(porosity)/real(ny*nz*(nxmax-4))
write(*,*) "Porosity = ", porosity
difftot(:,:,:)=Dg
open(10, file="DIFFTOT.dat")
do k=1,nz
do j=1,ny
do i=1,nxmax
read(10,*) difftot(i,j,k)
enddo
enddo
enddo
close(10)
difftot(0:2,:,:)=Dg
difftot(nxmax-1:nxmax+1,:,:)=Dg
endif
call MPI_BCAST (porosity,1,MPI_DOUBLE_PRECISION,0,nallgrp,ierr)
lentot=(nxmax+2)*(ny+2)*(nz+2)
call MPI_BCAST(lstot,lentot,MPI_integer,0,nallgrp,ierr)
do k=0,nz+1
do j=0,ny+1
do i=0,nx+1
ii=nx*idx+i
ls(i,j,k)=lstot(ii,j,k)
enddo
enddo
enddo
lentot=(nxmax+2)*(ny+2)*(nz+2)
call MPI_BCAST(difftot,lentot,MPI_DOUBLE_PRECISION,0,nallgrp,ierr)
do k=0,nz+1
do j=0,ny+1
do i=0,nx+1
ii=nx*idx+i
diff(i,j,k)=difftot(ii,j,k)
enddo
enddo
enddo
wallg=.false.
walle=.false.
wallPt=.false.
do k=0,nz+1
do j=0,ny+1
do i=0,nx+1
if(ls(i,j,k).eq.0) then
diff(i,j,k)=0.5d0+diff(i,j,k)/Dg*(tau_p-0.5d0)
endif
if(ls(i,j,k).ne.0) wallg(i,j,k)=.true.
if(ls(i,j,k).ne.1) walle(i,j,k)=.true.
if(ls(i,j,k).ne.2) wallPt(i,j,k)=.true.
enddo
enddo
enddo
return
end subroutine
Subroutine Cal_RoughFactor
use start
use mpi
integer istatus(MPI_STATUS_SIZE)
double precision,dimension(0:nx+1,0:ny+1,0:nz+1)::ttemp
!===================== Obtain the Real Roughness Factor =========================
ttemp=0.d0
do k=1,nz
do j=1,ny
do i=1,nx
mm=0
if(Naf_o2(i,j,k).gt.1.d-6) then
if(ls(i+1,j,k).eq.2) mm=mm+1
if(ls(i-1,j,k).eq.2) mm=mm+1
if(ls(i,j+1,k).eq.2) mm=mm+1
if(ls(i,j-1,k).eq.2) mm=mm+1
if(ls(i,j,k+1).eq.2) mm=mm+1
if(ls(i,j,k-1).eq.2) mm=mm+1
endif
ttemp(i,j,k)=mm
enddo
enddo
enddo
RF=sum(ttemp(:,:,:))
call MPI_REDUCE(RF,SumRF,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr)
call MPI_BCAST (SumRF,1,MPI_DOUBLE_PRECISION,0,nallgrp,ierr)
if(myid.eq.0) then
open(11,file="UpdateRoughFactor_1.txt")
SumRF=SumRF/real(ny*nz)
write(11,*) "RoughFactor= (m Pt 2 * m MEA -2)",SumRF
ECSA=SumRF/real(nxmax-4)/21.45d0/1.0e6/dx
write(11,*) "Specific surface area over volumePt (m 2 * g -1)=",ECSA
close(11)
endif
return
End Subroutine
subroutine COLLISION_BGK
use start
use mpi
integer istatus(MPI_STATUS_SIZE)
double precision,dimension(0:q-1)::feq_p, feq_n
double precision::rho_p, rho_n
integer :: m_i, j_i, k_i
do k_i=1,nz
do j_i=1,ny
do i=1,nx
if(.not.wallg(i,j_i,k_i)) then
rho_p = 0.d0
do m_i=0,q-1
rho_p = rho_p + fp(m_i,i,j_i,k_i)
enddo
Por_o2(i,j_i,k_i) = rho_p
call equilibrium_D3Q19(feq_p, rho_p, 0.d0, 0.d0, 0.d0, w)
do m_i=0,q-1
ffp(m_i,i,j_i,k_i) = fp(m_i,i,j_i,k_i) - (fp(m_i,i,j_i,k_i) - feq_p(m_i)) / tau_p
enddo
else if(.not.walle(i,j_i,k_i)) then
rho_n = 0.d0
do m_i=0,q-1
rho_n = rho_n + fn(m_i,i,j_i,k_i)
enddo
Naf_o2(i,j_i,k_i) = rho_n
call equilibrium_D3Q19(feq_n, rho_n, 0.d0, 0.d0, 0.d0, w)
do m_i=0,q-1
ffn(m_i,i,j_i,k_i) = fn(m_i,i,j_i,k_i) - (fn(m_i,i,j_i,k_i) - feq_n(m_i)) / tau_n
enddo
endif
enddo
enddo
enddo
return
end subroutine
subroutine STREAM_D3Q19
use start
use mpi
integer istatus(MPI_STATUS_SIZE)
double precision, dimension(0:q-1,0:ny+1,0:nz+1):: tmp1,tmp2,tmp1n,tmp2n
integer left,right,tag1,tag2,tag3,tag4
integer ilen, ii, jj, kk, m_dir
integer :: opp_dir
integer :: i_loc, j_loc, k_loc
! MPI communication for boundary regions
if(nprocsx>1) then
ilen=(ny+2)*(nz+2)*q
right=mxp
left=mxm
if(idx==0) left=MPI_PROC_NULL
if(idx==nprocsx-1) right=MPI_PROC_NULL
tmp1(:,:,:)=ffp(:,nx,:,:)
tmp2=0.d0
tag1=10
call MPI_SENDRECV(tmp1(0,0,0),ilen,MPI_DOUBLE_PRECISION,right&
&,tag1,tmp2(0,0,0),ilen,MPI_DOUBLE_PRECISION,left,tag1,nallgrp,istatus,ierr)
ffp(:,0,:,:)=tmp2(:,:,:)
tmp1(:,:,:)=ffp(:,1,:,:)
tmp2=0.d0
tag2=20
call MPI_SENDRECV(tmp1(0,0,0),ilen,MPI_DOUBLE_PRECISION,left&
&,tag2,tmp2(0,0,0),ilen,MPI_DOUBLE_PRECISION,right,tag2,nallgrp,istatus,ierr)
ffp(:,nx+1,:,:)=tmp2(:,:,:)
tmp1n(:,:,:)=ffn(:,nx,:,:)
tmp2n=0.d0
tag3=30
call MPI_SENDRECV(tmp1n(0,0,0),ilen,MPI_DOUBLE_PRECISION,right&
&,tag3,tmp2n(0,0,0),ilen,MPI_DOUBLE_PRECISION,left,tag3,nallgrp,istatus,ierr)
ffn(:,0,:,:)=tmp2n(:,:,:)
tmp1n(:,:,:)=ffn(:,1,:,:)
tmp2n=0.d0
tag4=40
call MPI_SENDRECV(tmp1n(0,0,0),ilen,MPI_DOUBLE_PRECISION,left&
&,tag4,tmp2n(0,0,0),ilen,MPI_DOUBLE_PRECISION,right,tag4,nallgrp,istatus,ierr)
ffn(:,nx+1,:,:)=tmp2n(:,:,:)
endif
! Periodic boundaries in y and z directions
ffp(:,:,0,:)=ffp(:,:,ny,:)
ffp(:,:,ny+1,:)=ffp(:,:,1,:)
ffp(:,:,:,0)=ffp(:,:,:,nz)
ffp(:,:,:,nz+1)=ffp(:,:,:,1)
ffn(:,:,0,:)=ffn(:,:,ny,:)
ffn(:,:,ny+1,:)=ffn(:,:,1,:)
ffn(:,:,:,0)=ffn(:,:,:,nz)
ffn(:,:,:,nz+1)=ffn(:,:,:,1)
! Perform streaming
do k_loc=1,nz
do j_loc=1,ny
do i_loc=1,nx
do m_dir=0,q-1
ii = i_loc - fcx(m_dir+1)
jj = j_loc - fcy(m_dir+1)
kk = k_loc - fcz(m_dir+1)
! C-style periodic boundary conditions
ii = mod(ii - 1 + nx, nx) + 1
jj = mod(jj - 1 + ny, ny) + 1
kk = mod(kk - 1 + nz, nz) + 1
fp(m_dir,i_loc,j_loc,k_loc) = ffp(m_dir,ii,jj,kk)
fn(m_dir,i_loc,j_loc,k_loc) = ffn(m_dir,ii,jj,kk)
enddo
enddo
enddo
enddo
call Nafion_Pore_D3Q19
call Nafion_Pt_D3Q19
! Bounce-back for solid boundaries
do k_loc=1,nz
do j_loc=1,ny
do i_loc=1,nx
if(wallg(i_loc,j_loc,k_loc)) then
do m_dir=0,q-1
opp_dir = find_opposite(m_dir)
ffp(m_dir,i_loc,j_loc,k_loc) = fp(opp_dir,i_loc,j_loc,k_loc)
enddo
endif
if(walle(i_loc,j_loc,k_loc)) then
do m_dir=0,q-1
opp_dir = find_opposite(m_dir)
ffn(m_dir,i_loc,j_loc,k_loc) = fn(opp_dir,i_loc,j_loc,k_loc)
enddo
endif
enddo
enddo
enddo
return
end subroutine
integer function find_opposite(dir)
integer, intent(in) :: dir
integer :: opp
select case(dir)
case(0)
opp = 0
case(1)
opp = 2
case(2)
opp = 1
case(3)
opp = 4
case(4)
opp = 3
case(5)
opp = 6
case(6)
opp = 5
case(7)
opp = 10
case(8)
opp = 9
case(9)
opp = 8
case(10)
opp = 7
case(11)
opp = 14
case(12)
opp = 13
case(13)
opp = 12
case(14)
opp = 11
case(15)
opp = 18
case(16)
opp = 17
case(17)
opp = 16
case(18)
opp = 15
case default
opp = 0
end select
find_opposite = opp
end function
subroutine LBM_BOUNDARY_D3Q19
use start
use mpi
integer istatus(MPI_STATUS_SIZE)
integer :: m_dir, opp_dir
integer :: j_loc, k_loc
if(idx==0) then
i=1
! Inlet boundary - constant concentration
do m_dir=0,q-1
if(fcx(m_dir+1) > 0) then
opp_dir = find_opposite(m_dir)
do k_loc=1,nz
do j_loc=1,ny
fp(m_dir,i,j_loc,k_loc) = EtaD*chgin - fp(opp_dir,i,j_loc,k_loc)
enddo
enddo
endif
enddo
! Nafion phase at inlet
do m_dir=0,q-1
do k_loc=1,nz
do j_loc=1,ny
fn(m_dir,i,j_loc,k_loc) = ffn(m_dir,i,j_loc,k_loc)
enddo
enddo
enddo
endif
if(idx==nprocsx-1) then
i=nx
! Outlet boundary - zero gradient
do m_dir=0,q-1
if(fcx(m_dir+1) < 0) then
opp_dir = find_opposite(m_dir)
do k_loc=1,nz
do j_loc=1,ny
fp(m_dir,i,j_loc,k_loc) = fp(opp_dir,i,j_loc,k_loc)
enddo
enddo
endif
enddo
do m_dir=0,q-1
do k_loc=1,nz
do j_loc=1,ny
fn(m_dir,i,j_loc,k_loc) = ffn(m_dir,i,j_loc,k_loc)
enddo
enddo
enddo
endif
return
end subroutine
Subroutine output
use start
use mpi
double precision::act0,lamda0
integer istatus(MPI_STATUS_SIZE)
integer :: IU
IU=myid+101
write(name,'(I3.3)') IU
fname='Rate_1_'//TRIM(ADJUSTL(name))//'.DAT'
open (301, file=fname)
do i=1,nx
do j=1,ny
do k=1,nz
write (301,*) 4.d0*farady*chmcal(i,j,k)/dt
enddo
enddo
enddo
close(301)
fname='Cone_1_'//TRIM(ADJUSTL(name))//'.DAT'
open (301, file=fname)
do i=1,nx
do j=1,ny
do k=1,nz
write (301,*) Naf_o2(i,j,k)
enddo
enddo
enddo
close(301)
fname='Conc_1_'//TRIM(ADJUSTL(name))//'.DAT'
open (301, file=fname)
do i=1,nx
do j=1,ny
do k=1,nz
write (301,*) Por_o2(i,j,k)
enddo
enddo
enddo
close(301)
fname='ChPt_1_'//TRIM(ADJUSTL(name))//'.DAT'
open (301, file=fname)
do i=1,nx
do j=1,ny
do k=1,nz
write (301,*) chPt(i,j,k)
enddo
enddo
enddo
close(301)
return
End Subroutine
subroutine MACRO_BOUNDARY
use start
use mpi
integer istatus(MPI_STATUS_SIZE)
integer :: m_dir
do k=1,nz
do j=1,ny
do i=1,nx
if(.not.wallg(i,j,k)) then
Por_o2(i,j,k)=0.d0
do m_dir=0,q-1
Por_o2(i,j,k)=Por_o2(i,j,k)+fp(m_dir,i,j,k)
enddo
Por_o2(i,j,k)=max(Por_o2(i,j,k),0.d0)
Naf_o2(i,j,k)=0.d0
else if(.not.walle(i,j,k)) then
Naf_o2(i,j,k)=0.d0
do m_dir=0,q-1
Naf_o2(i,j,k)=Naf_o2(i,j,k)+fn(m_dir,i,j,k)
enddo
Naf_o2(i,j,k)=max(Naf_o2(i,j,k),0.d0)
Por_o2(i,j,k)=0.d0
else
Por_o2(i,j,k)=0.d0
Naf_o2(i,j,k)=0.d0
endif
enddo
enddo
enddo
return
end subroutine
subroutine Nafion_Pore_D3Q19
use start
use mpi
integer istatus(MPI_STATUS_SIZE)
integer :: m_dir, ii, jj, kk, opp
do k=1,NZ
do j=1,NY
do i=1,NX
if(.not.walle(i,j,k)) then
do m_dir=0,q-1
ii = i + fcx(m_dir+1)
jj = j + fcy(m_dir+1)
kk = k + fcz(m_dir+1)
if(ii >= 1 .and. ii <= nx .and. jj >= 1 .and. jj <= ny .and. kk >= 1 .and. kk <= nz) then
if(.not.wallg(ii,jj,kk)) then
! Interface between Nafion and pore
opp = find_opposite(m_dir)
ffn(m_dir,ii,jj,kk) = Nup1_pn/GNa_pn*ffn(opp,i,j,k) + Nup2_pn/GNa_pn*ffp(m_dir,ii,jj,kk)
endif
endif
enddo
endif
if(.not.wallg(i,j,k)) then
do m_dir=0,q-1
ii = i + fcx(m_dir+1)
jj = j + fcy(m_dir+1)
kk = k + fcz(m_dir+1)
if(ii >= 1 .and. ii <= nx .and. jj >= 1 .and. jj <= ny .and. kk >= 1 .and. kk <= nz) then
if(.not.walle(ii,jj,kk)) then
! Interface between pore and Nafion
opp = find_opposite(m_dir)
ffp(m_dir,ii,jj,kk) = Gup1_pn/GNa_pn*ffp(opp,i,j,k) + Gup2_pn/GNa_pn*ffn(m_dir,ii,jj,kk)
endif
endif
enddo
endif
enddo
enddo
enddo
return
end subroutine
subroutine Nafion_Pt_D3Q19
use start
use mpi
implicit none
integer istatus(MPI_STATUS_SIZE)
double precision::che,cmol
integer::nchPt, opp
integer::ii, jj, kk, m_dir
SumRate=0.d0; chmcal=0.d0; chPt=0.d0
do k=1,nz
do j=1,ny
do i=1,nx
nchPt=0
if(.not.walle(i,j,k)) then
do m_dir=0,q-1
ii = i + fcx(m_dir+1)
jj = j + fcy(m_dir+1)
kk = k + fcz(m_dir+1)
if(ii >= 1 .and. ii <= nx .and. jj >= 1 .and. jj <= ny .and. kk >= 1 .and. kk <= nz) then
if(.not.wallPt(ii,jj,kk)) then
! Interface between Nafion and Pt
opp = find_opposite(m_dir)
che=2.d0*ffn(opp,i,j,k)/two
cmol=(2.d0*ffn(opp,i,j,k)*Rpt0+che)/(EtaD*Rpt0+1.d0)
ffn(m_dir,ii,jj,kk)=EtaD*cmol-ffn(opp,i,j,k)
SumRate=SumRate+(ffn(opp,i,j,k)-ffn(m_dir,ii,jj,kk))
chmcal(i,j,k)=chmcal(i,j,k)+(ffn(opp,i,j,k)-ffn(m_dir,ii,jj,kk))
chPt(i,j,k)=chPt(i,j,k)+che
nchPt=nchPt+1
endif
endif
enddo
endif
chPt(i,j,k)=chPt(i,j,k)/(float(nchPt)+1.E-20)
enddo
enddo
enddo
call MPI_REDUCE(SumRate,TotRate,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr)
call MPI_BCAST (TotRate,1,MPI_DOUBLE_PRECISION,0,nallgrp,ierr)
TotRate=TotRate*dx/dt
return
end subroutine
end PROGRAM MAIN
最新发布