Spring2008:Bradshaw
From GeoMod
Contents |
[edit]
01/28/2008 Arrays
from visual import *
from random import uniform,seed
scene.forward = (0, -0.5,1)
scene.up = (0,0,1)
scene.forward = (0, 1,0)
scene.lights = [vector(0, -10, 0)]
xaxis = curve(pos=[ (20, 0, 0), (0, 0,0)], color=(1,0,0))
yaxis = curve(pos=[ (0, 10, 0), (0, 0,0)], color=(0,1,0))
zaxis = curve(pos=[ (0, 0, 10), (0, 0,0)], color=(0,0,1))
balls=[]
seed(4)
nballs=5
vx=zeros((nballs,))
vz=zeros((nballs,))
for i in range(nballs):
rad=uniform(0.5,2)
vx[i]=uniform(1,2)
balls.append(sphere(radius=rad,pos=(i,0,10)))
'''gravity constant'''
g = -9.8
'''change in time'''
dt = 0.01
'''initial z velocity'''
velz = 0.0
'''initial time'''
t = 0.0
print vx,vz
while 1:
rate(60)
'''time'''
t = t + dt
vz=vz+g*dt
for n,i in enumerate(balls):
i.pos.z = i.pos.z + vz[n] * dt #+ 0.5 * g * dt * dt
i.pos.x = i.pos.x + vx[n] * dt
if i.pos.z < i.radius:
vz[n] = -vz[n] * 0.75
vx[n] = vx[n] * 0.75
i.pos.z = i.radius
print vz[3:5]
[edit]
01/30/2008 Hit
from visual import *
from random import uniform,seed
def contact(ball1,ball2):
hit=false
#calculate distance
a=ball1.pos.z-ball2.pos.z
b=ball1.pos.x-ball2.pos.x
d=ball1.pos.y-ball1.pos.y
c=((a*a)+(b*b)+(d*d))**0.5
#sum radii
r=ball1.radius+ball2.radius
#determine if there is a hit
if c <= r:
hit=true
return hit
scene.forward = (0, -0.5,1)
scene.up = (0,0,1)
scene.forward = (0, 1,0)
scene.lights = [vector(0, -10, 0)]
xaxis = curve(pos=[ (20, 0, 0), (0, 0,0)], color=(1,0,0))
yaxis = curve(pos=[ (0, 10, 0), (0, 0,0)], color=(0,1,0))
zaxis = curve(pos=[ (0, 0, 10), (0, 0,0)], color=(0,0,1))
balls=[]
seed(4)
nballs=5
vx=zeros((nballs,))
vz=zeros((nballs,))
for i in range(nballs):
rad=uniform(0.5,2)
vx[i]=uniform(1,2)
balls.append(sphere(radius=rad,pos=(i,0,10)))
'''gravity constant'''
g = -9.8
'''change in time'''
dt = 0.01
'''initial z velocity'''
velz = 0.0
'''initial time'''
t = 0.0
print vx,vz
while t==0.0:
rate(60)
'''time'''
t = t + dt
vz=vz+g*dt
for n,i in enumerate(balls):
i.pos.z = i.pos.z + vz[n] * dt #+ 0.5 * g * dt * dt
i.pos.x = i.pos.x + vx[n] * dt
if i.pos.z < i.radius:
vz[n] = -vz[n] * 0.75
vx[n] = vx[n] * 0.75
i.pos.z = i.radius
for n in range(nballs-1):
for m in range (n+1,nballs):
ball1=balls[n]
ball2=balls[m]
print 'contact',n,m
if contact(ball1,ball2):
print n,m
[edit]
02/01/2008 Hit and Move
from visual import *
from random import uniform,seed
def contact(ball1,ball2):
hit=false
#calculate distance
a=ball1.pos.z-ball2.pos.z
b=ball1.pos.x-ball2.pos.x
d=ball1.pos.y-ball1.pos.y
c=((a*a)+(b*b)+(d*d))**0.5
#sum radii
r=ball1.radius+ball2.radius
#determine if there is a hit
if c <= r:
hit=true
return hit
scene.forward = (0, -0.5,1)
scene.up = (0,0,1)
scene.forward = (0, 1,0)
scene.lights = [vector(0, -10, 0)]
xaxis = curve(pos=[ (20, 0, 0), (0, 0,0)], color=(1,0,0))
yaxis = curve(pos=[ (0, 10, 0), (0, 0,0)], color=(0,1,0))
zaxis = curve(pos=[ (0, 0, 10), (0, 0,0)], color=(0,0,1))
balls=[]
seed(4)
nballs=5
vx=zeros((nballs,))
vz=zeros((nballs,))
for i in range(nballs):
rad=uniform(0.5,2)
vx[i]=uniform(1,5)
balls.append(sphere(radius=rad,pos=(i*2,0,10)))
'''gravity constant'''
g = -9.8
'''change in time'''
dt = 0.01
'''initial z velocity'''
velz = 0.0
'''initial time'''
t = 0.0
print vx,vz
while 1:
rate(60)
'''time'''
t = t + dt
vz=vz+g*dt
for n,i in enumerate(balls):
i.pos.z = i.pos.z + vz[n] * dt #+ 0.5 * g * dt * dt
i.pos.x = i.pos.x + vx[n] * dt
if i.pos.z < i.radius:
vz[n] = -vz[n] * 0.75
vx[n] = vx[n] * 0.75
i.pos.z = i.radius
for n in range(nballs-1):
for m in range (n+1,nballs):
ball1=balls[n]
ball2=balls[m]
print
if contact(ball1,ball2):
i.pos.x=i.pos.x+0.5
i.pos.z=i.pos.z+0.5
print
[edit]
02/20/2008 Explicit Groundwater Equation in Excel
First the equation for transient groundwater flow was derived from dv/dt = sum Q. The final equation is:
hit+1 = head at new time step hi,hi-1,hi+1 = head at old time step Δt = change in time K = hydraulic conductivity Φ = porosity hit+1 = hi + ((-2hi^2 + hi-1^2 + hi+1^2)*Δt*K) / Φ
An Excel spreadsheet was used to make a 1-D model. Eleven nodes were used, and the equation was run until the model reached a nearly steady state (t=540). The parameters were set as follows:
for t = 0, hx = 5 for t > 0, h(at x=0) = 30, h(at x=11) = 5 K = 0.001 Φ = 0.3 Δt = 3
[edit]
02/22/2008 Explicit Groundwater Model - VPython
The same formula used in the Excel model was incorporated in to VPython. The code is as follows:
from visual import *
from visual.graph import *
def plot_line(data_array):
glist = []
ct = -1
for i in data_array:
ct+=1
glist.append((ct,i))
c = gcurve(pos=glist)
'''set up parameters'''
K = 0.001
porosity = 0.3
dt = 3
s = K * dt / porosity
nsteps = 350#number of timesteps
dhdt_limit = 0.001
'''initial conditions'''
h = ones((11,)) * 5.0
t = 0
for t in range(1,nsteps+1):
t += 1
print '''timestep =''', t
'''set boundary conditions'''
h[0] = 30.0
h[10] = 5.0
'''solve for heads'''
#print h[1:10]
h[1:10] = h[1:10] + (-2.0*h[1:10]*h[1:10] + h[0:9]*h[0:9] + h[2:11]*h[2:11]) * s
print h
plot_line(h)
[edit]
02/27/2008 Explicit and Implicit Temperature Model
- Create two graphs with both implicit and explicit 1d solutions for the transient temperature model where
Compare implicit and explicit solutions dt = 1000 K = 0.01 dx = 10 t = 1000
- Create two graphs with both implicit and explicit 1d solutions for the transient temperature model where
dt = 1000 K = 0.01 dx = 10 t = 20000
- Do both methods give the same result?
At t = 1000 the explicit model has slightly lower values than the implicit model between nodes 1 - 4. The data is as follows:
node: 1 2 3 4 5 6 7 8 9 10 11
temp ex: 30 5 5 5 5 5 5 5 5 5 5
temp im: 30 7.098 5.176 5.015 5 5 5 5 5 5 5
At t = 20000 the values for the explicit model and the implicit model are almost exactly the equal. The data is as follows:
node: 1 2 3 4 5 6 7 8 9 10 11
temp ex: 30 20.29 12.8 8.28 6.14 5.33 5.08 5.02 5.002 5 5
temp im: 30 20.20 12.82 8.45 6.33 5.46 5.14 5.04 5.01 5.002 5
- Change the timestep used in both models and determine what is the stability criteria for both models.
Explicit model: the model is stable when dt is less than or equal to 4500 Implicit model: the model is stable at all timesteps
- Determine how many timesteps will be necessary to get to steady state using the explicit solution.
I used 94 timesteps at dt = 4500 to reach steady state in the explicit model. I defined steady state as: the change in temperature from one timestep to the next was equal to 0.01
- What are the advantages and disadvantages of using the implicit and explicit models.
Advantages:
The implicit model can handle a timestep of more than 4500.
The implicit model reaches steady state faster.
The explicit model is simple and easy to use.
Disadvantages:
The implicit model has a complicated formula.
The implicit model solves its formula with multiple iterations.
The explicit model reaches numerical instability at a timestep greater than 4500.
[edit]
