<jupyter_start><jupyter_text>Transient Heat Conduction in Various GeometriesTransport of heat in a solid is described by the familiar thermal diffusion model\begin{align*}\rho C_p\frac{\partial T}{\partial t} & = \nabla\cdot(k\nabla T)\end{align*} RescalingWe'll assume the thermal conductivity $k$ is a constant, and define thermal diffusivity in the conventional way\begin{align*}\alpha & = \frac{k}{\rho C_p}\end{align*}We will further assume symmetry with respect to all spatial coordinates except $r$ where $r$ extends from $-R$ to $+R$. The boundary conditions are\begin{align*}T(t,R) & = T_{\infty} & \forall t > 0 \\\nabla T(t,0) & = 0 & \forall t \geq 0 \end{align*}where we have assumed symmetry with respect to $r$ and uniform initial conditions $T(0, r) = T_0$ for all $0 \leq r \leq R$. Following standard scaling procedures, we introduce the dimensionless variables\begin{align*}T' & = \frac{T - T_0}{T_\infty - T_0} \\r' & = \frac{r}{R} \\t' & = t \frac{\alpha}{R^2}\end{align*} Dimensionless ModelUnder these conditions the problem reduces to \begin{align*}\frac{\partial T'}{\partial t'} & = \nabla^2 T'\end{align*}with auxiliary conditions\begin{align*}T'(0, r') & = 0 & \forall 0 \leq r' \leq 1\\T'(t', 1) & = 1 & \forall t' > 0\\\nabla T'(t', 0) & = 0 & \forall t' \geq 0 \\\end{align*}which we can specialize to specific geometries.  Preliminary Code<jupyter_code>!pip install -q pyomo
!wget -N -q "https://ampl.com/dl/open/ipopt/ipopt-linux64.zip"
!unzip -o -q ipopt-linux64

ipopt_executable='/content/ipopt'
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

def model_plot(m):
    r = sorted(m.r)
    t = sorted(m.t)

    rgrid = np.zeros((len(t), len(r)))
    tgrid = np.zeros((len(t), len(r)))
    Tgrid = np.zeros((len(t), len(r)))

    for i in range(0, len(t)):
        for j in range(0, len(r)):
            rgrid[i,j] = r[j]
            tgrid[i,j] = t[i]
            Tgrid[i,j] = m.T[t[i], r[j]].value

    fig = plt.figure(figsize=(10,6))
    ax = fig.add_subplot(1, 1, 1, projection='3d')
    ax.set_xlabel('Distance r')
    ax.set_ylabel('Time t')
    ax.set_zlabel('Temperature T')
    p = ax.plot_wireframe(rgrid, tgrid, Tgrid)<jupyter_output><empty_output><jupyter_text>Planar CoordinatesSuppressing the prime notation, for a slab geometry the model specializes to \begin{align*}\frac{\partial T}{\partial t} & = \frac{\partial^2 T}{\partial r^2}\end{align*}with auxiliary conditions\begin{align*}T(0, r) & = 0 & \forall 0 \leq r \leq 1 \\T(t, 1) & = 1 & \forall t > 0\\\frac{\partial T}{\partial r} (t, 0) & = 0 & \forall t \geq 0 \\\end{align*}<jupyter_code>!ls -a /content
from pyomo.environ import *
from pyomo.dae import *

m = ConcreteModel()

m.r = ContinuousSet(bounds=(0,1))
m.t = ContinuousSet(bounds=(0,2))

m.T = Var(m.t, m.r)

m.dTdt   = DerivativeVar(m.T, wrt=m.t)
m.dTdr   = DerivativeVar(m.T, wrt=m.r)
m.d2Tdr2 = DerivativeVar(m.T, wrt=(m.r, m.r))

@m.Constraint(m.t, m.r)
def pde(m, t, r):
    if t == 0:
        return Constraint.Skip
    if r == 0 or r == 1:
        return Constraint.Skip
    return m.dTdt[t,r] == m.d2Tdr2[t,r]

m.obj = Objective(expr=1)

m.ic  = Constraint(m.r, rule=lambda m, r:    m.T[0,r] == 0 if r > 0 and r < 1 else Constraint.Skip)
m.bc1 = Constraint(m.t, rule=lambda m, t:    m.T[t,1] == 1)
m.bc2 = Constraint(m.t, rule=lambda m, t: m.dTdr[t,0] == 0)

TransformationFactory('dae.finite_difference').apply_to(m, nfe=50, scheme='FORWARD', wrt=m.r)
TransformationFactory('dae.finite_difference').apply_to(m, nfe=50, scheme='FORWARD', wrt=m.t)
SolverFactory('ipopt', executable=ipopt_executable).solve(m, tee=True).write()
model_plot(m)<jupyter_output>Ipopt 3.12.8: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.8, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:    30347
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:    10299
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total[...]<jupyter_text>Cylindrical CoordinatesSuppressing the prime notation, for a cylindrical geometry the model specializes to \begin{align*}\frac{\partial T}{\partial t} & = \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)\end{align*}Expanding,\begin{align*}\frac{\partial T}{\partial t} & = \frac{\partial^2 T}{\partial t^2} + \frac{1}{r}\frac{\partial T}{\partial r}\end{align*}with auxiliary conditions\begin{align*}T(0, r) & = 0 & \forall 0 \leq r \leq 1\\T(t, 1) & = 1 & \forall t > 0\\\frac{\partial T}{\partial r} (t, 0) & = 0 & \forall t \geq 0 \\\end{align*}<jupyter_code>from pyomo.environ import *
from pyomo.dae import *

