Coax heat transfer analysis

coax_thermal.jl

using DifferentialEquations
using BoundaryValueDiffEq
using DiffEqCallbacks
using ForwardDiff
using Plots

struct Material
    k_coefs::Vector
end

k(m::Material,T) = 10^sum([coef * log10(T)^(i - 1) for (i, coef) ∈ enumerate(m.k_coefs)])

# From the NIST cryogenic material properties
stainless316 = Material([-1.4087, 1.3982, 0.2543, -0.6260, 0.2334, 0.4256, -0.4658, 0.1650, -0.0199])

function temp_profile(material, T_hot, T_cold, L, dL)
    # Avoid inexact errors by forcing things to be floats
    bv = [Float64(T_hot), Float64(T_cold)]
    therm_cond(T) = k(material, T)
    therm_cond′(T) = ForwardDiff.derivative(therm_cond, T)
    # Define problem
    function temp_on_beam!(du, u, _, _)
        T = u[1]
        dT = u[2]
        du[1] = dT
        du[2] = -1 / therm_cond(T) * therm_cond′(T) * dT^2
    end
    # Set boundary conditions
    function bc!(residual, u, _, _)
        residual[1] = u[1][1] - bv[1]
        residual[2] = u[end][1] - bv[2]
    end
    # Solve
    bvp = TwoPointBVProblem(temp_on_beam!, bc!, bv, (0.0, L))
    sol = solve(bvp, MIRK4(); dt=dL)
    # Return solution
    sol.t, map(first, sol.u)
end

x, temp = temp_profile(stainless316, 300, 77, 5.90, 0.1)

plot(x,temp,label="DE Solution",xlabel="Position (cm)",ylabel="Temperature (K)",title="Temperature Profile")
plot!([0,5.9],[300,77],label="Linear Solution")