"# Continuous Casting with Conic Supply Vessel\n",
"\n",
"In the following example, we consider a continuous steel casting production plant that is supplied with a conic vessel holding molten steel.\n",
"To ensure operations of the production plant, we require a througput of 6000 kg/min of molten steel and changing the vessels takes 60 seconds.\n",
"The height of the liquid steel inside the vessel is 2 metres initially and the radius at the bottom is $R_1 = 0.9$, the radius at the top is $R_2 = 1.2$ m.\n",
"The vessel has an outlet of radius $r$ and we need to optimise the radius such that we can maintain the required throughput and include the 60 second change-over time as well.\n",
"N.B. This is the same example as discussed in the lecture \"Numerische Modellierung in der Prozesstechnik\" and we discuss this example here to apply the methods of this course to the same example.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "JvESDWsUhe94"
},
"source": [
"## Torrielli's Law\n",
"\n",
"We model this using Torricelli's law, assuming\n",
"\n",
"- laminar and continuous flow\n",
"- ignorign the effects of viscosity\n",
"- incompressible liqued with uniform density\n",
"- $r << R_1, R_2$\n",
"- no effects from the geometry of the outlet orifice, etc.\n",
"\n",
"Then, Torricelli's law gives the velocity of the liquid as:\n",
"$$ v = \\sqrt{2gh}$$\n",
"where $g$ is the constant of gravity, and $h$ is the height of the liquid.\n",
"\n",
"The volumetric flow rate trough the hole at the bottom is given by\n",
"$$Q = av$$\n",
"with $a = \\pi r^2$.\n",
"Due to the conic shape of the vessel, the volume of a slice through the cone depends on the height, leading to:\n",
"$$ \\frac{dV}{dt} = - Q = - a \\sqrt{2 g h} = A(h) \\frac{dh}{dt}$$\n",
"\n",
"\n",
"\n",
"Hence, the differential equation we have to solve (for each $r$) is given by:\n",
"$$\n",
"\\frac{dh}{dt} = -\\frac{a}{\\pi \\left[ R_1 + s h \\right]^2} \\sqrt{2 g h}\n",
"$$\n",
"\n",
"where:\n",
"\n",
"$$\n",
"s = \\frac{R_2 - R_1}{H}\n",
"$$\n",
"\n",
"Substituting $ s $ into the equation:\n",
"\n",
"$$\n",
"\\frac{dh}{dt} = -\\frac{a}{\\pi \\left[ R_1 + \\left( \\dfrac{R_2 - R_1}{H} \\right) h \\right]^2} \\sqrt{2 g h}\n",
"# constant of gravity, take the pre-defined one instead of defining our own\n",
"import scipy\n",
"\n",
"\n",
"import torch\n",
"from torchdiffeq import odeint\n",
"\n",
"\n",
"from datetime import datetime"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "EiNuMdQfkauk"
},
"source": [
"# Plot Height vs Time\n",
"\n",
"As a first step, we visualise how the height of the molten steel in the conic vessel changes with time.\n",
"\n",
"We will use $r=2$ cm for this example, though the radius $r$ is the quantity we want to optimise later to make sure we can maintain the required throughput, including the change-over time.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "AcIz5EjClHSh"
},
"outputs": [],
"source": [
"# Geometry of the truncated cone\n",
"H = 2.0 # Total height of the vessel (m)\n",
"R1 = 0.9 # Radius at the bottom (m)\n",
"R2 = 1.2 # Radius at the top (m)\n",
"s = (R2 - R1) / H # Slope of the radius as a function of height\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "ccZ5IeSUwf4j"
},
"outputs": [],
"source": [
"# assume a fixed radius of the orifice\n",
"a = np.pi * (0.02)**2 # Area of the hole at the bottom (m^2)\n",
" # assuming r=2cm"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "OP5j1PjUlMN8"
},
"outputs": [],
"source": [
"# Time parameters\n",
"t0 = 0.0 # Initial time (s)\n",
"t_end = 3000.0 # End time (s)\n",
"dt = 0.1 # Time step (s)\n",
"t_values = torch.arange(t0, t_end + dt, dt)\n",
"\n",
"# Initial condition\n",
"h0 = torch.tensor([H], dtype=torch.float32) # Initial height of the liquid (full vessel)\n",
"\n",
"# Physical Constants\n",
"g = torch.tensor(# YOUR CODE HERE, dtype=torch.float32)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ATNqzi02nMgB"
},
"source": [
"This is the differential equation from above, assuming Torricelli's law.\n",
"\n",
"We need to encode this as ```torch``` tensors, so that we can make use of the advanced compute capabilities of PyTorch.\n",
"\n",
"We use [torch.clamp](https://pytorch.org/docs/stable/generated/torch.clamp.html) to make sure that the height of the steel in the vessel does not become negative.\n",
"\n",
"**N.B.** The function ```def dhdt(t,h)``` also depends on the parameter ```a``` (the area of the outlet orifice) and ```g```. Here, we have specified these as a global variables above.\n",
"\n",
"Since the function ```odeint``` does not take additional arguments, we have to specify these outside the scope of the function.\n",
"This is not an issue at this stage, but we will see later that this requires us to write the code in a specific way when we want to optimise the radius of the orifice and, hence, the area ```a```."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "RRO687otlxU5"
},
"outputs": [],
"source": [
"# Function defining the differential equation dh/dt\n",
"We will now find the optimal radius of the output orifice $r$.\n",
"As discussed above, we need to maintain a specific throughput, as well as consider the time it takes to change vessels for continuous operations.\n",
"\n",
"This means that we need to find a radius of the output orifice such that the vessel is empty at the target time ```t_target = 428.32``` s.\n",
"\n",
"N.B. in the following, we will denote the radius of the output orifice for each step of the optimisation as $r_0$ to avoid confusion with values inside each optimisation step.\n",
"We now define the function that computes the height of the remaining molten steel inside the vessel at the target time ```t_target``` for a given radius $r_0$ of the output orifice.\n",
"\n",
"We then want to find the value $r_0$ where the height is zero at the end, i.e. when the vessel is empty.\n",
"\n",
"\n",
"This is implemented in the following way: For each radius, we solve the differential equation, calculate the height at various times $t$, and, in particular the time at ```t_target```.\n",
"\n",
"Since we want to find the value of $r_0$ for which the function is zero, we want to find the *root* of the function of the heights as we change $r_0$.\n",
"\n",
"Note:\n",
"Unlike earlier, the area of the orifice now depends on the radius $r_0$ that we are changing. Ideally, we would add another argument to the function with the differential equation.\n",
"However, the function ```odeint``` of [torchdiffeq](https://github.com/rtqichen/torchdiffeq) does not support this. Therefore, we define the function for the differential equation within the function computing the height at time ```t_target```. Then, for each round, the value of the area $a$ is defined and we, kind-off, side-step the problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "GJG0nbEtsSof"
},
"outputs": [],
"source": [
"# Compute the height of the molten steel for a given\n",
"# radius of the outlet orifice r_0 after t_target has elapsed\n",
"\n",
"def compute_h_at_tend(r_o):\n",
" # Ensure r_o is a tensor with requires_grad=True\n",
" h_tend = h_values[-1, 0] # Height at t_target\n",
"\n",
" return h_tend, r_o"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "bj2N5KIL0CB9"
},
"source": [
"Before we find the optimal value, we can plot the result to get a feeling for our expectation.\n",
"\n",
"N.B. due to the line ```h = h.clamp(min=0.0)``` we cannot have negative heights in the vessel."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "1dEYfP_s0jWq"
},
"source": [
"**Exercise**\n",
"\n",
"Plot the function of the height of the remaining steel in the vessel after the target time ```t_target``` has elapsed for the radius in the interval $(0.02, 0.05)$ for 100 steps.\n",
"Note that we need PyTorch tensors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "V7OcDeyg0umW"
},
"outputs": [],
"source": [
"# Generate x_space and compute y_space\n",
"x_space = torch.linspace(0.02, 0.05, steps=50) # Orifice radii from 0.02 m to 0.05 m\n",
"Instead of using a package to solve the differential equation, we can also use an iterative approach and discretise the differential equation.\n",
"The following function computes the height at a specific target time ```desired_time``` as a function of the radius ```r``` of the output orifice.\n",
"\n",
"The remainnig parameters (```R1, R2, initial_height```) specifiy the geometry and initial height of the molten steel in the vessel. We can make these parameters argument to the function here as we are not constrained by the calling function like ```odeint``` earlier.\n",
"\n",
"N.B. We need to work with PyTorch tensors (or types that are easily converted) to be able to call Newton's method efficiently later on."
" # Update current_height for the next iteration\n",
" current_height = next_height\n",
"\n",
" # The height at the desired time\n",
" height_at_desired_time = current_height\n",
"\n",
" return height_at_desired_time"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "kkMJ_cT9E-H4"
},
"source": [
"**Exercise**\n",
"\n",
"Use the above function to call the function ```newtons_method``` we defined earlier and find the optimal radius of the output orifice."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "ybu1ouSSFQhG"
},
"outputs": [],
"source": [
"##\n",
"## your code goes here\n",
"##"
]
}
],
"metadata": {
"colab": {
"provenance": []
},
"kernelspec": {
"display_name": "Python 3",
"name": "python3"
},
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
%% Cell type:markdown id: tags:
# Continuous Casting with Conic Supply Vessel
In the following example, we consider a continuous steel casting production plant that is supplied with a conic vessel holding molten steel.
To ensure operations of the production plant, we require a througput of 6000 kg/min of molten steel and changing the vessels takes 60 seconds.
The height of the liquid steel inside the vessel is 2 metres initially and the radius at the bottom is $R_1 = 0.9$, the radius at the top is $R_2 = 1.2$ m.
The vessel has an outlet of radius $r$ and we need to optimise the radius such that we can maintain the required throughput and include the 60 second change-over time as well.
N.B. This is the same example as discussed in the lecture "Numerische Modellierung in der Prozesstechnik" and we discuss this example here to apply the methods of this course to the same example.
%% Cell type:markdown id: tags:
## Torrielli's Law
We model this using Torricelli's law, assuming
- laminar and continuous flow
- ignorign the effects of viscosity
- incompressible liqued with uniform density
- $r << R_1, R_2$
- no effects from the geometry of the outlet orifice, etc.
Then, Torricelli's law gives the velocity of the liquid as:
$$ v = \sqrt{2gh}$$
where $g$ is the constant of gravity, and $h$ is the height of the liquid.
The volumetric flow rate trough the hole at the bottom is given by
$$Q = av$$
with $a = \pi r^2$.
Due to the conic shape of the vessel, the volume of a slice through the cone depends on the height, leading to:
$$ \frac{dV}{dt} = - Q = - a \sqrt{2 g h} = A(h) \frac{dh}{dt}$$
Hence, the differential equation we have to solve (for each $r$) is given by:
$$
\frac{dh}{dt} = -\frac{a}{\pi \left[ R_1 + s h \right]^2} \sqrt{2 g h}
$$
where:
$$
s = \frac{R_2 - R_1}{H}
$$
Substituting $ s $ into the equation:
$$
\frac{dh}{dt} = -\frac{a}{\pi \left[ R_1 + \left( \dfrac{R_2 - R_1}{H} \right) h \right]^2} \sqrt{2 g h}
$$
%% Cell type:code id: tags:
```
# make sure we have the required libraries
# ! pip install torchdiffeq
```
%% Cell type:code id: tags:
```
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# constant of gravity, take the pre-defined one instead of defining our own
import scipy
import torch
from torchdiffeq import odeint
from datetime import datetime
```
%% Cell type:markdown id: tags:
# Plot Height vs Time
As a first step, we visualise how the height of the molten steel in the conic vessel changes with time.
We will use $r=2$ cm for this example, though the radius $r$ is the quantity we want to optimise later to make sure we can maintain the required throughput, including the change-over time.
%% Cell type:code id: tags:
```
# Geometry of the truncated cone
H = 2.0 # Total height of the vessel (m)
R1 = 0.9 # Radius at the bottom (m)
R2 = 1.2 # Radius at the top (m)
s = (R2 - R1) / H # Slope of the radius as a function of height
```
%% Cell type:code id: tags:
```
# assume a fixed radius of the orifice
a = np.pi * (0.02)**2 # Area of the hole at the bottom (m^2)
# assuming r=2cm
```
%% Cell type:code id: tags:
```
# Time parameters
t0 = 0.0 # Initial time (s)
t_end = 3000.0 # End time (s)
dt = 0.1 # Time step (s)
t_values = torch.arange(t0, t_end + dt, dt)
# Initial condition
h0 = torch.tensor([H], dtype=torch.float32) # Initial height of the liquid (full vessel)
# Physical Constants
g = torch.tensor(# YOUR CODE HERE, dtype=torch.float32)
```
%% Cell type:markdown id: tags:
This is the differential equation from above, assuming Torricelli's law.
We need to encode this as ```torch``` tensors, so that we can make use of the advanced compute capabilities of PyTorch.
We use [torch.clamp](https://pytorch.org/docs/stable/generated/torch.clamp.html) to make sure that the height of the steel in the vessel does not become negative.
**N.B.** The function ```def dhdt(t,h)``` also depends on the parameter ```a``` (the area of the outlet orifice) and ```g```. Here, we have specified these as a global variables above.
Since the function ```odeint``` does not take additional arguments, we have to specify these outside the scope of the function.
This is not an issue at this stage, but we will see later that this requires us to write the code in a specific way when we want to optimise the radius of the orifice and, hence, the area ```a```.
%% Cell type:code id: tags:
```
# Function defining the differential equation dh/dt
def dhdt(t, h):
# Ensure h is non-negative
h = # YOUR CODE HERE
# Compute r(h)
r = R1 + s * h
# Compute A(h)
A = torch.pi * r ** 2
# Compute dh/dt
sqrt_term = torch.sqrt(2 * g * h)
dh = - (a / A) * sqrt_term
# Handle h = 0 to prevent division by zero
dh = torch.where(h > 0.0, dh, torch.tensor(0.0))
return dh
```
%% Cell type:code id: tags:
```
# Solve the ODE using torchdiffeq's odeint function
We will now find the optimal radius of the output orifice $r$.
As discussed above, we need to maintain a specific throughput, as well as consider the time it takes to change vessels for continuous operations.
This means that we need to find a radius of the output orifice such that the vessel is empty at the target time ```t_target = 428.32``` s.
N.B. in the following, we will denote the radius of the output orifice for each step of the optimisation as $r_0$ to avoid confusion with values inside each optimisation step.
We now define the function that computes the height of the remaining molten steel inside the vessel at the target time ```t_target``` for a given radius $r_0$ of the output orifice.
We then want to find the value $r_0$ where the height is zero at the end, i.e. when the vessel is empty.
This is implemented in the following way: For each radius, we solve the differential equation, calculate the height at various times $t$, and, in particular the time at ```t_target```.
Since we want to find the value of $r_0$ for which the function is zero, we want to find the *root* of the function of the heights as we change $r_0$.
Note:
Unlike earlier, the area of the orifice now depends on the radius $r_0$ that we are changing. Ideally, we would add another argument to the function with the differential equation.
However, the function ```odeint``` of [torchdiffeq](https://github.com/rtqichen/torchdiffeq) does not support this. Therefore, we define the function for the differential equation within the function computing the height at time ```t_target```. Then, for each round, the value of the area $a$ is defined and we, kind-off, side-step the problem.
%% Cell type:code id: tags:
```
# Compute the height of the molten steel for a given
# radius of the outlet orifice r_0 after t_target has elapsed
Before we find the optimal value, we can plot the result to get a feeling for our expectation.
N.B. due to the line ```h = h.clamp(min=0.0)``` we cannot have negative heights in the vessel.
%% Cell type:markdown id: tags:
**Exercise**
Plot the function of the height of the remaining steel in the vessel after the target time ```t_target``` has elapsed for the radius in the interval $(0.02, 0.05)$ for 100 steps.
Note that we need PyTorch tensors.
%% Cell type:code id: tags:
```
# Generate x_space and compute y_space
x_space = torch.linspace(0.02, 0.05, steps=50) # Orifice radii from 0.02 m to 0.05 m
plt.title('Height of Liquid in Truncated Cone Over Time (Optimal r_o via Newton\'s Method)')
plt.xlabel('Time (s)')
plt.ylabel('Height (m)')
plt.grid(True)
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
# Optimising Outlet Radius - Iterative Approach
Instead of using a package to solve the differential equation, we can also use an iterative approach and discretise the differential equation.
The following function computes the height at a specific target time ```desired_time``` as a function of the radius ```r``` of the output orifice.
The remainnig parameters (```R1, R2, initial_height```) specifiy the geometry and initial height of the molten steel in the vessel. We can make these parameters argument to the function here as we are not constrained by the calling function like ```odeint``` earlier.
N.B. We need to work with PyTorch tensors (or types that are easily converted) to be able to call Newton's method efficiently later on.
