Message Boards

WOLFRAM COMMUNITY

7361 Views

4 Replies

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Chemistry Engineering Mathematics

solving systems of equations with given initial conditions

dylan keefe

Posted 10 years ago

I would like to solve {na, nb, nw, T, P} at z=20 given these differential equations and initial conditions where these values are given d := .686; L := 6.73; U := 812; Te := 800; eP := .6; OS := .45; Dp := .0046; nao := 1.345; nbo := 1.412; mwa := .064; mwb := .032; mwi := .014; G := 4(naomwa + nbomwb + nimwi)/(Pid^2); dens := (4GP[z])/(Pidd8.314T[z](na[z] + nb[z] + nw[z] + ni)); H298 := 2(-94470) - (2(-70950)); aa := 5.697; ab := 6.713; aw := 12.13; ba := 0.016; bb := -0.0000000879; bw := .00812; ya := -0.00001185; yb := .00000417; yw := 0; da := .000000003172; db := -.000000002544; dw := 0; H[T[z]] := (H298 ) + (2aw - 2aa - 2ab)(T[z] - 298) + ((2bw - 2ba - bb)/2)(T[z]^2 - 298^2) + ((2yw - 2ya - yb)/3)(T[z]^3 - 298^3) + ((2dw - 2da - db)/4)(T[z]^4 - 298^4) NDSolve[{na'[z] == (-2((Pidd)/4))r, nb'[z] == (-(Pidd)/4)r, nw'[z] == 2((Pidd)/4)r, T'[z] == 1((PidU(Te - T[z])) - (((Pidd)/4)rH))/(na[z]Cpa + nb[z]Cpb + nw[z]Cpw + niCpi), P'[z] == -1((1 - eP)/ eP^3)(G^2/(densOSDp9.81))(((150(1 - eP).32)/(OSDpG)) + 1.75), na[0] == 1.345, nb[0] == 1.412, nw[0] == 0, T[0] == 727.6, P[0] == 2.5}, {na, nb, nw, T, P}, {z, 20}] is what i thought would work but it is not Attachments:

I would like to solve {na, nb, nw, T, P} at z=20 given these differential equations and initial conditions where these values are given

d := .686;
L := 6.73;
U := 812;
Te := 800;
eP := .6;
OS := .45;
Dp := .0046;
nao := 1.345;
nbo := 1.412;
mwa := .064;
mwb := .032;
mwi := .014;
G := 4*(nao*mwa + nbo*mwb + ni*mwi)/(Pi*d^2);
dens := (4*G*P[z])/(Pi*d*d*8.314*T[z]*(na[z] + nb[z] + nw[z] + ni));
H298 := 2*(-94470) - (2*(-70950));
aa := 5.697;
ab := 6.713;
aw := 12.13;
ba := 0.016;
bb := -0.0000000879;
bw := .00812;
ya := -0.00001185;
yb := .00000417;
yw := 0;
da := .000000003172;
db := -.000000002544;
dw := 0;
H[T[z]] := (H298 ) + (2*aw - 2*aa - 2*ab)*(T[z] - 
     298) + ((2*bw - 2*ba - bb)/2)*(T[z]^2 - 
     298^2) + ((2*yw - 2*ya - yb)/3)*(T[z]^3 - 
     298^3) + ((2*dw - 2*da - db)/4)*(T[z]^4 - 298^4)


NDSolve[{na'[z] == (-2*((Pi*d*d)/4))*r, nb'[z] == (-(Pi*d*d)/4)*r, 
  nw'[z] == 2*((Pi*d*d)/4)*r, 
  T'[z] == 1*((Pi*d*U*(Te - T[z])) - (((Pi*d*d)/4)*r*H))/(na[z]*Cpa + 
       nb[z]*Cpb + nw[z]*Cpw + ni*Cpi), 
  P'[z] == -1*((1 - eP)/
      eP^3)*(G^2/(dens*OS*Dp*9.81))*(((150*(1 - eP)*.32)/(OS*Dp*G)) + 
      1.75), na[0] == 1.345, nb[0] == 1.412, nw[0] == 0, 
  T[0] == 727.6, P[0] == 2.5}, {na, nb, nw, T, P}, {z, 20}]

is what i thought would work but it is not

POSTED BY: dylan keefe

4 Replies

Sort By:

Bill Simpson

Posted 10 years ago

1: The help page for NDSolve shows NDSolve[{eqns},{vars},{z,min,max}] and you have only {z,num}, but maybe that is ok. Ignoring that for the moment, evaluating your notebook, note that leading sols= and trailing [[1]] that I added, gives In[45]:= Off[NDSolve::ndsz] In[46]:= sols = NDSolve[{na'[z] == (-2((Pidd)/4))r, nb'[z] == (-(Pidd)/4)r, nw'[z] == 2((Pidd)/4)r, T'[z] == 1((PidU(Te - T[z])) - (((Pidd)/4)rH))/(na[z]Cpa + nb[z]Cpb + nw[z]Cpw + niCpi), P'[z] == -1((1 - eP)/eP^3)(G^2/(densOSDp9.81))(((150(1 - eP).32)/(OSDpG)) + 1.75), na[0] == 1.345, nb[0] == 1.412, nw[0] == 0, T[0] == 727.6, P[0] == 2.5}, {na, nb, nw, T, P}, {z, 20}][[1]] Out[46]= {na -> InterpolatingFunction[{{0.`, 2.0480816649058247`^-10}}, <<<snip>>>} The help page for InterpolatingFunction shows Interpolatingfunction[{domain},table] so it appears that your solution has a domain from zero to 210^-10. Let's see if we can get results from that. In[47]:= naf = na /. sols Out[47]= InterpolatingFunction[{{0.`, 2.0480816649058247`^-10}},<>] In[48]:= naf[210^-10] Out[48]= 1.34499 + 0. I In[49]:= Plot[naf[z], {z, 0, 210^-10}] Out[49]= <<<ApparentlyPerfectlyReasonablePlotSnipped>>> So it appears that within the domain you are getting results and can plot them. Even if I correct the syntax to use {z,0,20} inside NDSolve to tell it you want the solution within that domain the result is unchanged. If I disable your Off[NDSolve::ndsz] so I can see the error messages begin generated by NDSolve it tells me NDSolve::ndsz: At z == 2.0480816649058247`^-10, step size is effectively zero; singularity or stiff system suspected. >> that your system blew up inside NDSolve at 210^-10 and exited. So, can you look at your system, perhaps look at the denominators, check anything else and perhaps see why it would blow up at that point? Almost always, disabling warnings and errors doesn't make the problem go away, it just hides what is needed to understand them.

POSTED BY: Bill Simpson

dylan keefe

Posted 10 years ago

here is the modified code... it seems to be something else perhaps? Attachments:

POSTED BY: dylan keefe

Bill Simpson

Posted 10 years ago

NDSolve demands that all variables be assigned values. r,H,Cpa,Cpb,Cpw,ni,Cpi, have no values. Unfortunately NDSolve doesn't provide an error message clearly stating undefined variables are included and eliminate your H[T[z]] :=... Then see what you get.

POSTED BY: Bill Simpson

dylan keefe

Posted 10 years ago

..

POSTED BY: dylan keefe

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback