diff --git a/Oven Brick Model.ipynb b/Oven Brick Model.ipynb
index d09fa32..e7152b1 100644
--- a/Oven Brick Model.ipynb
+++ b/Oven Brick Model.ipynb
@@ -129,6 +129,309 @@
"print(solution.gradient(bc=[{\"derivative\": 1}, {\"value\": 0}]).data)"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "57d2748e",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\partial_t u = - \\frac{1}{2} \\left| \\nabla u \\right|^2 - \\nabla^2 u - \\nabla^4 u \n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 68,
+ "id": "87ac9698",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from pde import PDEBase, ScalarField, UnitGrid\n",
+ "\n",
+ "\n",
+ "class KuramotoSivashinskyPDE(pde.PDEBase):\n",
+ " \"\"\"Implementation of the normalized Kuramoto–Sivashinsky equation\"\"\"\n",
+ "\n",
+ " def evolution_rate(self, state, t=0):\n",
+ " \"\"\"implement the python version of the evolution equation\"\"\"\n",
+ " state_lap = state.laplace(bc=\"auto_periodic_neumann\")\n",
+ " state_lap2 = state_lap.laplace(bc=\"auto_periodic_neumann\")\n",
+ " state_grad = state.gradient(bc=\"auto_periodic_neumann\")\n",
+ " return -state_grad.to_scalar(\"squared_sum\") / 2 - state_lap - state_lap2\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c20ac3d8",
+ "metadata": {},
+ "source": [
+ "# Heat diffusion equation with Coke Oven Brick Properties"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fcd17e59",
+ "metadata": {},
+ "source": [
+ "## For constant $k$, $C_p$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e87a0bbf",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\partial_t T = \\alpha \\Delta T\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "116d5038",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\partial_t T = \\frac{k}{\\rho C_p} \\Delta T\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dc87d038",
+ "metadata": {},
+ "source": [
+ "## For varying $k$, $C_p$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d8d70c8f",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\partial_t T = \\frac{1}{\\rho C_p} \\nabla \\cdot \\left( k \\nabla T \\right)\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f9ed2160",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\partial_t T = \\frac{1}{\\rho C_p} \\left(\\nabla k \\cdot \\nabla T + k \\nabla^2 T \\right)\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "511804d4",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "k = c_1 T + c_0\n",
+ "$$\n",
+ "\n",
+ "$$\n",
+ "\\nabla k = c_1 \\nabla T\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ab5e4123",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "k = c_1 T + c_0\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 97,
+ "id": "5a43dfeb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "class CokeOvenBrickHeatEqn(pde.PDEBase):\n",
+ " \"\"\"Implementation of the normalized Kuramoto–Sivashinsky equation\"\"\"\n",
+ "\n",
+ " def __init__ (self, bc=\"auto_periodic_neumann\"):\n",
+ " self.bc = bc\n",
+ " self.rho = 1900 # kg / m3\n",
+ " self.kCoef0 = 0.93 # W / m / K\n",
+ " self.kCoef1 = 0.698e-3 # W / m / K2\n",
+ " self.cpCoef0 = 837.2 # J / kg / K\n",
+ " self.cpCoef1 = 251.2e-3 # J / kg / K2\n",
+ "\n",
+ " def k (self, T):\n",
+ " return T * kCoef1 + kCoef0 \n",
+ " \n",
+ " def cp (self, T):\n",
+ " return T * cpCoef1 + cpCoef0\n",
+ " \n",
+ " def update_bc (self, gradT_chamber=None, T_oven=None):\n",
+ " bc0, bc1 = self.bc\n",
+ " if gradT_chamber:\n",
+ " self.bc[0] = {\"derivative\": gradT_chamber} \n",
+ " if T_oven:\n",
+ " self.bc[1] = {\"value\": T_oven}\n",
+ "\n",
+ " def evolution_rate(self, state, t=0):\n",
+ " \"\"\"implement the python version of the evolution equation\"\"\"\n",
+ " state_lap = state.laplace(bc=self.bc)\n",
+ " state_grad = state.gradient(bc=self.bc)\n",
+ " \n",
+ " k = self.kCoef1 * state + self.kCoef0\n",
+ " cp = self.cpCoef1 * state + self.cpCoef0\n",
+ " \n",
+ " k_grad = self.kCoef1 * state_grad\n",
+ " \n",
+ " return (state_grad.dot(k_grad) + k * state_lap) / self.rho / cp\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 98,
+ "id": "8f2a53b4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "T0 = 800 + 273.15\n",
+ "simulation_time = 24 * 60 * 60 # 24 hours in seconds\n",
+ "dt = 1 # timestep 1 second\n",
+ "bc_update_period = 60 # update bc every 60 seconds\n",
+ "brick_thickness = 0.14 # m\n",
+ "n_grid_brick = 32\n",
+ "\n",
+ "grid = pde.CartesianGrid([[0, brick_thickness]], n_grid_brick, periodic=False)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 91,
+ "id": "d28a0363",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "3761b63479294b0a83c952d90335ba4a",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ " 0%| | 0/1.0 [00:00\n",
+ " \n",
+ " Figure\n",
+ "
\n",
+ " 
\n",
+ " \n",
+ " "
+ ],
+ "text/plain": [
+ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[-4.93749698e+01 -9.66271511e+01 -9.43865529e+01 -9.21586376e+01\n",
+ " -8.99488392e+01 -8.77625494e+01 -8.56051044e+01 -8.34817711e+01\n",
+ " -8.13977342e+01 -7.93580833e+01 -7.73678001e+01 -7.54317463e+01\n",
+ " -7.35546517e+01 -7.17411020e+01 -6.99955284e+01 -6.83221965e+01\n",
+ " -6.67251959e+01 -6.52084306e+01 -6.37756098e+01 -6.24302389e+01\n",
+ " -6.11756113e+01 -6.00148009e+01 -5.89506549e+01 -5.79857873e+01\n",
+ " -5.71225729e+01 -5.63631422e+01 -5.57093761e+01 -5.51629023e+01\n",
+ " -5.47250916e+01 -5.43970549e+01 -5.41796406e+01 -2.45345502e+05]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "solution = eq.solve(solution, t_range=10000, dt=0.01)\n",
+ "solution.plot()\n",
+ "plt.ylim(1000, 1100)\n",
+ "print(solution.gradient(bc=[{\"derivative\": 1}, {\"value\": 0}]).data)"
+ ]
+ },
{
"cell_type": "markdown",
"id": "3164bd16",