Arjen Markus (3 march 2006) I have been wanting to bring in some physics into the Wiki for some time now. The problem is that solving physical puzzles always carries a lot of context with it. The script below shows how to solve one particular puzzle that is not too difficult.
The actual problem is artistically based on things I deal with for a living: it considers the spreading of some pollutant in a river system. The river is a slowly flowing body of water, 100 km long. It has a main channel and a shallow bank along it (which is almost stagnant water). The pollutant is spilled at time zero (so I can use the initial condition to simulate it). Then the pollutant is followed over ten days and the concentration at three points along the river is monitored.
There are plenty of things to do before this becomes something even remotely useful:
Still, it was fun to implement and test :)
# wqmodel.tcl --
# Initial set-up for a toy water quality model:
# - Define a schematisation with finite volumes
# - Define the substances to be modelled
# - Solve the advection-diffusion equations
#
# Note:
# Using arrays may not be very efficient, it
# does keep the code simple for the moment
#
namespace eval ::WQModel {
variable segments {}
variable substances {}
variable volume
variable surface
variable exchanges
variable conc
variable mass
variable deriv
namespace export segment substances initialise timeframe exchange \
initial-values boundary-values getconc nextstep
}
# segment --
# Define a segment by its ID and volume and surface area
#
# Arguments:
# id (Numerical) ID of the segment
# vol Volume of the segment
# surf Surface area (to compute the depth for instance)
#
# Result:
# None
#
proc ::WQModel::segment {id vol surf} {
variable segments
variable volume
variable surface
lappend segments $id
set volume($id) $vol
set surface($id) $surf
}
# substances --
# Define the names (and order) of the substances
#
# Arguments:
# args List of substance names
#
# Result:
# None
#
proc ::WQModel::substances {args} {
variable substances
set substances $args
}
# initial-values --
# Define the initial values of the substances per segment
#
# Arguments:
# id Segment ID
# values List of initial values
#
# Result:
# None
#
proc ::WQModel::initial-values {id values} {
variable conc
set conc($id) $values
}
# boundary-values --
# Define the boundary values of the substances per
# boundary segment
#
# Arguments:
# id Segment ID (starting with "B")
# values List of boundary values
#
# Result:
# None
#
proc ::WQModel::boundary-values {id values} {
variable segments
variable conc
variable deriv
set conc($id) $values
set deriv($id) $values ;# They must simply exist
}
# exchange --
# Define the exchanges between segments
#
# Arguments:
# from ID of segment "from"
# to ID of segment "to"
# flow Flow rate
# disp Dispersion coefficient
# area Area of exchange surface
# distance Representative distance between centres
#
# Result:
# None
#
# Note:
# Boundary segments should have an ID starting
# with "B"
#
proc ::WQModel::exchange {from to flow disp area distance} {
variable exchanges
lappend exchanges $from $to $flow $disp $area $distance
}
# timeframe --
# Define the start and stop time and time step
#
# Arguments:
# start Start time
# stop Stop time
# timestep Time step
#
# Result:
# None
#
# Note:
# Time in days (flow rate in m3/s)
#
proc ::WQModel::timeframe {start stop timestep} {
variable start_time
variable stop_time
variable time_step
set start_time $start
set stop_time $stop
set time_step $timestep
}
# initialise_derivs --
# Initialise the derivatives
#
# Arguments:
# None
#
# Result:
# None
#
proc ::WQModel::initialise_derivs {} {
variable segments
variable substances
variable deriv
set dv {}
foreach s $substances {
lappend dv 0.0
}
foreach id $segments {
set deriv($id) $dv
}
}
# initialise --
# Initialise the computation
#
# Arguments:
# None
#
# Result:
# None
#
proc ::WQModel::initialise {} {
variable segments
variable substances
variable conc
variable volume
variable mass
variable time
variable start_time
set time $start_time
foreach id $segments {
set mass($id) {}
set vol $volume($id)
foreach s $substances c $conc($id) {
lappend mass($id) [expr {$c*$vol}]
}
}
}
# new_concentrations --
# Compute the new concentrations
#
# Arguments:
# None
#
# Result:
# None
#
proc ::WQModel::new_concentrations {} {
variable segments
variable substances
variable conc
variable mass
variable volume
variable deriv
variable time_step
set dv {}
foreach s $substances {
lappend dv 0.0
}
foreach id $segments {
set newconc {}
set newmass {}
foreach m $mass($id) dv $deriv($id) v $volume($id) {
set m [expr {$m + $dv * $time_step}]
lappend newmass $m
lappend newconc [expr {$m / $v}]
}
set conc($id) $newconc
set mass($id) $newmass
}
}
# transport_processes --
# Compute the derivatives due to transport
#
# Arguments:
# None
#
# Result:
# None
#
proc ::WQModel::transport_processes {} {
variable segments
variable conc
variable deriv
variable exchanges
foreach {from to flow disp area distance} $exchanges {
#
# Straightforward backward differences ... not accurate, but simple!
#
if { $flow > 0.0 } {
set qfrom [expr {86400.0*$flow}]
set qto 0.0
} else {
set qfrom 0.0
set qto [expr {-86400.0*$flow}]
}
set disprate [expr {86400.0*$disp*$area/$distance}]
set newfrom {}
set newto {}
foreach cfrom $conc($from) dvfrom $deriv($from) \
cto $conc($to) dvto $deriv($to) {
lappend newfrom [expr {$dvfrom - $cfrom * $qfrom + $cto *$qto - $disprate * ($cfrom-$cto)}]
lappend newto [expr {$dvto + $cfrom * $qfrom - $cto *$qto - $disprate * ($cto-$cfrom)}]
}
set deriv($from) $newfrom
set deriv($to) $newto
}
}
# nextstep --
# Compute the concentrations for the next step
#
# Arguments:
# None
#
# Result:
# 1 if there is a next step, 0 if the computation has ended
#
proc ::WQModel::nextstep {} {
variable time
variable time_step
variable stop_time
#
# Prepare the step
#
initialise_derivs
#
# Deal with water quality processes
#
# TODO
#
# Add the waste loads
#
# TODO
#
# Transport processes
#
transport_processes
#
# Compute the new concentrations
#
new_concentrations
set time [expr {$time + $time_step}]
set ::time $time
return [expr {$time <= $stop_time}]
}
# getconc --
# Return the concentrations for a particular segment
#
# Arguments:
# id The ID of the segment
#
# Result:
# List of concentrations
#
proc ::WQModel::getconc {id} {
variable conc
return $conc($id)
}
# main --
# A stretch of river with a calamitous spill
# For more realism: shallow sides with almost stagnant water
#
namespace import ::WQModel::*
substances pesticide
#
# The segments and their volumes
#
set length 1000.0 ;# segment is 1 kilometer long
set width 200.0 ;# 200 m wide
set depth 5.0 ;# 5 m deep
for { set i 0 } { $i < 200 } { incr i 2 } {
set j [expr {$i+1}]
segment $i [expr {$length*$width*$depth}] [expr {$length*$width}]
segment $j [expr {0.1*$length*$width*$depth}] [expr {0.3*$length*$width}]
initial-values $i 0.0
initial-values $j 0.0
}
#
# The exchanges between main channel and shallow sides
#
set area [expr {$length*0.1/0.3}]
set distance [expr {0.5*($width+0.3*$width)}]
set flow 0.0 ;# m3/s
set disp 0.001 ;# m2/s
for { set i 0 } { $i < 200 } { incr i 2 } {
set j [expr {$i+1}]
exchange $i $j $flow $disp $area $distance
}
#
# The exchanges along main channel
#
set area [expr {$width*$depth}]
set distance $length
set flow 100.0 ;# m3/s
set disp 0.001 ;# m2/s
for { set i -2 } { $i < 200 } { incr i 2 } {
set j [expr {$i+2}]
set from $i
set to $j
if { $i == -2 } {
set from "B1"
}
if { $j >= 200 } {
set to "B2"
}
exchange $from $to $flow $disp $area $distance
}
#
# Define the spill
#
boundary-values B1 0.0
boundary-values B2 0.0
initial-values 10 1.0
#
# Now compute the solution
#
global time
timeframe 0.0 10.0 0.01
initialise
while { [nextstep] } {
if { abs(10.0*$time-int(10.0*$time+0.0001)) < 0.0001 } {
puts "$time [getconc 40] [getconc 80] [getconc 160]"
}
}