model.f90

Go to the documentation of this file.

module model
    use constants
    use globals
    use foaming_globals_m
    use fmodena
    use modenastuff
    implicit none
    private
    public molar_balance,kinmodel,odesystem,dim_var
contains
!********************************BEGINNING*************************************
subroutine  odesystem (neq, t, y, ydot)
    use phys_prop, only:rderiv
    integer, intent(in) :: neq
    real(dp), intent(in) :: t
    real(dp), intent(in) :: y(neq)
    real(dp), intent(out) :: ydot(neq)
    integer :: i,j
    real(dp) :: z,zw,ze,zww,zee,lamw,lame,cw,ce,cww,cee,c,dcw,dce,dil,bll
    call dim_var(t,y)
    call molar_balance(y)
    ydot=0
    if (kin_model==4) then
        call kinmodel(y)
        do i=1,size(kineq)
            ydot(kineq(i))=kinsource(i)
        enddo
        if (w0>1e-8_dp) then
            ydot(xweq)=-ydot(kineq(5))/w0*1.0e3_dp
        else
            ydot(xweq)=0
        endif
        ydot(xoheq)=-(ydot(kineq(2))+ydot(kineq(3)))/oh0*1.0e3_dp
    else
        ydot(xoheq) = aoh*exp(-eoh/rg/temp)*(1-y(xoheq))*&
            (nco0-2*w0*y(xweq)-oh0*y(xoheq)) !polyol conversion
        if (kin_model==3) then
            ! gelling influence on kinetics
            if (y(xoheq)>0.5_dp .and. y(xoheq)<0.87_dp) then
                ydot(xoheq)=ydot(xoheq)*&
                    (-2.027_dp*y(xoheq)+2.013_dp)
            endif
            if (y(xoheq)>0.87_dp) then
                ydot(xoheq)=ydot(xoheq)*(3.461_dp*y(xoheq)-2.761_dp)
            endif
        endif
        if (w0>1e-3) then
            ! water conversion
            ! ydot(xWeq) = AW*exp(-EW/Rg/temp)*(1-y(xWeq))*&
            !     (NCO0-2*W0*y(xWeq)-OH0*y(xOHeq)) 2nd order
            ydot(xweq) = aw*exp(-ew/rg/temp)*(1-y(xweq)) !1st order
        endif
        if (dilution) then
            if (co2_pos==1) then
                bll=mb(2)/vsh*mbl(2)/rhop
            else
                bll=mb(1)/vsh*mbl(1)/rhop
            endif
            dil=1/(1+rhop/rhobl*bll)
            ydot(xoheq)=ydot(xoheq)*dil
            ydot(xweq)=ydot(xweq)*dil
        endif
    endif
    !temperature (enthalpy balance)
    ydot(teq) = -dhoh*oh0/(rhop*cp)*ydot(xoheq)-dhw*w0/(rhop*cp)*ydot(xweq)
    do i=1,ngas
        if (geometry=="3D") then
            ydot(teq) = ydot(teq) - dhv(i)*12*pi*mbl(i)*d(i)*radius**4/&
                (rhop*cp*vsh)*(y(fceq+i-1)-kh(i)*y(fpeq+i-1))/(dz(1)/2)
        elseif (geometry=="2D") then
            ydot(teq) = ydot(teq) - dhv(i)*4*pi*mbl(i)*d(i)*radius**2/&
                (rhop*cp*vsh)*(y(fceq+i-1)-kh(i)*y(fpeq+i-1))/(dz(1)/2)
        endif
    enddo
    if (kin_model==4) then
        ydot(kineq(19))=ydot(teq)
    endif
    if (firstrun) then
        if (inertial_term) then
            if (geometry=="3D") then
                ydot(req) = y(req+1)    !radius (momentum balance)
                ydot(req+1) = (sum(y(fpeq:lpeq)) + pair - pamb - &
                    2*sigma/radius - 4*eta*y(req+1)/radius - &
                    3._dp/2*y(req+1)**2)/(radius*rhop)
            elseif (geometry=="2D") then
                ydot(req) = y(req+1)    !radius (momentum balance)
                ydot(req+1) = (sum(y(fpeq:lpeq)) + pair - pamb - &
                    sigma/radius - 2*eta*y(req+1)/radius - &
                    3._dp/2*y(req+1)**2)/(radius*rhop)
            endif
        else
            if (geometry=="3D") then
                ydot(req) = (sum(y(fpeq:lpeq)) + pair - pamb - &
                    2*sigma/radius)*radius/(4*eta)   !radius (momentum balance)
            elseif (geometry=="2D") then
                ydot(req) = (sum(y(fpeq:lpeq)) + pair - pamb - &
                    sigma/radius)*radius/(2*eta)   !radius (momentum balance)
            endif
        endif
        !partial pressure (molar balance)
        do i=fpeq,lpeq
            if (geometry=="3D") then
                ydot(i) = -3*y(i)*ydot(req)/radius + y(i)/temp*ydot(teq) + &
                    9*rg*temp*d(i-fpeq+1)*radius*&
                    (y(fceq+i-fpeq)-kh(i-fpeq+1)*y(i))/(dz(1)/2)
            elseif (geometry=="2D") then
                ydot(i) = -2*y(i)*ydot(req)/radius + y(i)/temp*ydot(teq) + &
                    4*rg*temp*d(i-fpeq+1)*&
                    (y(fceq+i-fpeq)-kh(i-fpeq+1)*y(i))/(dz(1)/2)
            endif
        enddo
    else
        !partial pressure (molar balance)
        do i=fpeq,lpeq
            if (geometry=="3D") then
                ydot(i) = -3*y(i)*rderiv(time)/radius + y(i)/temp*ydot(teq) + &
                    9*rg*temp*d(i-fpeq+1)*radius*&
                    (y(fceq+i-fpeq)-kh(i-fpeq+1)*y(i))/(dz(1)/2)
            elseif (geometry=="2D") then
                ydot(i) = -2*y(i)*rderiv(time)/radius + y(i)/temp*ydot(teq) + &
                    4*rg*temp*d(i-fpeq+1)*&
                    (y(fceq+i-fpeq)-kh(i-fpeq+1)*y(i))/(dz(1)/2)
            endif
        enddo
    endif
    do j=1,ngas
        do i=1,p+1
            if (i==1) then !bubble boundary
                zw=0e0_dp
                z=dz(i)/2
                ze=dz(i)
                zee=ze+dz(i+1)/2
                lame=(ze-z)/(zee-z)
                c=y(fceq+(i-1)*ngas+j-1)
                cee=y(fceq+i*ngas+j-1)
                cw=kh(j)*y(fpeq+j-1)
                ce=cee*lame+c*(1-lame)
                dcw=(c-cw)/(z-zw)
                dce=(cee-c)/(zee-z)
            elseif(i==p+1) then !outer boundary
                zww=z
                zw=ze
                z=zee
                ze=ze+dz(i)
                lamw=(zw-zww)/(z-zww)
                cww=y(fceq+(i-2)*ngas+j-1)
                c=y(fceq+(i-1)*ngas+j-1)
                cw=c*lamw+cww*(1-lamw)
                ce=c
                dcw=(c-cww)/(z-zww)
                dce=0e0_dp
            else
                zww=z
                zw=ze
                z=zee
                ze=ze+dz(i)
                zee=ze+dz(i+1)/2
                lamw=(zw-zww)/(z-zww)
                lame=(ze-z)/(zee-z)
                cww=y(fceq+(i-2)*ngas+j-1)
                c=y(fceq+(i-1)*ngas+j-1)
                cee=y(fceq+i*ngas+j-1)
                cw=c*lamw+cww*(1-lamw)
                ce=cee*lame+c*(1-lame)
                dcw=(c-cww)/(z-zww)
                dce=(cee-c)/(zee-z)
            endif
            !concentration (molar balance)
            if (geometry=="3D") then
                ydot(fceq+(i-1)*ngas+j-1) = &
                    9*d(j)*((ze+radius**3)**(4._dp/3)*dce -&
                    (zw+radius**3)**(4._dp/3)*dcw)/dz(i)
                !reaction source
                if (j==co2_pos) ydot(fceq+(i-1)*ngas+j-1) = &
                    ydot(fceq+(i-1)*ngas+j-1) + w0*ydot(xweq)
            elseif (geometry=="2D") then
                ydot(fceq+(i-1)*ngas+j-1) = &
                    4*d(j)*((ze+radius**2)*dce -&
                    (zw+radius**2)*dcw)/dz(i)
                !reaction source
                !WARNING: reaction source exaggerated to provide reasonable
                !growth
                if (j==co2_pos) ydot(fceq+(i-1)*ngas+j-1) = &
                    ydot(fceq+(i-1)*ngas+j-1) + w0*ydot(xweq)*1e4
            endif
        enddo
    enddo
