kill(all)$

maxima:numerical phenomena,stirred tank heater,version 1.0
peter.vlasschaert@gmail.com,01/03/2023

info :Material & Energy Balances.

load packages :
==========

load(lrats)$

Basic Thermodynamics : intro H,enthalpy,U = internal energy .
----------------------------------------------------------------------------------------
see also :wikipedia , first law of thermodynamics,cycles,state...........

total energy : 
==========
1e) tot[E]   = total energy
2e) totb[E] = total energy per unit mass
3e) enthalpy , H

enthalpy :used general for chemical thermodynamics ,ΔH,T,ΔS,ΔG,ΔU
H = enthalpy                        , H = U+p*V, "thermo dynamic function"
T = temperature                   , SI units , Kelvin.
S = Entropy                          ,  chaos in the system under consideration.
G = Gibbs Free Energy         , G =H-T*S,"thermo dynamic function"
U = internal Energy              , chemical process , here for this model only contribution for
                                                 total energy, potential energy and kinetic energy are neglected.
see : wikipedia , Gibbs Free Energy.,explain : Δ G = 0,equilibrium state.

rem :remenber simple gas law , p.V =R*n*T→ p = R.n/V*T =R*ρ*T
R = gas cte = 8,31 J/(mol*K)

→ 1e) tot[E]   = total energy
=====================

s1a:print("tot[E]=internalenergy +kineticenergy + potentialenergy")$

s1a1:eval_string(args([s1a])[1]);

info : Get formula from a "print",statement,use "args"

s1a2:lhs(s1a1) =U+1/2*m*v^2+m*g*h;

s1a3:[part(rhs(s1a1),3) =part(rhs(s1a2),3),part(rhs(s1a1),1)=part(rhs(s1a2),2),
    part(rhs(s1a1),2)=part(rhs(s1a2),1)];

info : first in list ,s1a3[1] or first[s1a3],.....

print("internal energy","U =","unit = Joule=J","[N*m]")$
print("velocity"            ,"v=",    "unit=","[m/s]")$
print("acceleration"    ,"a =",  "unit =","[m^2/s]")$
print("height"              ,"h="  ,"unit=","[m]")$
print("mass"               ,"m="  ,"unit=","[Kg]")$

to make tableform : use args (use list)-> eval_string →matrix → table .
info : use args (use list)-→ eval_string ,see above :s1a1

info : a=g ,acceleration on earth surface.
----------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------
"newton second law",F = m*a = Force = mass*acceleration , from N =  [kg]*[m/s^2]           =[kg*m/s^2    ]
"energy                    " E = Energy = Force *distance             , from  Joule =[kg*m/s^2]*[m]  =[Kg*m^2/s^2]
-----------------------------------------------------------------------------------------------------------------------------------------

→ 2e) totb[E] = total energy per unit mass
=================================

 Energy ,per unit mass ?

info :use molecular weight mass  can be converted  to "mole"
used : in general by chemical people.

s1b1:totb[E] =expand(rhs(s1a2)/m);

info : totb[E] = totalbar Energy

set: Ub = U/m , b=bar

s1b2:subst([U/m=Ub],s1b1);

sometimes : F = force , be also F =[m^3/s],flowrate,specify by the problem.

other symbols : we will use for model above

p = pressure         [N/m^2]
Q= heat                [N*m]
F = flowrate         [m^3/s] ,volumeflow ,volume flowrate
ρ = density           [Kg/m^3]
W = work             [N*m]
T = Temperature [K],[°C],[°F],[°R]*
t  = time                [s]
V = volume           [m^3/s]
H = enthalpy         [J]             , H=U+p*V
?  = massrate         [Kg/s],massflow,  mass flowrate
?  = energyrate      [J/s] ,energyflow, energy flowrate

→ 3e) enthalpy , H
===============

→ enthalpy :

s1c:H =U+p*V;/*enthalpy*/

→ enthalpybar , Hb, per unit mass , mole,..... ;b=bar

s1c1:Hb=expand(rhs(s1c/m));

s1d :reverse(ρ =m/V);/*density*/

s1c2:expand(lratsubst([U/m=Ub],s1c1))$
s1c2:expand(lratsubst([s1d],s1c2));

Hb = enthalpyperunitmass                         [J/Kg]             , Hb=Ub+p/ρ
Ub = internalenergyperunitmass               [J/Kg]              , Ub =U/m

info :
-------

------------------------------------------------------------------------------------

*K =Kelvin ,C=Celsius,F= Fharenheit ;R =Ramure,see:wikipedia
? not specific symbol
------------------------------------------------------------------------------------

info : index notation ,i =input, o = output
------------------------------------------------------

print("F[i] = input flowrate → ",F[i]);

use also:F[i] =Fi,others, Fo ,po,pi,........................

Build model :(s)tirred (t)ank (h)eater. (STH).
----------------------------------------------------------

→ short discription :
===============

---------------------------------------------------------------------------------------------------------------------------
Z     = heating system  ( =Q = heat added to system per unit time :[J/s])
F[i] = input flowrate [m^3/s]
T[i] = input temperature 
V     = Volume
----------------------------------------------------------------------------------------------------------------------------
condition system :perfect mixed stirred tank 
system                 : SISO , single input ,single output
temperature        : T=T(t), need find → dT/dt = ? ("ode : (o)rdinary (d)ifferential (e)quations")
----------------------------------------------------------------------------------------------------------------------------

theory:

→ ACCUM = IN -OUT +Work (on the system) +Heat (on the system)

A) Material Balance :
---------------------------
Work (on the system) =0
Heat (on the system)  = 0

