Fick's_law_of_diffusion.f90

!  Here's an example of a Fortran program that solves a challenging mass-transfer problem involving diffusion in a one-dimensional domain under transient conditions. The problem we will model is the diffusion of a substance in a one-dimensional solid, where the concentration of the substance changes over time due to diffusion.

! We'll use the finite difference method to discretize the governing partial differential equation (PDE), which is Fick's law of diffusion:


program mass_transfer_diffusion
    implicit none

    ! Declare variables
    integer, parameter :: nx = 101, nt = 1000   ! Number of spatial and time points
    integer :: i, j
    real(8), parameter :: L = 1.0d0, T = 0.1d0  ! Length of the domain and total time
    real(8) :: D = 1.0d-5d0      ! Diffusion coefficient
    real(8) :: dx, dt            ! Spatial and time steps
    real(8) :: C(nx), C_new(nx)  ! Concentration array
    real(8) :: x, time
    real(8) :: alpha              ! Stability criterion

    ! Initialize space and time step sizes
    dx = L / (nx - 1)
    dt = T / (nt - 1)
    
    ! Stability criterion for the explicit method
    alpha = D * dt / (dx**2)

    ! Check for stability
    if (alpha > 0.5d0) then
        print *, "Warning: The solution may be unstable due to large alpha."
    end if

    ! Initialize concentration profile
    C = 0.0d0  ! Start with all concentrations zero
    C(1) = 1.0d0  ! Boundary condition: concentration at x = 0 is C0

    ! Time-stepping loop
    do j = 1, nt
        ! Copy current concentration to the new array
        C_new = C
        
        ! Compute new concentration using explicit finite difference scheme
        do i = 2, nx - 1
            C_new(i) = C(i) + alpha * (C(i-1) - 2.0d0 * C(i) + C(i+1))
        end do
        
        ! Apply boundary conditions
        C_new(1) = 1.0d0    ! C0 at x = 0
        C_new(nx) = 0.0d0   ! Zero concentration at x = L (can be adjusted)

        ! Update concentration for the next time step
        C = C_new

        ! Output concentration profile at selected times
        if (mod(j, 100) == 0) then
            time = j * dt
            call print_concentration(C, nx, time)
        end if
    end do

end program mass_transfer_diffusion


! Subroutine to print concentration profile
subroutine print_concentration(C, nx, time)
    implicit none
    real(8), dimension(nx) :: C
    integer :: nx, i
    real(8) :: time

    print *, "Time = ", time
    do i = 1, nx
        print *, "x = ", (i-1) * (1.0d0 / (nx-1)), " C(x,t) = ", C(i)
    end do
end subroutine print_concentration