💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Variable Conductivity k = f(T)
Core Numerical Engine in Fortran 90 • 48 total downloads
dependent_k.f90
! =========================================================================
! Source File: dependent_k.f90
! =========================================================================
program dependent_k
implicit none
! Constants
real, parameter :: PI = 3.141592653589793
! Inputs
integer :: geom ! 1 = Plane, 2 = Cylinder, 3 = Sphere
real :: L_thick, Area ! Plane wall parameters (thickness, area)
real :: r1, r2, H_len ! Cylinder/Sphere parameters (radii, length)
integer :: model ! 1 = Linear, 2 = Polynomial, 3 = Tabulated
real :: k0, beta ! Linear parameters
real :: poly_a, poly_b, poly_c ! Polynomial parameters: k = a + b*T + c*T^2
integer :: num_pts ! Number of tabulated points
real, allocatable :: tab_T(:), tab_k(:), tab_theta(:) ! Tabulated arrays
real :: T1, T2 ! Boundary temperatures (C)
! Outputs & Solved variables
real :: theta1, theta2
real :: Q, k_eff
real :: T_node(11), x_node(11), theta_node(11)
real :: Q_const, k_const, T_const(11)
! Temporary/Iterative variables
integer :: i, j, temp_idx
real :: temp_val, k_val, m_val, dT
real :: r_val, pos_ratio
real :: f_val, df_val, T_guess
! 1. Read Inputs from stdin
read(*,*) geom
if (geom == 1) then
read(*,*) L_thick
read(*,*) Area
else if (geom == 2) then
read(*,*) r1
read(*,*) r2
read(*,*) H_len
else if (geom == 3) then
read(*,*) r1
read(*,*) r2
end if
read(*,*) model
if (model == 1) then
read(*,*) k0
read(*,*) beta
else if (model == 2) then
read(*,*) poly_a
read(*,*) poly_b
read(*,*) poly_c
else if (model == 3) then
read(*,*) num_pts
allocate(tab_T(num_pts))
allocate(tab_k(num_pts))
allocate(tab_theta(num_pts))
do i = 1, num_pts
read(*,*) tab_T(i)
read(*,*) tab_k(i)
end do
! Sort tabulated points in increasing order of T (simple bubble sort)
do i = 1, num_pts - 1
do j = i + 1, num_pts
if (tab_T(i) > tab_T(j)) then
temp_val = tab_T(i)
tab_T(i) = tab_T(j)
tab_T(j) = temp_val
temp_val = tab_k(i)
tab_k(i) = tab_k(j)
tab_k(j) = temp_val
end if
end do
end do
! Compute cumulative integrals (theta) at tabulated temperatures
tab_theta(1) = 0.0
do i = 2, num_pts
tab_theta(i) = tab_theta(i-1) + 0.5 * (tab_k(i-1) + tab_k(i)) * (tab_T(i) - tab_T(i-1))
end do
end if
read(*,*) T1
read(*,*) T2
! 2. Compute Surface Thetas (theta1 and theta2)
theta1 = get_theta(T1)
theta2 = get_theta(T2)
! 3. Compute Heat Transfer Rate (Q) & Mean Effective Conductivity (k_eff)
if (geom == 1) then
Q = (Area / L_thick) * (theta1 - theta2)
k_eff = (theta1 - theta2) / (T1 - T2)
else if (geom == 2) then
Q = (2.0 * PI * H_len / log(r2 / r1)) * (theta1 - theta2)
k_eff = (theta1 - theta2) / (T1 - T2)
else if (geom == 3) then
Q = (4.0 * PI * r1 * r2 / (r2 - r1)) * (theta1 - theta2)
k_eff = (theta1 - theta2) / (T1 - T2)
end if
! 4. Generate Node Profiles (11 nodes from inner to outer boundary)
do i = 1, 11
pos_ratio = real(i - 1) / 10.0
if (geom == 1) then
x_node(i) = pos_ratio * L_thick
theta_node(i) = theta1 - pos_ratio * (theta1 - theta2)
else
r_val = r1 + pos_ratio * (r2 - r1)
x_node(i) = r_val
if (geom == 2) then
theta_node(i) = theta1 - (theta1 - theta2) * (log(r_val / r1) / log(r2 / r1))
else
theta_node(i) = theta1 - (theta1 - theta2) * ((1.0/r1 - 1.0/r_val) / (1.0/r1 - 1.0/r2))
end if
end if
T_node(i) = get_temp_from_theta(theta_node(i))
end do
! 5. Comparison with Constant-k Solver (evaluating at average temperature)
k_const = get_k_at_temp(0.5 * (T1 + T2))
if (geom == 1) then
Q_const = (k_const * Area / L_thick) * (T1 - T2)
do i = 1, 11
T_const(i) = T1 - (real(i - 1) / 10.0) * (T1 - T2)
end do
else if (geom == 2) then
Q_const = (2.0 * PI * k_const * H_len / log(r2 / r1)) * (T1 - T2)
do i = 1, 11
r_val = r1 + (real(i - 1) / 10.0) * (r2 - r1)
T_const(i) = T1 - (T1 - T2) * (log(r_val / r1) / log(r2 / r1))
end do
else if (geom == 3) then
Q_const = (4.0 * PI * k_const * r1 * r2 / (r2 - r1)) * (T1 - T2)
do i = 1, 11
r_val = r1 + (real(i - 1) / 10.0) * (r2 - r1)
T_const(i) = T1 - (T1 - T2) * ((1.0/r1 - 1.0/r_val) / (1.0/r1 - 1.0/r2))
end do
end if
! 6. Output Engineering Report
write(*,*) "=========================================================================="
write(*,*) " STEADY CONDUCTION WITH VARIABLE PROPERTIES k = f(T) "
write(*,*) "=========================================================================="
if (geom == 1) then
write(*, '(A)') " Geometry: PLANE WALL"
write(*, '(A, F10.4, A, F10.4, A)') " Dimensions: Thickness = ", L_thick, " m, Area = ", Area, " m2"
else if (geom == 2) then
write(*, '(A)') " Geometry: CYLINDER"
write(*, '(A, F10.4, A, F10.4, A, F10.4, A)') " Dimensions: r1 = ", r1, " m, r2 = ", r2, " m, Length = ", H_len, " m"
else if (geom == 3) then
write(*, '(A)') " Geometry: SPHERE"
write(*, '(A, F10.4, A, F10.4, A)') " Dimensions: r1 = ", r1, " m, r2 = ", r2, " m"
end if
if (model == 1) then
write(*, '(A)') " Property Model: Linear k(T) = k0*(1 + beta*T)"
write(*, '(A, F10.4, A, E13.4, A)') " Parameters: k0 = ", k0, " W/(m-K), beta = ", beta, " 1/K"
else if (model == 2) then
write(*, '(A)') " Property Model: Polynomial k(T) = a + b*T + c*T^2"
write(*, '(A, F10.4, A, E13.4, A, E13.4, A)') " Parameters: a = ", poly_a, ", b = ", poly_b, ", c = ", poly_c
else if (model == 3) then
write(*, '(A)') " Property Model: Tabulated / Interpolated points"
write(*, '(A, I3)') " Tabulated points: ", num_pts
do i = 1, num_pts
write(*, '(A, I2, A, F8.2, A, F8.3)') " Point ", i, ": T = ", tab_T(i), " C, k = ", tab_k(i)
end do
end if
write(*, '(A, F8.2, A, F8.2, A)') " Boundary Temperatures: T1 = ", T1, " C, T2 = ", T2, " C"
write(*,*) "--------------------------------------------------------------------------"
write(*,*) " THERMAL PERFORMANCE ANALYSIS"
write(*,*) "--------------------------------------------------------------------------"
write(*, '(A, F14.2, A)') " Heat Transfer Rate (Kirchhoff Q): ", Q, " W"
write(*, '(A, F14.2, A)') " Heat Transfer Rate (Constant k_avg Q):", Q_const, " W"
write(*, '(A, F14.2, A)') " Absolute Difference in Q: ", abs(Q - Q_const), " W"
write(*, '(A, F14.4, A)') " Mean Effective Conductivity (k_eff): ", k_eff, " W/(m-K)"
write(*, '(A, F14.4, A)') " Constant Conductivity (k_avg): ", k_const, " W/(m-K)"
write(*,*) "--------------------------------------------------------------------------"
write(*,*) " SPATIAL PROFILE DATA"
write(*,*) "--------------------------------------------------------------------------"
if (geom == 1) then
write(*,*) " Node Position x(m) Transformed theta Temp T(C) Temp Constant-k(C)"
else
write(*,*) " Node Radius r(m) Transformed theta Temp T(C) Temp Constant-k(C)"
end if
do i = 1, 11
write(*, '(I4.2, F16.4, F20.4, F12.2, F21.2)') &
i, x_node(i), theta_node(i), T_node(i), T_const(i)
end do
write(*,*) "=========================================================================="
contains
! Helper: evaluates theta(T) = \int_0^T k(T') dT'
real function get_theta(temp)
real, intent(in) :: temp
integer :: k_idx
real :: m_local, dT_local
if (model == 1) then
get_theta = k0 * (temp + 0.5 * beta * temp**2)
else if (model == 2) then
get_theta = poly_a * temp + 0.5 * poly_b * temp**2 + (poly_c / 3.0) * temp**3
else if (model == 3) then
! For tabulated model, T_ref = tab_T(1)
if (temp <= tab_T(1)) then
! Extrapolate using first interval gradient
m_local = (tab_k(2) - tab_k(1)) / (tab_T(2) - tab_T(1))
dT_local = temp - tab_T(1)
get_theta = tab_k(1) * dT_local + 0.5 * m_local * dT_local**2
else if (temp >= tab_T(num_pts)) then
! Extrapolate using last interval gradient
m_local = (tab_k(num_pts) - tab_k(num_pts-1)) / (tab_T(num_pts) - tab_T(num_pts-1))
dT_local = temp - tab_T(num_pts)
get_theta = tab_theta(num_pts) + tab_k(num_pts) * dT_local + 0.5 * m_local * dT_local**2
else
! Find interval
k_idx = 1
do while (tab_T(k_idx+1) < temp .and. k_idx < num_pts - 1)
k_idx = k_idx + 1
end do
m_local = (tab_k(k_idx+1) - tab_k(k_idx)) / (tab_T(k_idx+1) - tab_T(k_idx))
dT_local = temp - tab_T(k_idx)
get_theta = tab_theta(k_idx) + tab_k(k_idx) * dT_local + 0.5 * m_local * dT_local**2
end if
end if
end function get_theta
! Helper: evaluates k(T) at a given temperature T
real function get_k_at_temp(temp)
real, intent(in) :: temp
integer :: k_idx
real :: m_local
if (model == 1) then
get_k_at_temp = k0 * (1.0 + beta * temp)
else if (model == 2) then
get_k_at_temp = poly_a + poly_b * temp + poly_c * temp**2
else if (model == 3) then
if (temp <= tab_T(1)) then
m_local = (tab_k(2) - tab_k(1)) / (tab_T(2) - tab_T(1))
get_k_at_temp = tab_k(1) + m_local * (temp - tab_T(1))
else if (temp >= tab_T(num_pts)) then
m_local = (tab_k(num_pts) - tab_k(num_pts-1)) / (tab_T(num_pts) - tab_T(num_pts-1))
get_k_at_temp = tab_k(num_pts) + m_local * (temp - tab_T(num_pts))
else
k_idx = 1
do while (tab_T(k_idx+1) < temp .and. k_idx < num_pts - 1)
k_idx = k_idx + 1
end do
m_local = (tab_k(k_idx+1) - tab_k(k_idx)) / (tab_T(k_idx+1) - tab_T(k_idx))
get_k_at_temp = tab_k(k_idx) + m_local * (temp - tab_T(k_idx))
end if
end if
end function get_k_at_temp
! Helper: inverts theta to retrieve T
real function get_temp_from_theta(theta_val)
real, intent(in) :: theta_val
real :: disc, root
integer :: iter, k_idx
real :: m_local, T_val, f_val, df_local
if (model == 1) then
if (beta == 0.0) then
get_temp_from_theta = theta_val / k0
else
disc = 1.0 + 2.0 * beta * theta_val / k0
if (disc >= 0.0) then
! Standard physical root
root = sqrt(disc)
if (1.0 + beta * 0.5 * (T1 + T2) >= 0.0) then
get_temp_from_theta = (-1.0 + root) / beta
else
get_temp_from_theta = (-1.0 - root) / beta
end if
else
get_temp_from_theta = theta_val / k0 ! fallback
end if
end if
else if (model == 2) then
! Solve a*T + (b/2)*T^2 + (c/3)*T^3 - theta = 0 using Newton's method
T_val = 0.5 * (T1 + T2) ! initial guess
iter = 0
do while (iter < 100)
iter = iter + 1
f_val = poly_a * T_val + 0.5 * poly_b * T_val**2 + (poly_c / 3.0) * T_val**3 - theta_val
df_local = poly_a + poly_b * T_val + poly_c * T_val**2
if (abs(df_local) < 1e-6) exit
dT = f_val / df_local
T_val = T_val - dT
if (abs(dT) < 1e-5) exit
end do
get_temp_from_theta = T_val
else if (model == 3) then
! Tabulated inversion
if (theta_val <= 0.0) then
! Extrapolate below tab_T(1)
m_local = (tab_k(2) - tab_k(1)) / (tab_T(2) - tab_T(1))
if (m_local == 0.0) then
get_temp_from_theta = tab_T(1) + theta_val / tab_k(1)
else
disc = tab_k(1)**2 + 2.0 * m_local * theta_val
if (disc >= 0.0) then
get_temp_from_theta = tab_T(1) + (-tab_k(1) + sqrt(disc)) / m_local
else
get_temp_from_theta = tab_T(1) + theta_val / tab_k(1)
end if
end if
else if (theta_val >= tab_theta(num_pts)) then
! Extrapolate above tab_T(num_pts)
m_local = (tab_k(num_pts) - tab_k(num_pts-1)) / (tab_T(num_pts) - tab_T(num_pts-1))
if (m_local == 0.0) then
get_temp_from_theta = tab_T(num_pts) + (theta_val - tab_theta(num_pts)) / tab_k(num_pts)
else
disc = tab_k(num_pts)**2 + 2.0 * m_local * (theta_val - tab_theta(num_pts))
if (disc >= 0.0) then
get_temp_from_theta = tab_T(num_pts) + (-tab_k(num_pts) + sqrt(disc)) / m_local
else
get_temp_from_theta = tab_T(num_pts) + (theta_val - tab_theta(num_pts)) / tab_k(num_pts)
end if
end if
else
! Find interval
k_idx = 1
do while (tab_theta(k_idx+1) < theta_val .and. k_idx < num_pts - 1)
k_idx = k_idx + 1
end do
m_local = (tab_k(k_idx+1) - tab_k(k_idx)) / (tab_T(k_idx+1) - tab_T(k_idx))
if (m_local == 0.0) then
get_temp_from_theta = tab_T(k_idx) + (theta_val - tab_theta(k_idx)) / tab_k(k_idx)
else
disc = tab_k(k_idx)**2 + 2.0 * m_local * (theta_val - tab_theta(k_idx))
if (disc >= 0.0) then
get_temp_from_theta = tab_T(k_idx) + (-tab_k(k_idx) + sqrt(disc)) / m_local
else
get_temp_from_theta = tab_T(k_idx) + (theta_val - tab_theta(k_idx)) / tab_k(k_idx)
end if
end if
end if
end if
end function get_temp_from_theta
end program dependent_k
Solver Description
Calculate steady heat conduction with variable thermal conductivity k=f(T) across plane walls, cylinders, and spheres using the Kirchhoff transformation.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 dependent_k.f90 -o dependent_k_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
dependent_k_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry (1=Slab, 2=Cylinder, 3=Sphere)
1
! Slab thickness L [m]
0.1
! Slab area A [m2]
1.0
! Conductivity model (1=Linear, 2=Quadratic, 3=Tabular)
1
! Reference conductivity k0 [W/m-K]
1.5
! Temperature coefficient beta [1/K]
0.003
! Boundary 1 temperature T1 [°C]
300.0
! Boundary 2 temperature T2 [°C]
50.0
1
! Slab thickness L [m]
0.1
! Slab area A [m2]
1.0
! Conductivity model (1=Linear, 2=Quadratic, 3=Tabular)
1
! Reference conductivity k0 [W/m-K]
1.5
! Temperature coefficient beta [1/K]
0.003
! Boundary 1 temperature T1 [°C]
300.0
! Boundary 2 temperature T2 [°C]
50.0