end subroutine odesystem
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine molar_balance(y)
    real(dp), intent(in) :: y(:)
    integer :: i,j
    !numerical integration
    mb=0e0_dp
    !rectangle rule
    do i=1,p+1
        do j=1,ngas
            mb(j)=mb(j)+y(fceq+j-1+(i-1)*ngas)*dz(i)
        enddo
    enddo
    if (geometry=="3D") then
        mb=mb*4*pi/3 !moles in polymer
    elseif (geometry=="2D") then
        mb=mb*pi !moles in polymer
    endif
    do i=1,ngas
        if (geometry=="3D") then
            mb2(i)=pressure(i)*radius**3*4*pi/(3*rg*temp) !moles in bubble
        elseif (geometry=="2D") then
            mb2(i)=pressure(i)*radius**2*pi/(rg*temp) !moles in bubble
        endif
    enddo
    mb3=mb+mb2 !total moles
end subroutine molar_balance
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine dim_var(t,y)
    use phys_prop, only:physical_properties,rb
    real(dp), intent(in) :: t
    real(dp), intent(in) :: y(:)
    integer :: i
    time=t
    if (firstrun) then
        radius=y(req) ! calculate bubble radius
    else
        radius=rb(time) ! use calculated bubble radius
    endif
    temp=y(teq)
    conv=y(xoheq)
    call physical_properties(temp,conv,radius)
    eqconc=y(fpeq)*kh(1)  !only first gas
    do i=1,ngas
        pressure(i)=y(fpeq+i-1)
    enddo
    do i=1,ngas
        wblpol(i)=mb(i)*mbl(i)/(rhop*vsh)
    enddo
    avconc=mb/vsh
    if (geometry=="3D") then
        laplace_pres=2*sigma/radius
        pair=pair0*r0**3/radius**3*temp/temp0
        porosity=radius**3/(radius**3+s0**3-r0**3)
        st=(s0**3+radius**3-r0**3)**(1._dp/3)-radius !thickness of the shell
    elseif (geometry=="2D") then
        laplace_pres=sigma/radius
        pair=pair0*r0**2/radius**2*temp/temp0
        porosity=radius**2/(radius**2+s0**2-r0**2)
        st=(s0**2+radius**2-r0**2)**(1._dp/2)-radius !thickness of the shell
    endif
    rhofoam=(1-porosity)*rhop
    bub_pres=sum(pressure)+pair-pamb
end subroutine dim_var
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine kinmodel(y)
    real(dp), intent(in) :: y(:)
    integer :: i
    if (kin_model==4) then
        do i=1,size(kineq)
            call modena_inputs_set(kininputs, kininputspos(i+1), y(kineq(i)))
        enddo
    endif
    call modena_inputs_set(kininputs, kininputspos(1), time)
    call modena_inputs_set(kininputs, kininputspos(20), y(teq)-273.15_dp)
    ret = modena_model_call(kinmodena, kininputs, kinoutputs)
    if(ret /= 0) then
        call exit(ret)
    endif
    if (kin_model==4) then
        kinsource=0
        do i=1,size(kineq)
            kinsource(i) = modena_outputs_get(kinoutputs, kinoutputspos(i))
        enddo
    endif
end subroutine kinmodel
!***********************************END****************************************
end module model