Microchannel Heat Sink Sizer (MCHE)

program microchannel_heat_exchanger_mche
    implicit none
    integer :: iostat_val, coolant_type, substrate_mat
    double precision :: Q_chip_W, W_footprint_mm, L_footprint_mm, W_channel_um, H_channel_um, W_fin_um
    double precision :: V_dot_mlpm, Tin_C, k_substrate, t_base_mm
    double precision :: rho_c, cp_c, mu_c, k_c, Pr_c, m_dot_kgs, N_channels, W_c_m, H_c_m, W_w_m
    double precision :: D_h_m, alpha_c, A_flow_tot, V_fluid, Re_ch, Po_num, f_darcy, Nu_ch, h_conv
    double precision :: m_fin, eta_fin, A_base_m2, A_eff_m2, R_conv, R_cond, R_cal, R_total_K_W
    double precision :: T_junction_C, dp_channel_kPa, heat_flux_Wcm2
    character(len=32) :: cooling_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) coolant_type    ! 1=Deionized Water, 2=50/50 Water-Glycol, 3=Dielectric Fluorinert (FC-72)
    read(*,*,iostat=iostat_val) substrate_mat   ! 1=Copper (k=398 W/m.K), 2=Silicon (k=148 W/m.K), 3=Aluminum (k=205 W/m.K)
    read(*,*,iostat=iostat_val) Q_chip_W        ! Heat Dissipation [W] (e.g. 350.0)
    read(*,*,iostat=iostat_val) W_footprint_mm  ! Chip Width W [mm] (e.g. 20.0)
    read(*,*,iostat=iostat_val) L_footprint_mm  ! Chip Length L [mm] (e.g. 20.0)
    read(*,*,iostat=iostat_val) W_channel_um    ! Microchannel Width [um] (e.g. 150.0)
    read(*,*,iostat=iostat_val) H_channel_um    ! Microchannel Height [um] (e.g. 600.0)
    read(*,*,iostat=iostat_val) W_fin_um        ! Microfin Pitch/Width [um] (e.g. 100.0)
    read(*,*,iostat=iostat_val) V_dot_mlpm      ! Coolant Flow Rate [mL/min] (e.g. 800.0)
    read(*,*,iostat=iostat_val) Tin_C           ! Coolant Inlet Temp [deg C] (e.g. 25.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for microchannel heat sink calculation.'
        stop
    end if

    if (Q_chip_W <= 0.0d0 .or. W_channel_um <= 0.0d0 .or. H_channel_um <= 0.0d0 .or. V_dot_mlpm <= 0.0d0) then
        write(*,*) 'ERROR: Thermal power, microchannel dimensions and flow rate must be positive.'
        stop
    end if

    ! Coolant properties
    if (coolant_type == 1) then ! Water
        rho_c = 997.0d0; cp_c = 4180.0d0; mu_c = 0.00089d0; k_c = 0.61d0
    else if (coolant_type == 2) then ! 50/50 Water-Glycol
        rho_c = 1060.0d0; cp_c = 3300.0d0; mu_c = 0.00350d0; k_c = 0.38d0
    else ! Dielectric FC-72
        rho_c = 1680.0d0; cp_c = 1100.0d0; mu_c = 0.00064d0; k_c = 0.057d0
    end if

    if (substrate_mat == 1) then
        k_substrate = 398.0d0 ! Copper
    else if (substrate_mat == 2) then
        k_substrate = 148.0d0 ! Silicon
    else
        k_substrate = 205.0d0 ! Aluminum
    end if

    t_base_mm = 1.0d0
    W_c_m = W_channel_um * 1.0d-6
    H_c_m = H_channel_um * 1.0d-6
    W_w_m = W_fin_um * 1.0d-6

    N_channels = floor((W_footprint_mm * 1.0d-3) / (W_c_m + W_w_m))
    D_h_m = (2.0d0 * W_c_m * H_c_m) / (W_c_m + H_c_m)
    alpha_c = H_c_m / W_c_m

    m_dot_kgs = (V_dot_mlpm * 1.0d-6 / 60.0d0) * rho_c
    A_flow_tot = N_channels * (W_c_m * H_c_m)
    V_fluid = (V_dot_mlpm * 1.0d-6 / 60.0d0) / max(1.0d-9, A_flow_tot)

    Re_ch = (rho_c * V_fluid * D_h_m) / mu_c
    Pr_c = (cp_c * mu_c) / k_c

    ! Laminar rectangular channel Poiseuille and Nusselt numbers (Kwak correlation)
    Po_num = 96.0d0 * (1.0d0 - 1.3553d0/alpha_c + 1.9467d0/(alpha_c**2) - 1.7012d0/(alpha_c**3))
    f_darcy = Po_num / max(1.0d0, Re_ch)

    Nu_ch = 8.235d0 * (1.0d0 - 2.0421d0/alpha_c + 3.0853d0/(alpha_c**2) - 2.4765d0/(alpha_c**3))
    Nu_ch = max(4.0d0, Nu_ch)

    h_conv = Nu_ch * k_c / D_h_m

    ! Fin efficiency
    m_fin = sqrt((2.0d0 * h_conv) / (k_substrate * W_w_m))
    eta_fin = tanh(m_fin * H_c_m) / (m_fin * H_c_m)

    A_base_m2 = (W_footprint_mm * 1.0d-3) * (L_footprint_mm * 1.0d-3)
    A_eff_m2 = N_channels * (L_footprint_mm * 1.0d-3) * (W_c_m + 2.0d0 * eta_fin * H_c_m)

    ! Thermal Resistances [K/W]
    R_conv = 1.0d0 / (h_conv * A_eff_m2)
    R_cond = (t_base_mm * 1.0d-3) / (k_substrate * A_base_m2)
    R_cal = 1.0d0 / (m_dot_kgs * cp_c)
    R_total_K_W = R_conv + R_cond + R_cal

    T_junction_C = Tin_C + Q_chip_W * R_total_K_W
    heat_flux_Wcm2 = Q_chip_W / ((W_footprint_mm/10.0d0) * (L_footprint_mm/10.0d0))

    ! Pressure drop across microchannels [kPa]
    dp_channel_kPa = (f_darcy * ((L_footprint_mm * 1.0d-3) / D_h_m) * (rho_c * V_fluid**2 / 2.0d0)) / 1000.0d0

    if (T_junction_C < 70.0d0) then
        cooling_status = 'EXCELLENT (TJ < 70°C)'
    else if (T_junction_C < 85.0d0) then
        cooling_status = 'SAFE OPERATING LIMIT'
    else
        cooling_status = 'WARNING: THERMAL THROTTLING RISK'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — MICROCHANNEL HEAT EXCHANGER (MCHE)'
    write(*,'(A,F8.1,A,I5,A,I5,A)') 'Channels Dimensions          = ', W_channel_um, 'x', &
        int(H_channel_um), ' um (', int(N_channels), ' channels)'
    write(*,'(A,F10.1,A,F10.2,A)') 'Hydraulic Diam / Re Number   = ', D_h_m*1.0d6, ' um / ', Re_ch, ' (Laminar)'
    write(*,'(A,F10.1,A)')         'Heat Flux Density            = ', heat_flux_Wcm2, ' W/cm2'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,A)')       'CHIP JUNCTION TEMPERATURE Tj = ', T_junction_C, ' deg C | ', trim(cooling_status)
    write(*,'(A,F10.4,A)')         'Total Thermal Resistance Rth = ', R_total_K_W, ' K/W (C/W)'
    write(*,'(A,F10.1,A)')         'Convective Heat Coeff h_conv = ', h_conv, ' W/(m2.K)'
    write(*,'(A,F10.2,A)')         'Microchannel Pressure Drop   = ', dp_channel_kPa, ' kPa'
    write(*,'(A)') '============================================================'

end program microchannel_heat_exchanger_mche