!============================================================================== ! ThermoFluidCalc — Calculator #20 : FVM Time-Step Estimator !============================================================================== ! Physics : In the Finite-Volume Method the time step must satisfy ! ! Convective limit : Δt_conv ≤ V / Σ_f |v_{n,f} · S_f| ! Diffusive limit : Δt_diff ≤ V / (2 · α · Σ_f (S_f / d_f)) ! Combined : Δt = min(Δt_conv, Δt_diff) × safety ! ! where ! V = cell volume (m³ or m² in 2-D) ! S_f = face area (m² or m in 2-D) ! v_{n,f} = velocity component normal to face f (m/s) ! d_f = distance from cell centroid to neighbour centroid through f (m) ! α = diffusivity (m²/s) ! safety = user-specified safety factor (0 < safety ≤ 1) ! ! Modes: ! 1 = Quick single-face : Δt = V / (|v|·S) × safety ! 2 = Convective (multi) : Δt = V / Σ|v_n·S_f| × safety ! 3 = Diffusive (multi) : Δt = V / (2·α·Σ(S_f/d_f)) × safety ! 4 = Combined (multi) : min of modes 2 and 3 ! ! Reference : Versteeg & Malalasekera, "An Introduction to CFD — The FVM" ! Moukalled, Mangani & Darwiche, "The FVM in CFD", §7 ! ! Build: ! gfortran -O2 -o fvm_timestep fvm_timestep.f90 ! ! Input (stdin): ! Line 1 : mode safety ! Mode 1 line 2 : V v S ! Mode 2-4 line 2 : V nfaces alpha ! then nfaces lines : S_f v_nf d_f ! ! Output (stdout): structured key=value + per-face data !============================================================================== program fvm_timestep implicit none ! --- Double precision kind ------------------------------------------------- integer, parameter :: dp = selected_real_kind(15, 307) integer, parameter :: MAX_FACES = 20 ! --- Variables ------------------------------------------------------------- integer :: mode, nfaces, i real(dp) :: safety, V, v_single, S_single real(dp) :: alpha real(dp) :: Sf(MAX_FACES), vnf(MAX_FACES), df(MAX_FACES) real(dp) :: flux_f, sum_flux, sum_diff_ratio real(dp) :: dt_conv, dt_diff, dt_combined character(len=40) :: mode_name, limiting ! --- Read mode and safety -------------------------------------------------- read(*,*) mode, safety ! Clamp safety if (safety <= 0.0_dp) safety = 0.1_dp if (safety > 1.0_dp) safety = 1.0_dp ! ── Mode 1 : Quick single-face ──────────────────────────────────────────── if (mode == 1) then read(*,*) V, v_single, S_single if (V <= 0.0_dp .or. S_single <= 0.0_dp) then write(*,'(A)') 'ERROR=Volume and face area must be positive.' stop end if mode_name = 'Quick Single-Face' sum_flux = abs(v_single) * S_single dt_conv = 0.0_dp dt_diff = 0.0_dp dt_combined = 0.0_dp sum_diff_ratio = 0.0_dp if (sum_flux > 0.0_dp) then dt_conv = V / sum_flux * safety end if write(*,'(A,I2)') 'MODE=', mode write(*,'(A,A)') 'MODE_NAME=', trim(mode_name) write(*,'(A,ES15.8)') 'CELL_VOLUME=', V write(*,'(A,ES15.8)') 'SAFETY=', safety write(*,'(A,I2)') 'NFACES=', 1 write(*,'(A,ES15.8,A,ES15.8,A,ES15.8)') & 'FACE_1: S_f=', S_single, ', v_n=', v_single, ', flux=', sum_flux write(*,'(A,ES15.8)') 'SUM_FLUX=', sum_flux write(*,'(A,ES15.8)') 'DT_CONV=', dt_conv write(*,'(A)') 'DT_DIFF=N/A' write(*,'(A,ES15.8)') 'DT_RESULT=', dt_conv write(*,'(A)') 'LIMITING=convective' ! ── Modes 2-4 : Multi-face ──────────────────────────────────────────────── else if (mode >= 2 .and. mode <= 4) then read(*,*) V, nfaces, alpha if (V <= 0.0_dp) then write(*,'(A)') 'ERROR=Cell volume must be positive.' stop end if if (nfaces < 1 .or. nfaces > MAX_FACES) then write(*,'(A,I3,A)') 'ERROR=Number of faces must be 1-', MAX_FACES, '.' stop end if ! Read face data do i = 1, nfaces read(*,*) Sf(i), vnf(i), df(i) end do ! Compute sums sum_flux = 0.0_dp sum_diff_ratio = 0.0_dp ! Print header select case (mode) case (2); mode_name = 'Convective (Multi-Face)' case (3); mode_name = 'Diffusive (Multi-Face)' case (4); mode_name = 'Combined (Multi-Face)' end select write(*,'(A,I2)') 'MODE=', mode write(*,'(A,A)') 'MODE_NAME=', trim(mode_name) write(*,'(A,ES15.8)') 'CELL_VOLUME=', V write(*,'(A,ES15.8)') 'ALPHA=', alpha write(*,'(A,ES15.8)') 'SAFETY=', safety write(*,'(A,I3)') 'NFACES=', nfaces write(*,'(A)') 'FACES_START' do i = 1, nfaces flux_f = abs(vnf(i)) * Sf(i) sum_flux = sum_flux + flux_f if (df(i) > 0.0_dp) then sum_diff_ratio = sum_diff_ratio + Sf(i) / df(i) end if write(*,'(I3,A,ES15.8,A,ES15.8,A,ES15.8,A,ES15.8)') & i, ',', Sf(i), ',', vnf(i), ',', df(i), ',', flux_f end do write(*,'(A)') 'FACES_END' write(*,'(A,ES15.8)') 'SUM_FLUX=', sum_flux write(*,'(A,ES15.8)') 'SUM_DIFF_RATIO=', sum_diff_ratio ! Compute Δt values dt_conv = 0.0_dp dt_diff = 0.0_dp if (sum_flux > 0.0_dp) then dt_conv = V / sum_flux * safety end if if (alpha > 0.0_dp .and. sum_diff_ratio > 0.0_dp) then dt_diff = V / (2.0_dp * alpha * sum_diff_ratio) * safety end if select case (mode) case (2) write(*,'(A,ES15.8)') 'DT_CONV=', dt_conv write(*,'(A)') 'DT_DIFF=N/A' dt_combined = dt_conv limiting = 'convective' case (3) write(*,'(A)') 'DT_CONV=N/A' write(*,'(A,ES15.8)') 'DT_DIFF=', dt_diff dt_combined = dt_diff limiting = 'diffusive' case (4) write(*,'(A,ES15.8)') 'DT_CONV=', dt_conv write(*,'(A,ES15.8)') 'DT_DIFF=', dt_diff if (dt_conv > 0.0_dp .and. dt_diff > 0.0_dp) then dt_combined = min(dt_conv, dt_diff) if (abs(dt_conv - dt_diff) < 1.0e-15_dp) then limiting = 'equal' else if (dt_conv < dt_diff) then limiting = 'convective' else limiting = 'diffusive' end if else if (dt_conv > 0.0_dp) then dt_combined = dt_conv limiting = 'convective' else if (dt_diff > 0.0_dp) then dt_combined = dt_diff limiting = 'diffusive' else dt_combined = 0.0_dp limiting = 'none' end if end select write(*,'(A,ES15.8)') 'DT_RESULT=', dt_combined write(*,'(A,A)') 'LIMITING=', trim(limiting) else write(*,'(A)') 'ERROR=Invalid mode (must be 1-4).' stop end if end program fvm_timestep