1
1
#=
2
- # Dynamic Flow in simple Gas Network
2
+ # Dynamic Flow in Simple Gas Network
3
3
4
- This Example is based on the paper
4
+ This example is based on the paper
5
5
6
6
> Albertus J. Malan, Lukas Rausche, Felix Strehle, Sören Hohmann,
7
7
> Port-Hamiltonian Modelling for Analysis and Control of Gas Networks,
8
8
> IFAC-PapersOnLine, Volume 56, Issue 2, 2023, https://doi.org/10.1016/j.ifacol.2023.10.193.
9
9
10
- and tries replicate a simple simulation of flow in a 3-node gas network.
10
+ and tries to replicate a simple simulation of flow in a 3-node gas network.
11
11
12
12
This example can be dowloaded as a normal Julia script [here](@__NAME__.jl). #md
13
13
@@ -27,13 +27,13 @@ CairoMakie.activate!(type="svg") #hide
27
27
#=
28
28
## Node Models
29
29
30
- In this example, we use the equation based modeling using `ModelingToolkit.jl`.
31
- To verify the equations on a basic level we also provide units to eveything to
30
+ In this example, we use equation- based modeling using `ModelingToolkit.jl`.
31
+ To verify the equations on a basic level, we also provide units to everything to
32
32
perform dimensionality checks.
33
33
34
34
There are 2 node models used in the paper. The first node type has a constant
35
35
pressure.
36
- Additionally, we ad some "internal" state `q̃_inj` which we want to plot later (see also [Observables](@ref)).
36
+ Additionally, we add some "internal" state `q̃_inj` which we want to plot later (see also [Observables](@ref)).
37
37
=#
38
38
@mtkmodel ConstantPressureNode begin
39
39
@parameters begin
@@ -55,7 +55,7 @@ nothing #hide
55
55
The second node model is a variable pressure node. It has one output state (the pressure) and one input state,
56
56
the aggregated flows from the connected pipes.
57
57
As an internal state we have the injected flow from our source/load.
58
- The source/load behaviour itself is provided via a time dependent function.
58
+ The source/load behavior itself is provided via a time- dependent function.
59
59
=#
60
60
@mtkmodel VariablePressureNode begin
61
61
@structural_parameters begin
@@ -79,10 +79,10 @@ nothing #hide
79
79
#=
80
80
## Pipe Model
81
81
82
- The pipe is modeld as a first order ODE for the volumetric flow at the `dst` end. It has two inputs:
83
- the pressure at the source and and the pressure at the destination end.
82
+ The pipe is modeled as a first- order ODE for the volumetric flow at the `dst` end. It has two inputs:
83
+ the pressure at the source and the pressure at the destination end.
84
84
Later on, we'll specify the model to be antisymmetric, thus the flow is calculated explicitly for the
85
- destination end, but the source end will just recive just that times (-1).
85
+ destination end, but the source end will just receive that times (-1).
86
86
=#
87
87
@mtkmodel Pipe begin
88
88
@parameters begin
@@ -157,7 +157,6 @@ T = 278u"K" # simulation temperature
157
157
r = 0.012 u " mm" # pipe roughness
158
158
D = 0.6 u " m" # pipe diameter
159
159
160
- # # TODO : here is switched the lenths l12 and l13. The results are better. Is this a mistake in the paper?
161
160
L₁₂ = 90 u " km"
162
161
L₁₃ = 80 u " km"
163
162
L₂₃ = 100 u " km"
@@ -175,7 +174,7 @@ sinθ₂₃ = 0.0
175
174
nothing # hide
176
175
177
176
#=
178
- Lastly, we need to calculate the compressibility factor, the speed of sound and
177
+ Lastly, we need to calculate the compressibility factor, the speed of sound, and
179
178
the density at standard conditions:
180
179
=#
181
180
Z̃ = 1 - 3.52 * p̃/ pc * exp (- 2.26 * (T̃/ Tc)) + 0.274 * (p̃/ pc)^ 2 * exp (- 1.878 * (T̃/ Tc)) # (5)
@@ -188,16 +187,16 @@ nothing # hide
188
187
The equivalent "pressure capacity" at the nodes is calculated as a sum of the connected
189
188
pipe parameters according to (28).
190
189
191
- Here use defintions based on the speed and "standard" conditions.
190
+ Here we use definitions based on the speed and "standard" conditions.
192
191
=#
193
192
C₂ = L₁₂* A/ (2 * ρ̃* c̃^ 2 ) + L₂₃* A/ (2 * ρ̃* c̃^ 2 ) # (28)
194
193
C₃ = L₁₃* A/ (2 * ρ̃* c̃^ 2 ) + L₂₃* A/ (2 * ρ̃* c̃^ 2 ) # (28)
195
194
nothing # hide
196
195
197
196
#=
198
- Alternatively, we could calculate `Z2` and `Z3` based on the actuel pressure and simulation temperature.
199
- Then we could calculated the speed of sound for the "correct" conditions at the node.
200
- It seems to have very little effect on the actual results so I kept it simple.
197
+ Alternatively, we could calculate `Z2` and `Z3` based on the actual pressure and simulation temperature.
198
+ Then we could calculate the speed of sound for the "correct" conditions at the node.
199
+ It seems to have very little effect on the actual results, so I kept it simple.
201
200
=#
202
201
nothing # hide
203
202
@@ -212,12 +211,12 @@ from tracing it. To do so, we use `@register_symbolic` to declare it as a symbol
212
211
as a blackbox.
213
212
214
213
Additionally, we need to tell ModelingToolkit about the units of this object. This is just used
215
- for the static unit check during construction of the model. Later one , when we generate the
216
- Julia code from the symbolic reepresentation all units will be stripped.
214
+ for the static unit check during construction of the model. Later on , when we generate the
215
+ Julia code from the symbolic representation, all units will be stripped.
217
216
218
217
!!! note "Discontinuities in RHS"
219
- The picewise linear interpolated function creates discontinuities in the RHS of the system.
220
- However since we know the times exactly, we can handle this by simply giving a list of explicit
218
+ The piecewise linear interpolated function creates discontinuities in the RHS of the system.
219
+ However, since we know the times exactly, we can handle this by simply giving a list of explicit
221
220
tstops to the solve command, to make sure those are hit exactly.
222
221
=#
223
222
load2 = LinearInterpolation (- 1 * Float64[20 , 30 , 10 , 30 , 20 ], [0 , 4 , 12 , 20 , 24 ]* 3600.0 ; extrapolation= ExtrapolationType. Constant)
@@ -229,13 +228,13 @@ nothing #hide
229
228
#=
230
229
## Building the Network
231
230
232
- To bild the Network we first need to define the components. This is a two step process:
231
+ To build the network, we first need to define the components. This is a two- step process:
233
232
234
233
- first create the symbolic `ODESystem` using ModelingToolkit
235
234
- secondly build a NetworkDynamics component model ([`VertexModel`](@ref)/[`EdgeModel`](@ref)) based on the symbolic system.
236
235
237
236
In the first step we can use the keyword arguments to pass "default" values for our parameters and states.
238
- Those values will be automaticially transfered to the metadata of the component model the second step.
237
+ Those values will be automatically transferred to the metadata of the component model in the second step.
239
238
240
239
The second step requires to define the interface variables, i.e. what are the "input" states of your
241
240
component model and what are the "output" states.
@@ -253,10 +252,10 @@ v2 = VertexModel(v2_mtk, [:q̃_nw], [:p]; name=:v2, vidx=2)
253
252
v3 = VertexModel (v3_mtk, [:q̃_nw ], [:p ]; name= :v3 , vidx= 3 )
254
253
255
254
#=
256
- For the edge Model we have two inputs: the pressure on both source and destination end.
257
- There is a single output state: the volumetric flow. However we also need to tell
255
+ For the edge model we have two inputs: the pressure on both source and destination end.
256
+ There is a single output state: the volumetric flow. However, we also need to tell
258
257
NetworkDynamics about the coupling type. In this case we use `AntiSymmetric`, which
259
- meas that the source end will recieve the same flow, just inverted sign.
258
+ means that the source end will receive the same flow, just with inverted sign.
260
259
=#
261
260
@named e12_mtk = Pipe (; L= L₁₂, sinθ= sinθ₁₂, D, A, γ, η, r, g, T, Tc, pc, Rs, c̃, ρ̃, p̃)
262
261
@named e13_mtk = Pipe (; L= L₁₃, sinθ= sinθ₁₃, D, A, γ, η, r, g, T, Tc, pc, Rs, c̃, ρ̃, p̃)
@@ -269,10 +268,10 @@ e23 = EdgeModel(e23_mtk, [:p_src], [:p_dst], AntiSymmetric([:q̃]); name=:e23, s
269
268
#=
270
269
To build the network object we just need to pass the vertices and edges to the constructor.
271
270
272
- Note that we've used the `vidx` and `src`/`dst` keywords in the constructors to define for
273
- each component to which "part" of the network it belongs.
271
+ Note that we've used the `vidx` and `src`/`dst` keywords in the constructors to define
272
+ for each component to which "part" of the network it belongs.
274
273
275
- This means, the constructor can automaticially construct a graph based on those informations and
274
+ This means the constructor can automatically construct a graph based on that information and
276
275
we don't need to pass it explicitly.
277
276
=#
278
277
@@ -284,17 +283,17 @@ structurally different, because both functions capure a unique loadprofile funct
284
283
285
284
## Finding a Steady State
286
285
287
- To simulate the systme , we first need to find a steadystate . As a "guess" for that
286
+ To simulate the system , we first need to find a steady state . As a "guess" for that,
288
287
we create a `NWState` object from the network.
289
288
This will allocate flat arrays for states `u` and parameters `p` and fill them with the
290
289
default values.
291
290
=#
292
291
uguess = NWState (nw)
293
292
294
293
#=
295
- This is not a steadystate of the system however. To find a true steadystate we want to ensure
296
- that the lhs of the system is zero.
297
- We can solve for a steady state numerically by defining a Nonlinear Rootfind problem.
294
+ This is not a steady state of the system however. To find a true steady state, we want to ensure
295
+ that the LHS of the system is zero.
296
+ We can solve for a steady state numerically by defining a Nonlinear Rootfinding problem.
298
297
299
298
To do so, we need to wrap the Network object in a closure.
300
299
=#
@@ -307,15 +306,15 @@ initsol = solve(initprob)
307
306
308
307
#=
309
308
We can create a new `NWState` object by wrapping the solution from the nonlinear problem and the
310
- original prameters in a new `NWState` object.
309
+ original parameters in a new `NWState` object.
311
310
=#
312
311
u0 = NWState (nw, initsol. u, uguess. p)
313
312
314
313
#=
315
314
## Solving the ODE
316
315
317
316
Using this as our initial state we can create the actual `ODEProblem`.
318
- Since the ode allways operates on flat state and aprameter arrays we use `uflat` and `pflat` to extract them.
317
+ Since the ODE always operates on flat state and parameter arrays, we use `uflat` and `pflat` to extract them.
319
318
=#
320
319
prob = ODEProblem (nw, uflat (u0), (0.0 ,24 * 3600 ), copy (pflat (u0)))
321
320
sol = solve (prob, Tsit5 (), tstops= [0 ,4 ,12 ,20 ,24 ]* 3600 )
@@ -356,8 +355,8 @@ Notably, the "internal" states defined in the symbolic models are not "states" i
356
355
For example, we captured the load profile in the `q̃_inj` state of the `VariablePressureNode`.
357
356
The only dynamic state of the model however is `p`.
358
357
Using the "observables" mechanism from SciML, which is implemented by NetworkDynamics, we can reconstruct
359
- those "optimized" states which have been removed symbolicially .
360
- Here we plot the reconstructed load profile of nodes 2 and 3. Also, we know that node 1 is infinetly stiff,
358
+ those "optimized" states which have been removed symbolically .
359
+ Here we plot the reconstructed load profile of nodes 2 and 3. Also, we know that node 1 is infinitely stiff,
361
360
acting as an infinite source of volumetric flow. We can reconstruct this flow too.
362
361
=#
363
362
fig = begin
0 commit comments