→ kg/s  =[kg/m^3*m^3/s] =[density * flowrate] =[density*velocity*surface tank ] =[Kg/m^3*m/s*m^2]

z1a:ACCUM =IN -OUT;

z1b:[ACCUM ='diff(m,t),IN =F[i]*ρ[i],OUT = F*ρ];

z1c:subst(z1b,z1a);

→use  s1d

z1d:solve(s1d,m)[1];

z1c1:subst([z1d],z1c);

B) Energy Balance.
---------------------------

z2a:ACCUM = IN -OUT +Work +Heat;

first assumption :

→ suppose  :[potentialenergy = 0,kineticenergy = 0]

use :s1a3

zb1:makelist(reverse(s1a3[i]),i,1,length(s1a3));

zb1:subst([potentialenergy = 0,kineticenergy = 0],zb1);

use :s1a1

zb2:subst([potentialenergy = 0,kineticenergy = 0],s1a1);

zb2:lhs(s1b2) =lhs((zb1[1]))/m;

lhs; ACCUM,energy balance
..............................................

zb3:ACCUM ='diff(U,t);

rhs; energy balance
................................

zc1: [IN =ρ[i]*F[i]*Ub[i],OUT =ρ*F*Ub];

use : z2a

zc2:subst(zc1,z2a);

 ? Heat :

zc2:subst([Heat=Q],zc2);

? Work :

zch:Work =W[shaft]+W[i]-W[o];

explain : Work(on the system)
W[i] -W[o] = "result energy because difference in the pressure"
energy to bring fluid in and drive them out after heating and mixing

zch1:[W[i] =p[i]*F[i],W[o] =F*p];/*W[i],W[o] =[m^3/s*N/m^2] =[J/s]*/

zc2:lratsubst([zch],zc2);

zc2:subst([zch1],zc2);

zc2;

final : Energy Balanc.
....................................

zd1:subst([zb3],zc2);

use :enthalpy ,s1c2

ze1:reverse(s1c2);

zea:listofvars(ze1);

→ index

total Energy Balance : U → H
..................................................

zeai:makelist([i],j1,length(zea));

zebi:makelist(f[i],f,zea);

zeaii:map("=",zea,zebi);

use :ze1

ze2:subst(zeaii,ze1);

haa:lratsubst([ze1],zd1);

→ enthalpy : rhs

haa:lratsubst([ze2],haa);

→ enthalpy :lhs

hxx:H=U+p*V;

hxxv:listofvars(hxx);

hxd:diff(hxx);

sh1: solve(hxd,del(U));

l1:subst([hk=H],del(hk)='diff(hk,t));

use :sh1

shll:makelist(subst([hk=hxxv[i]],del(hk)='diff(hk,t)),i,1,length(zea));

shkk:subst(shll,sh1);

assumption :p.V , (= neglected ),fluid
1e) ('diff(p,t,1)) = 0, Δp not to large, equation → V=0
2e) ('diff(V,t,1)) = 0,V =cte                                → p=0

shkkx:subst([V=0,p=0],shkk);

use :haa

shkU:haa;

use : fluid with assumption , shkkx

shkU:subst(shkkx,shkU);

assumption : W[shaft] =0

shkshaft:subst([W[shaft] =0],shkU);/*no mechanical work*/

rem : H rhs , enthalpy per unit of mass

what about : H =H(t),heat capacity =C[p]
------------------------------------------------------

google : heat capacity polynomial coefficients
→ C P = a + b × T + c × T 2 + d × T 3 + e × T 4 for temperature

approximation : heat capcity = constant

1e) total enthalpy
------------------------

z1Hb : H =V*ρ*Hb;/*V*ρ =m*/

z1Hbi:H[i]=V*ρ*Hb[i];

2e) no phase change
---------------------------

→ H

z1nofase1: Hb ='integrate(C[p],tt,T[ref],T);

→ H[i]

z1nofase2: Hb[i] ='integrate(C[p],tt,T[ref],T[i]);

3e) combine 1e) and 2e)
--------------------------------

total Energy Balance ,use 1e) & 2e)
---------------------------------------------------------

→ use : shkshaft , input ,z1Hb

zy1:lratsubst([z1Hb,z1Hbi],shkshaft);

zy1:subst([z1nofase1],lhs(zy1)) =subst([z1nofase2],rhs(zy1));

zy1:subst([z1nofase1],(zy1));

suppose :1e) density constant ρ[i] =ρ
                2e) volume constant F[i]=F , follows from above: Material Balance.

zy2:subst([ρ[i] =ρ,F[i]=F],zy1);

zy2:factorsum(expand(zy2));

because T =T(t) ,T[ref] =constant, differentation =0

depends(T,t);

zy3:ev(zy2,diff)$

zy3:expand(subst([diff(T[ref],t)=0],zy3));

zyh1:ratcoeff(lhs(zy3),diff(T,t));

zy3:expand(zy3/zyh1);

solve problem : 
--------------------
→specify                 : {V,ρ,C[p]}
→function time      : F,Q; T[i]  → F(t),Q(t),T[i](t)
initial condition : t=0 , T(t) =T(0) = some value

zy4h:[F=F(t),Q=Q(t),T[i]=T[i](t)];
zy4:subst(zy4h,zy3);

zy5:expand(lratsubst([F(t)/V=DE],zy4));
    

zy5:lhs(zy5)=collectterms(rhs(zy5),DE);

final expression : initial condition T(0),given.

zy5:subst([DE=F(t)/V],zy5);

'diff(T,t,1)=Q(t)/(V*C[p]*ρ)+(F(t)*(T[i](t)-T))/V

rem : collectterms only on one part of expression.

--------------------------------------------------------------------------------------------------