[scilab-Users] No reaction

[Jaundre Venter]

Hi

okay i have made the changes and the reaction curves are looking better but
for some reason i still get no reaction on the last three curves. It is just
showing a straight line. can it be that i have choosen to plot the wrong
type of graph? Maybe because of a time value or to much ODE"s?

On Wed, Oct 5, 2011 at 10:24 PM, Jaundre Venter <jaundreventer at gmail.com>wrote:

> Thank you very much Adrien.
>
> always nice if someone can explain to you where your problems are and why.
> Thanks
>
> Was there any other problems you saw that i have to be aware of?
>
>
> On Wed, Oct 5, 2011 at 5:38 PM, Adrien Vogt-Schilb <vogt at centre-cired.fr>wrote:
>
>>  Hi
>>
>> When you use ode, it's ok, if say, dx(1) depends on dx(4).
>> but you still have say that to scilab properly, something like:
>>
>> function dx = f(t,x)
>> dx(6)=((F+Fab+Floss)*(x(2))),  // culture Volume V
>> dx(1)=(((mu)*(x(1)))-(((x(1))/(x(6)))*((dx(6))))*(CO)*(x(2)))), //biomass
>> concentration X
>> and so one
>>
>> note that because i had to know dx(6) to compute dx(1) i just computed
>> dx(6) before dx(1): no problem. and note that i used x(2). The idea of the
>> ode is to compute dx from x!
>>
>> make sure you understand that using dx_6 instead of dx(6), your ODE solver
>> is not updating dx_6 at each time step, it is using the initial and only
>> dx_6 forever. That's why your last varaibles do not move, somehow their
>> speeds are never updated.
>> for instance, dx(6)=((F+Fab+Floss)*(HION)),  // culture Volume V is
>> constance in time (i guess)
>>
>>
>>
>> On 05/10/2011 15:00, Jaundre Venter wrote:
>>
>> Hi Adrien
>>
>> i am new to SCILAB! I just want to say that.
>>
>> Yes dx_1 is equal to dx1 but the only reason why i have programmed it like
>> that is becasue the ODE's looks as follows - (see word file attached that
>> can explain the ODE"s better. with regards to the HION and CO it actually
>> refers to [H+] and CO2 as you said. the only reason why n multiplied HION
>> and CO wit hsome ODE's is because some does have a influence on some ODE's.
>>
>> This is the first time i am working with SCILAB thus i am struggling to
>> understand how SCILAB wants the code so that all 9 ODE's are shown and so
>> that the ODE's that is having a effect on other does happen. I though if you
>> refer to dx(1) for example in a other ODE it means that SCILAb will know the
>> dx(1) has a influence on the other ODE.
>>
>> The main goal of my assesment is to deliver similar results obtained from
>> MATLAB on SCILAB. all i got was how the grpahs should look like and the
>> ODE's.
>>
>> On Wed, Oct 5, 2011 at 2:43 PM, Adrien Vogt-Schilb <vogt at centre-cired.fr>wrote:
>>
>>>  Hi
>>>
>>> Try to isalote your problem
>>> if i understood well, the following code
>>>
>>>
>>> //  initial values
>>> x0=[0.1, 1e-5, 0, 15, 1.16, 100,0,297,0.5]';
>>> t=0:0.005:400;
>>> y=ode(x0, 0, t, f);
>>>
>>>  returns y such that sum(y(6:9,:)>x0) == 0 ?
>>> if this is true, we do not need the plots to solve the problem
>>> can you check that ?
>>>
>>> I believe the f function is erroneous.
>>> It seems that dx_1 should be equal to dx(1) at each time step, and that
>>> HION should be equal to x(2) at each time step, etc.
>>>
>>> in other terms, some of your phisical variables seem to be represented by
>>> to variables (i am guessing HION=[H+] and x(2)=[H+] also) but scilab does
>>> not have any chance to know that.
>>> if my guess is right, you have to rewrite the f function in a way that
>>> eliminates all references to HION, dx_1, dx_6 and so on
>>>
>>>
>>> On 05/10/2011 14:27, Jaundre Venter wrote:
>>>
>>> Hi all
>>>
>>> Can someone please explain to me the following:
>>>
>>> I am busy with a project of simulation the production of penicillin in a
>>> bio reactor. Now i have 9 ODE's which i want to simulate.
>>>
>>> now for some reason the last three graphs i am getting doesn't show any
>>> response what so ever. i am using the following code.
>>>
>>> dx(1)=(((mu)*(x(1)))-(((x(1))/(x(6)))*((dx_6)))*(CO)*(HION)), //biomass
>>> concentration X
>>> dx(2)=((z*(((mu)*(x(1)))-(((F)*(x(1)))/(x(6)))))+(QQ)), //hydrogen ion
>>> concentration H+
>>> dx(3)=((((mupp)*(x(1)))-((K)*(x(3)))-((x(3))/(x(6)))*(dx_6))*(HION)),
>>> //Penicilin concentration P
>>> dx(4)=((-((mu)/(Yxs))*(x(1)))-(((mupp)/(Yps))*(x(1)))-((mx)*(x(1)))+((Fsf)/(x(6)))-((x(4)/(x(6)))*(dx_6))),
>>> //Substrate concentration S
>>> dx(5)=(-(((mu)/(Yxo))*(x(1)))-(((mupp)/(Ypo))*(x(1)))-(((mo))*(x(1)))+((kla)*(cll-(x(5))))-(((x(5))/(x(6)))*(dx_6))),
>>> //dissolved oxygen
>>> dx(6)=((F+Fab+Floss)*(HION)),  // culture Volume V
>>> dx(7)=(((rq1)*(dx_1)*(x(6)))+(rq2)*(x(1))*(x(6))), //Heat generation Qrxn
>>> dx(8)=((((F)/(sf))*(Tf-(x(8))))+(1/((x(6))*(pcp)))*(QT)),  //
>>> Temperature T
>>> dx(9)=(((a1)*(dx_1))+((a2)*(x(1)))+(a3)),  //  CO2 evolution, CO2
>>> endfunction
>>>
>>> now when i ask for plotting the graphs i am using the following.:
>>>
>>> //  initial values
>>> x0=[0.1, 1e-5, 0, 15, 1.16, 100,0,297,0.5]';
>>> t=0:0.005:400;
>>> y=ode(x0, 0, t, f);
>>>
>>> // the plots of each variable
>>> da.title.text="BIOMASS CONCENTRATION"
>>> da.x_label.text="Time, hours";
>>> da.y_label.text="X,g/l ";
>>> scf(1);clf; //Opens and clears figure 1
>>> plot(t,y(1,:))
>>>
>>> da.title.text="HYDROGEN ION H+ CONCENTRATION"
>>> da.y_label.text="H+,mol/l ";
>>> scf(2);clf; //Opens and clears figure 2
>>> plot(t,y(2,:))
>>>
>>> da.title.text="PENICILLIN CONCENTRATION"
>>> da.y_label.text="P,g/l ";
>>> scf(3);clf; //Opens and clears figure 3
>>> plot(t,y(3,:))
>>>
>>> da.title.text="SUBSTRATE CONCENTRATION"
>>> da.y_label.text="S,g/l ";
>>> scf(4);clf; //Opens and clears figure 4
>>> plot(t,y(4,:))
>>>
>>> da.title.text="DISSOLVED OXYGEN CONCENTRATION"
>>> da.y_label.text="C_l,g/l ";
>>> scf(5);clf; //Opens and clears figure 5
>>> plot(t,y(5,:))
>>>
>>> da.title.text="CULTURE VOLUME"
>>> da.y_label.text="V,l";
>>> scf(6);clf; //Opens and clears figure 6
>>> plot(t,y(6,:))
>>>
>>> da.title.text="HEAT OF REACTION"
>>> da.y_label.text="Qrxn,cal";
>>> scf(7);
>>> clf; //Opens and clears figure 7
>>> plot(t,y(7),:)
>>>
>>> da.title.text="TEMPERATURE"
>>> da.y_label.text="T,Kelvin";
>>> scf(8);
>>> clf; //Opens and clears figure 8
>>> plot(t,y(8),:)
>>>
>>> da.title.text="CO2 EVOLUTION"
>>> da.y_label.text="CO2,mmol/l/";
>>> scf(9);
>>> clf; //Opens and clears figure 9
>>> plot(t,y(9),:)
>>>
>>> Am i doing something wrong? before the ODE's i have just programmed the
>>> initial values and constants :
>>>
>>> funcprot(0);
>>> function dx = f(t,x)
>>> K1=1.0e-10         //mol/l
>>> K2=7.0e-05         //mol/l
>>> Kx=0.15            // Contois saturation constant, g/l
>>> Kox=2e-02          // oxygen limitation constant
>>> mux=0.092          // maitenance coefficient on subsrate
>>> p=3                //constant
>>> Kp=0.0002          //   inhibition constant
>>> Kop=2e-02      // oxygen limitation constant
>>> K=0.04     //  Penicillin hydrolysis constant, per h
>>> Yxs=0.45  //   yield constant,g biomass/g glucose = dimensionless
>>> Yps=0.90  //   yield constant, g pinicillin/ g glucose = dimensionless
>>> mx= 0.014  //   Maintenance coefficient on substrate, per h
>>> Yxo=0.04  //   yield constant, g biomass/g oxygen = dimensionless
>>> Ypo=0.20  //   yield constant, g penicillin/g oxygen= dimensionless
>>> mo= 0.467  //   maintenance coefficient of oxygen, per h
>>> mup=0.0005  // specific rate of penicilline production (per h)
>>> sf=600 // Feed substrate concentration, g/l
>>> kla=23     // function of agitation power input and oxugen flow rate,
>>> dimensional
>>> cll=1.16   //  dissolved oxygen concentration, g/l
>>> Cab=3      // concentrations in both solutions
>>> Fa=5     // acid flow rate, l/h  !!
>>> Fb=5      // base flow rate, l/h !!
>>> delta_t=0.01   //  time step in digital PID controller - arbitrary
>>> value!!!
>>> z=10e-5     // constant
>>> F=0.042       //  feed substrate flow rate l/h
>>> T0=273         //  temperature at freezing, K
>>> Tv=373       //  temperature at boiling
>>> T=298  //  feed temp of substrate
>>> h=(2.5e-4)     //  constant
>>> Floss=(x(6)*(h)*(exp(5)*((T-T0)/(Tv-T0))))
>>> Fab=Fa+Fb  // volume increase due to influx of acid Fa and base Fb
>>> Fsf=((sf)*(F))
>>> kg= 7e-3 //  Arrhenius constant for growth
>>> kd=10e33  //  Arrhenius constant for cell death
>>> Eg= 5100  //  Activation energy for growth, cal/mol
>>> Ed= 50000  //  Activation energy for cell death, cal/mol
>>> R= 1.987  //  gas constant, cal/mol k
>>> T= 297  //  Temperature
>>> RT= R*T
>>> alpa= 70  //  constant in Kla
>>> betha= 0.4  //  constant in Kla
>>> Pw= 30  //  Agitation power input, W
>>> fg= 8.6  //  Flow rate of oxygen
>>> V=100  //  Volume
>>> QE= ((kg*exp(-(Eg/RT)))-(kd*exp(-(Ed/RT))))
>>> kla= alpa*((sqrt(fg)*(Pw/x(6)))^betha)
>>> mu
>>> =(((mux)/(1+((K1)/(x(2)))+((x(2))/(K2))))*((x(3))/(((Kx)*(x(1)))+(x(3))))*((x(5))/(((Kox)*(x(1)))+(x(5))))*(QE))
>>> //  Specific growth rate
>>> mupp =
>>> ((mup)*((x(4))/((Kp)+(x(4))+(x(4)^2)/(K1)))*((x(5)^p)/((Kop)*(x(1)))+(x(5)^p)))
>>> // Specific penicillin production rate
>>> B =(((1e-14/x(2)-x(2))*x(6)-Cab*(Fa+Fb)*delta_t)/(x(6)+(Fa+Fb)*delta_t))
>>> QQ =((-B+sqrt(B^2+4e-14))/2-(x(2)))*(1/delta_t)
>>> dx_6 = (F+Fab+Floss) //Culture Volume V
>>> dx_1 = (((mu)*(x(1)))-((x(1))/(x(6)))*(dx_6)) //biomass concentration X
>>> rq1 = 60  //  yield of heat generation, cal/g biomass
>>> rq2 = 1.6783e-4  //  Constant, cal/g biomass h
>>> Tf = 296  //  substrate feed temperature, Kelvin
>>> a = 1000   //  heat transfer coefficient of cooling/heating liquid, cal/h
>>> degree C
>>> b = 0.60   //  constant
>>> Fc=0.1   //  Cooling water flow rate, not sure about value, l/h
>>> pcCpc = 1/2000   //  Density times heat capacity of cooling liquid, per l
>>> degree C
>>> pcp = 1/1500   //  density times heat capacity of medium
>>> QT = ((x(7)-(((a)*(Fc^b+1))/((Fc)+((a)*(Fc^b))/2*pcCpc))))
>>> a1=0.143  //  constant relating CO2 to growth, mmol CO2/g biomass
>>> a2=4e-7  //  Constant relating CO2 to mainteneance energy, mmol CO2/g
>>> biomass h
>>> a3=1e-4 //  Constant relating CO2 to penicillin production, mmol CO2/l h
>>> CO= (((a1)*(dx_1))+((a2)*(x(1)))+(a3)),  //  CO2 evolution, CO2
>>> HION=((z*(((mu)*(x(1)))-(((F)*(x(1)))/(x(6)))))+(QQ))
>>>
>>> Thanks.
>>>
>>>
>>>
>>>   --
>>> Adrien Vogt-Schilb (Cired)
>>> Tel: (+33) 1 43 94 *73 77*
>>>
>>
>>
>>
>> --
>> Adrien Vogt-Schilb (Cired)
>> Tel: (+33) 1 43 94 *73 77*
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20111006/a39e550f/attachment.htm>