m = ConcreteModel()

m.r = ContinuousSet(bounds=(0,1))
m.t = ContinuousSet(bounds=(0,2))

m.T = Var(m.t, m.r)

m.dTdt   = DerivativeVar(m.T, wrt=m.t)
m.dTdr   = DerivativeVar(m.T, wrt=m.r)
m.d2Tdr2 = DerivativeVar(m.T, wrt=(m.r, m.r))

m.pde = Constraint(m.t, m.r, rule=lambda m, t, r: m.dTdt[t,r] == m.d2Tdr2[t,r] + (1/r)*m.dTdr[t,r]
        if r > 0 and r < 1 and t > 0 else Constraint.Skip)

m.ic  = Constraint(m.r, rule=lambda m, r:    m.T[0,r] == 0)
m.bc1 = Constraint(m.t, rule=lambda m, t:    m.T[t,1] == 1 if t > 0 else Constraint.Skip)
m.bc2 = Constraint(m.t, rule=lambda m, t: m.dTdr[t,0] == 0)

TransformationFactory('dae.finite_difference').apply_to(m, nfe=20, wrt=m.r, scheme='CENTRAL')
TransformationFactory('dae.finite_difference').apply_to(m, nfe=50, wrt=m.t, scheme='BACKWARD')
SolverFactory('ipopt', executable=ipopt_executable).solve(m).write()

model_plot(m)<jupyter_output># ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: -inf
  Upper bound: inf
  Number of objectives: 1
  Number of constraints: 4060
  Number of variables: 4110
  Sense: unknown
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Message: Ipopt 3.12.8\x3a Optimal Solution Found
  Termination condition: optimal
  Id: 0
  Error rc: 0
  Time: 0.16563725471496582
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0<jupyter_text>Spherical CoordinatesSuppressing the prime notation, for a cylindrical geometry the model specializes to \begin{align*}\frac{\partial T}{\partial t} & = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial T}{\partial r}\right)\end{align*}Expanding,\begin{align*}\frac{\partial T}{\partial t} & = \frac{\partial^2 T}{\partial t^2} + \frac{2}{r}\frac{\partial T}{\partial r}\end{align*}with auxiliary conditions\begin{align*}T(0, r) & = 0 & \forall 0 \leq r \leq 1\\T(t, 1) & = 1 & \forall t > 0\\\frac{\partial T}{\partial r} (t, 0) & = 0 & \forall t \geq 0 \\\end{align*}<jupyter_code>from pyomo.environ import *
from pyomo.dae import *

m = ConcreteModel()

m.r = ContinuousSet(bounds=(0,1))
m.t = ContinuousSet(bounds=(0,2))

m.T = Var(m.t, m.r)

m.dTdt   = DerivativeVar(m.T, wrt=m.t)
m.dTdr   = DerivativeVar(m.T, wrt=m.r)
m.d2Tdr2 = DerivativeVar(m.T, wrt=(m.r, m.r))

m.pde = Constraint(m.t, m.r, rule=lambda m, t, r: m.dTdt[t,r] == m.d2Tdr2[t,r] + (2/r)*m.dTdr[t,r]
        if r > 0 and r < 1 and t > 0 else Constraint.Skip)

m.ic  = Constraint(m.r, rule=lambda m, r:    m.T[0,r] == 0)
m.bc1 = Constraint(m.t, rule=lambda m, t:    m.T[t,1] == 1 if t > 0 else Constraint.Skip)
m.bc2 = Constraint(m.t, rule=lambda m, t: m.dTdr[t,0] == 0)

TransformationFactory('dae.finite_difference').apply_to(m, nfe=20, wrt=m.r, scheme='CENTRAL')
TransformationFactory('dae.finite_difference').apply_to(m, nfe=50, wrt=m.t, scheme='BACKWARD')
SolverFactory('ipopt', executable=ipopt_executable).solve(m).write()

model_plot(m)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.23064699205448355, 'fraction_numerical': 0.024290578887627697, 'mean_word_length': 3.171875, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 27, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '380114', 'license': 'Apache-2.0', 'n_tokens_mistral': 3303, 'n_tokens_neox': 3102, 'n_words': 794, 'path': 'notebooks/simulation/Heat_Conduction_in_Various_Geometries.ipynb', 'repo_name': 'edgBR/ND-Pyomo-Cookbook', 'size': '8810'}