# Diff of /trunk/doc/cookbook/onedheatdiff_var001.py

revision 2396 by ahallam, Thu Apr 16 04:24:07 2009 UTC revision 2397 by ahallam, Fri Apr 24 04:24:28 2009 UTC
# Line 33  import os #This package is necessary to Line 33  import os #This package is necessary to
33
34  ##ESTABLISHING VARIABLES  ##ESTABLISHING VARIABLES
35  #PDE related  #PDE related
36    mx = 500 # model lenght
37    my = 100 # model width
38    ndx = 500 # steps in x direction
39    ndy = 1 # steps in y direction
40
41  q=0 #our heat source temperature is now zero  q=0 #our heat source temperature is now zero
42  Tref=50.e6 #the starting temperature of our iron bar  Tref=2273 # Kelvin #the starting temperature of our intrusion
43  rho=2.6e6  rho = 2750 #kg/m^{3} density
44  eta=75.  cp = 790 #j/kg specific heat
45  kappa=240.  rhocp = rho*cp  #DENSITY * SPECIFIC HEAT
47    kappa=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY
48  #Script/Iteration Related  #Script/Iteration Related
49  t=0 #our start time, usually zero  t=0. #our start time, usually zero
50  tend=100.#the time we want to end the simulation  tday=10*365. #the time we want to end the simulation in days
51  h=0.25 #size of time step  tend=tday*24*60*60
52    outputs = 400 # number of time steps required.
53    h=(tend-t)/outputs #size of time step
54
55    print "Expected Number of Output Files is: ", outputs
56    print "Step size is: ", h/(24.*60*60), "days"
57
print "Expected Number of Output Files is: ", (tend-t)/h
58
59  i=0 #loop counter  i=0 #loop counter
60  save_path = "data/onedheatdiff_var001" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work  save_path = "data/onedheatdiff_var001" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work
61
62  #... generate domain ...  #... generate domain ...
63  rod = Rectangle(l0=0.05,l1=.01,n0=500, n1=1)  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
64  # extract finite points  # extract finite points
65  x=rod.getX()  x=model.getX()
66  #... open PDE ...  #... open PDE ...
67  mypde=LinearPDE(rod)  mypde=LinearPDE(model)
68  mypde.setSymmetryOn()  mypde.setSymmetryOn()
69  mypde.setValue(A=kappa*kronecker(rod),D=rho/h,d=eta,y=eta*Tref)  mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
70
71  # ... set initial temperature ....  # ... set initial temperature ....
72  T= -1*Tref*whereNegative(x[0]-0.025)+0.01*Tref*wherePositive(x[0]-0.025)+0.99*Tref*wherePositive(x[0]-0.04)  bound = x[0]-mx/(ndx/250.)
73    T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)
74
75  saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)  saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)
76
# Line 66  saveVTK(os.path.join(save_path,"data%03d Line 78  saveVTK(os.path.join(save_path,"data%03d
78  while t<=tend:  while t<=tend:
79        i+=1        i+=1
80        t+=h        t+=h
81        mypde.setValue(Y=rho/h*T)        mypde.setValue(Y=rhocp/h*T)
82        T=mypde.getSolution()        T=mypde.getSolution()
print T
83        saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)        saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)
84
85    # iteration var 2
86    #Tl = 0
87    #Tr = Tref
88    #while Tl < Tr*0.8:
89          #mypde.setValue(Y=rhocp/h*T)
90          #T=mypde.getSolution()
91          #i+=1
92          #x=rod.getX()
93          #Tl= x[0]
94          #Tr= x[
95
96    #print "Finish temp balance in:", i*h/(24.*60*60), " days"

Legend:
 Removed from v.2396 changed lines Added in v.2397