{ "cells": [ { "cell_type": "markdown", "id": "4bc9f590-888f-42b9-81f7-1096612adde8", "metadata": {}, "source": [ "## \"Old-school\" PyFVTool: the heat equation in a 1D Cartesian system\n", "\n", "This particular Notebook demonstrates the \"old-school\" expert approach to setting up the finite-volume solution to the heat equation using PyFVTool. The discretized terms of the matrix equation are explicitly evaluated and cast into the final sparse matrix equation which is then solved numerically with `pf.solveMatrixPDE()`.\n", "\n", "This example illustrates the inner workings of the finite-volume method. It is functionally and numerically identical to `cartesian1D_heat_equation_neumann.ipynb`. We compare the FVM solution by PyFVTool to the analytic solutions." ] }, { "cell_type": "code", "execution_count": 1, "id": "106afa46-7bfb-4a1b-9fac-4ffb69e45277", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.special import erf\n", "import matplotlib.pyplot as plt\n", "import pyfvtool as pf" ] }, { "cell_type": "code", "execution_count": 2, "id": "ad3537ee-83db-4f81-b7f2-48638d6fd8bb", "metadata": {}, "outputs": [], "source": [ "# Spoiler alert!\n", "\n", "def T_analytic_dirichlet(x,t, alfa, T0, Ts):\n", " return (T0-Ts)*erf(x/np.sqrt(4*alfa*t))+Ts\n", "\n", "def T_analytic_neumann(x,t, alfa, T0, k, qs):\n", " return T0 + qs/k*np.sqrt(4*alfa*t/np.pi)*np.exp(-x**2/(4*alfa*t))\\\n", " - qs/k*x*(1-erf(x/np.sqrt(4*alfa*t)))" ] }, { "cell_type": "code", "execution_count": 3, "id": "df5f937e-551d-47b7-8604-0739f83b5758", "metadata": {}, "outputs": [], "source": [ "# Parameters\n", "L = 1.0 # [m] domain length\n", "k = 20.0 # 0.6 for water, 0.025 for air W/m/K\n", "rho = 8000.0 # kg/m^3\n", "c = 500.0 # J/kg/K (4200 for water, 1000 for air)\n", "alfa = k/(rho*c) # heat diffusion\n", "T0 = 300.0 # [K]\n", "Ts = 350.0 # [K]\n", "qs = 1000 # [W/m^2]\n", "t_sim = L**2/(20*alfa) # [s]\n", "time_steps = 50\n", "dt = t_sim/time_steps # \n", "Nx = 50 # number of cells" ] }, { "cell_type": "code", "execution_count": 4, "id": "ae1d51bc-4e13-43e9-88e9-3a7479f21610", "metadata": {}, "outputs": [], "source": [ "m = pf.Grid1D(Nx, L)" ] }, { "cell_type": "code", "execution_count": 5, "id": "cc6c5f88-467b-4f95-bcce-2fdbee91da4f", "metadata": {}, "outputs": [], "source": [ "# Boundary condition\n", "left_bc = \"Neumann\"\n", "\n", "BC = pf.BoundaryConditions(m)\n", "if left_bc == \"Dirichlet\":\n", " BC.left.a[:] = 0.0\n", " BC.left.b[:] = 1.0\n", " BC.left.c[:] = Ts\n", " T_analytic = lambda x,t: T_analytic_dirichlet(x, t, alfa, T0, Ts)\n", "else:\n", " BC.left.a[:] = k\n", " BC.left.b[:] = 0.0\n", " BC.left.c[:] = -qs\n", " T_analytic = lambda x,t: T_analytic_neumann(x, t, alfa, T0, k, qs)" ] }, { "cell_type": "code", "execution_count": 6, "id": "f8f33dcd-ac71-4526-90d8-46eeeacde51e", "metadata": {}, "outputs": [], "source": [ "# Initial condition\n", "T_init = pf.CellVariable(m, T0, BC) # initial condition" ] }, { "cell_type": "code", "execution_count": 7, "id": "3d2335bc-6b3a-4098-a3b5-2196af34c23c", "metadata": {}, "outputs": [], "source": [ "# physical parameters\n", "alfa_cell = pf.CellVariable(m, alfa, pf.BoundaryConditions(m))\n", "alfa_face = pf.harmonicMean(alfa_cell)" ] }, { "cell_type": "code", "execution_count": 8, "id": "f3bf6b40-3aa2-48a3-a2e9-5c0d64d30f52", "metadata": {}, "outputs": [], "source": [ "M_diff = pf.diffusionTerm(alfa_face)\n", "[M_bc, RHS_bc] = pf.boundaryConditionsTerm(BC)" ] }, { "cell_type": "code", "execution_count": 9, "id": "67534101-60dc-42f4-8624-2eeb40b1cc6d", "metadata": {}, "outputs": [], "source": [ "t=0\n", "while t" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(1)\n", "plt.clf()\n", "plt.title('Left BC: '+left_bc)\n", "plt.plot(x, T_an, x, T_num, 'o')\n", "plt.legend({'Analytic', 'PyFVTool FVM'})\n", "plt.xlabel('x [m]')\n", "plt.ylabel('T [K]')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 12, "id": "a56c0184-6fe6-4ae2-810b-a8fbe39f6e78", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.1994295875793644e-05\n" ] } ], "source": [ "er = np.sum(np.abs(T_num-T_an)/T_an)/Nx\n", "print(er)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.9" } }, "nbformat": 4, "nbformat_minor": 5 }