Hi,
I've been testing sympy for a phd project I'm working on and I've run
into some challenges that I need some advice on. Consider the first
part in the example code below, where the goal is to mix a function
and some calculated parameter, utilizing a thermodynamic package which
calculates some derivatives. First I define some constants/
expressions which are necessary for my calculations and then some
implicit functions representing what I need from the thermodynamic
package (to make sure that I get the derivatives correctly). Finally,
I combine them into the expression rNH3. Numerical values of this
expressions and its derivatives with respect to T,V,n are needed for
the calculations in my main routine (reactor calculations).
In "CODE PART 2" I've given an example of how I calculate the
quantity rNH3 and substitute in the values from the constants in "CODE
PART 1" in addition to the ones from the thermodynamic package. In
"CODE PART 3" I calculate the derivative of rNH3 with respect to T and
substitute in the values of the constants and the thermodynamic
package.
Questions:
1) Can the code below be simplified in a manner that I still achieve
my goal?
2) In "CODE PART 3" I need to make sure that I substitute in the
expressions of thermodynamic package derivatives before I substitute
in the non-derivatives of the thermodynamic package. Doing this in the
opposite sequence will result in wrong results. This is probably due
to ignorance, but does the subs function work as a pure textual search
and replace of the expressions?
3) If there an easier path for substituting in values into the
expressions? It would be nice to somehow tell a_NH3, a_N2, a_H2 their
numerical values of the expressions and derivatives. Example:
a_NH3 = 1.0 #
a_NH3.diff(T) = 3.4 #
a_NH3.diff(V) = 1.0 #
and when the evalf function is called these values should be picked
up. Comments?
Thanks in advance!
----- CODE PART 1 -------
rk0 = 8.849e14 # Rate constant for
Arrhenius equation [kmol/h]
kJprkcal = 4.1868 # Conversion factor -
kJ / kcal
rE = 40765.0 * kJprkcal # kJ / kmol
rR = 8.314 * 1e3 # kJ / kmol kK
rc = 0.5 # Constant -
varies from 0.5 - 0.75
k,k0,E,R,T,V,c = var('k k0 E R T V c')
n = var('n')
k = k0 * exp(-E / R / T) # Arrhenius equation
a_NH3 = Function('a_NH3')(T, V, n) # Activity of
ammonia. Retrieve from thermodynamic package
a_N2 = Function('a_N2') (T, V, n) # Activity of
Nitrogen. Retrieve from thermodynamic package
a_H2 = Function('a_H2') (T, V, n) # Activity of
Hydrogen. Retrieve from thermodynamic package
Ka = Function('Ka')(T) # Equilibrium
constant. Retrieve from thermodynamic package
rNH3 = 2*k*(Ka**2.0*a_N2*(a_H2**3.0/a_NH3**2.0)**c -
(a_NH3**2.0/a_H2**3.0)**(1.0-c))
---- CODE PART 2 -----
def rate(t):
substconsts = [(k0, rk0),(E, rE), (R, rR), (c, rc)]
a = getactivity()
Ka_ = getKaFunctionPtr()
rNH3 = rNH3.subs([(a_NH3, a[0]),(a_N2, a[1]), (a_H2,
a[2])])
rNH3 = rNH3.subs([(Ka, Ka_(t))])
rNH3 = rNH3.subs([(T, t)])
rNH3 = rNH3.subs(substconsts)
return 1e3*rNH3.evalf()
-- CODE PART 3 ----
def rate_x(t):
a = getactivity()
a_x = getactivity_x()
Ka_ = getKaFunctionPtr()
Ka_x = getKa_xFunctionPtr()
substconsts = [(k0, rk0),(E, rE), (R, rR), (c, rc)]
rNH3_x = rNH3.diff(T)
rNH3_x = rNH3_x.subs([(a_NH3.diff(T), a_x[0]),
(a_N2.diff(T), a_x[1]), (a_H2.diff(T), a_x[2])])
rNH3_x = rNH3_x.subs([(Ka.diff(T), Ka_x(t))])
rNH3_x = rNH3_x.subs([(a_NH3, a[0]),(a_N2, a[1]), (a_H2,
a[2])])
rNH3_x = rNH3_x.subs([(Ka, Ka_(t))])
rNH3_x = rNH3_x.subs([(T, t)])
rNH3_x = rNH3_x.subs(substconsts)
return 1e3 * rNH3_x.evalf()
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
