program natural_convection implicit none ! Inputs integer :: geom_type double precision :: dim1, dim2 double precision :: T_s, T_inf double precision :: rho, mu, k_f, Pr, Cp, beta ! Constants double precision, parameter :: g = 9.80665d0 double precision, parameter :: PI = 3.141592653589793d0 ! Physical properties double precision :: nu, L_c, A_s, T_f, dT double precision :: Gr, Ra, Nu_avg, h_avg, Q ! Local profile variables integer :: i, n_points double precision :: x, dx, local_Ra, local_Gr, local_Nu, local_h, local_delta character(len=50) :: geom_name integer :: iostat_val ! Read inputs read(*,*,iostat=iostat_val) geom_type if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid geometry type input.' stop end if read(*,*,iostat=iostat_val) dim1 read(*,*,iostat=iostat_val) dim2 read(*,*,iostat=iostat_val) T_s read(*,*,iostat=iostat_val) T_inf read(*,*,iostat=iostat_val) rho read(*,*,iostat=iostat_val) mu read(*,*,iostat=iostat_val) k_f read(*,*,iostat=iostat_val) Pr read(*,*,iostat=iostat_val) Cp read(*,*,iostat=iostat_val) beta if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all convection parameters.' stop end if ! Basic validations if (dim1 <= 0.0d0) then write(*,*) 'ERROR: Dimension 1 must be positive.' stop end if if (rho <= 0.0d0 .or. mu <= 0.0d0 .or. k_f <= 0.0d0 .or. Pr <= 0.0d0 .or. Cp <= 0.0d0) then write(*,*) 'ERROR: Fluid properties must be positive.' stop end if ! Compute film temperature and kinematic viscosity T_f = (T_s + T_inf) / 2.0d0 nu = mu / rho dT = abs(T_s - T_inf) ! Setup Geometry parameters select case (geom_type) case (1) geom_name = "Vertical Plate" L_c = dim1 A_s = dim1 * dim2 case (2) geom_name = "Horizontal Plate (Upper Hot / Lower Cold)" if (dim2 <= 0.0d0) dim2 = dim1 ! Default to square if W is zero L_c = (dim1 * dim2) / (2.0d0 * (dim1 + dim2)) A_s = dim1 * dim2 case (3) geom_name = "Horizontal Plate (Lower Hot / Upper Cold)" if (dim2 <= 0.0d0) dim2 = dim1 L_c = (dim1 * dim2) / (2.0d0 * (dim1 + dim2)) A_s = dim1 * dim2 case (4) geom_name = "Horizontal Cylinder" if (dim2 <= 0.0d0) dim2 = 1.0d0 ! Default to 1m length L_c = dim1 A_s = PI * dim1 * dim2 case (5) geom_name = "Sphere" L_c = dim1 A_s = PI * dim1**2 case default write(*,*) 'ERROR: Invalid geometry type selected.' stop end select ! Grashof and Rayleigh numbers Gr = (g * beta * dT * L_c**3) / (nu**2) Ra = Gr * Pr ! Nusselt Number Calculations select case (geom_type) case (1) ! Churchill and Chu (1975) for Vertical Plate (All Ra) Nu_avg = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / & ((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2 case (2) ! Horizontal Plate (Upper Hot / Lower Cold) if (Ra < 1.0d4) then Nu_avg = 0.54d0 * Ra**0.25d0 ! Extrapolated/lower-laminar else if (Ra <= 1.0d7) then Nu_avg = 0.54d0 * Ra**0.25d0 else if (Ra <= 1.0d11) then Nu_avg = 0.15d0 * Ra**(1.0d0/3.0d0) else Nu_avg = 0.15d0 * Ra**(1.0d0/3.0d0) ! Extrapolated/upper-turbulent end if case (3) ! Horizontal Plate (Lower Hot / Upper Cold) ! Valid for 10^5 to 10^11 Nu_avg = 0.27d0 * Ra**0.25d0 case (4) ! Churchill and Chu (1975) for Horizontal Cylinder (Ra <= 10^12) Nu_avg = (0.60d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / & ((1.0d0 + (0.559d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2 case (5) ! Churchill (1983) for Sphere (Ra <= 10^11 and Pr >= 0.7) Nu_avg = 2.0d0 + (0.589d0 * Ra**0.25d0) / & ((1.0d0 + (0.469d0/Pr)**(9.0d0/16.0d0))**(4.0d0/9.0d0)) end select ! Average heat transfer coefficient h_avg = Nu_avg * k_f / L_c ! Heat transfer rate Q = h_avg * A_s * dT ! Output results write(*,'(A)') '============================================================' write(*,'(A)') ' NATURAL CONVECTION HEAT TRANSFER ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A,A)') ' Geometry Name = ', trim(geom_name) write(*,'(A,I2)') ' Geometry Type Code = ', geom_type write(*,'(A,F12.4,A)') ' Characteristic Length = ', L_c, ' m' write(*,'(A,F12.4,A)') ' Surface Area (As) = ', A_s, ' m2' write(*,'(A,F12.2,A)') ' Surface Temperature = ', T_s, ' deg-C' write(*,'(A,F12.2,A)') ' Ambient Temperature = ', T_inf, ' deg-C' write(*,'(A,F12.2,A)') ' Film Temperature = ', T_f, ' deg-C' write(*,*) write(*,'(A)') '--- FLUID PROPERTIES ----------------------------------------' write(*,'(A,ES12.4,A)') ' Density (rho) = ', rho, ' kg/m3' write(*,'(A,ES12.4,A)') ' Dynamic Viscosity (mu) = ', mu, ' Pa.s' write(*,'(A,ES12.4,A)') ' Kinematic Viscosity (nu) = ', nu, ' m2/s' write(*,'(A,F12.4,A)') ' Thermal Conductivity = ', k_f, ' W/m.K' write(*,'(A,F12.4)') ' Prandtl Number (Pr) = ', Pr write(*,'(A,F12.2,A)') ' Specific Heat (Cp) = ', Cp, ' J/kg.K' write(*,'(A,ES12.4,A)') ' Expansion Coeff (beta) = ', beta, ' 1/K' write(*,*) write(*,'(A)') '--- CONVECTION RESULTS --------------------------------------' write(*,'(A,ES12.4)') ' Grashof Number (Gr) = ', Gr write(*,'(A,ES12.4)') ' Rayleigh Number (Ra) = ', Ra write(*,'(A,F12.4)') ' Average Nusselt Number = ', Nu_avg write(*,'(A,F12.4,A)') ' Average h = ', h_avg, ' W/m2.K' write(*,'(A,F12.4,A)') ' Heat Transfer Rate (Q) = ', Q, ' W' write(*,*) ! ============================================ ! LOCAL PROFILE DATA (along coordinate x) ! ============================================ write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------' write(*,'(A)') ' x [m] Ra_x Nu_x h_x [W/m2.K] d [mm]' write(*,'(A)') ' ----------------------------------------------------------' n_points = 40 dx = L_c / dble(n_points) do i = 1, n_points x = dble(i) * dx local_Gr = (g * beta * dT * x**3) / (nu**2) local_Ra = local_Gr * Pr if (dT <= 1.0d-10 .or. local_Gr <= 1.0d-10) then local_Nu = 0.0d0 local_delta = 0.0d0 else if (local_Ra < 1.0d9) then ! Laminar boundary layer correlations local_Nu = 0.508d0 * (Pr / (0.952d0 + Pr))**0.25d0 * local_Ra**0.25d0 local_delta = 3.93d0 * x * ((0.952d0 + Pr) / (Pr**2 * local_Gr))**0.25d0 else ! Turbulent boundary layer approximations local_Nu = 0.0292d0 * local_Ra**0.4d0 local_delta = 0.15d0 * x * local_Ra**(-0.1d0) end if end if if (x > 0.0d0) then local_h = local_Nu * k_f / x else local_h = 0.0d0 end if write(*,'(F8.4,2X,ES12.4,2X,F10.2,4X,F10.4,4X,F10.4)') & x, local_Ra, local_Nu, local_h, local_delta*1000.0d0 end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' select case (geom_type) case (1) write(*,'(A)') ' Vertical Plate: Churchill and Chu (1975) Correlation' case (2) write(*,'(A)') ' Horizontal Plate (Upper Hot): Nu = 0.54 Ra^0.25 (lam), Nu = 0.15 Ra^(1/3) (turb)' case (3) write(*,'(A)') ' Horizontal Plate (Lower Hot): Nu = 0.27 Ra^0.25' case (4) write(*,'(A)') ' Horizontal Cylinder: Churchill and Chu (1975) Correlation' case (5) write(*,'(A)') ' Sphere: Churchill (1983) Correlation' end select write(*,'(A)') ' Boundary Layer: Laminar similarity solutions for natural convection.' end program natural_convection