Message Boards Message Boards

0
|
1395 Views
|
1 Reply
|
0 Total Likes
View groups...
Share
Share this post:

NDSolve second-order differential eqns with boundary condition

Posted 10 years ago
Change the yo (3rd boundary condition: y[0]) to meet the target(4th boundary condition): funy'[1] + pey*bd*funy[1]
=0, and get the answer yo (0<yo<1)
Initial yo=0.3
but,the program cannot run, please give me some help, many thans------------------------------------Clear[fx, fy, cxo, cy1, h, pey, noofc, hox, ans];

(*Input data for the calculation of the outlet concentrations*)
{
  {Flowrate of Cu Raff (m3/h), fx = 5}
} ;
{
  {Flowrate of barren solvent (m3/hr), fy = 1}
} ;
{
  {Subscript[U, 3] Subscript[O, 8] concentration in Cu Raff (mg/L),
   cxo = 267}
} ;
{
  {Subscript[U, 3] Subscript[O, 8]
     concentration in barren solvent (mg/L), cy1 = 130}
} ;
{
  {Effective column height (m), h = 9}
} ;
{
  {Peclet number of organic phase , pey = .7}
} ;
{
  {no of compartments, noofc = 13}
} ;
{
  {Height of transfer unit (m), hox = 1.65855537601783531}
} ;

(* Equilibrium constants represented by Langmuir model, \
cx^*=(a*cy)/(b*cy+1) *)
a = 0.00982089325801922;
b = -0.0002962283775423278;

(* Press once you complete the input section *)

nox = h/hox;
uxuy = fx/fy; dc = h/noofc; bd =
h/dc; SuperStar = (a/cxo - b)^-1; SuperStar[cx1] = (a*cy1)/(
b*cy1 + 1);
b1 = cxo - SuperStar[cx1]; b2 = SuperStar - cy1; a1 =
a*(SuperStar - cy1); a2 = a*cy1; a3 = b*(SuperStar - cy1);
a4 = (b*cy1) + 1;
yo = .3; parta[yo_] :=
Module[{useless},
  ini = NDSolve[{-x' - (nox/
         b1)*(b1*x + SuperStar[cx1] - (a1*y + a2)/(
          a3*y + a4)) == 0,
     y'' +
       pey*bd*y' + ((nox*pey*bd*uxuy)/
         b2)*(b1*x + SuperStar[cx1] - (a1*y + a2)/(
          a3*y + a4)) == 0, y'[0] == 0, x[0] == 1,
     y[0] == yo}, {x, y}, {z, 0, 1}]; funx[z_] = ini[[1, 1, 2]];
   funy[z_] = ini[[1, 2, 2]]; funy'[1] + pey*bd*funy[1]];
try = FindRoot[parta == 0, {ans, 0, 1}, MaxIterations -> 350,
   AccuracyGoal -> 6];
Plot[{funx, funy}, {z, 0, 1}, Frame -> True,
FrameLabel -> {"Dimensionless length (Z) ",
   "Dimensionless conc (X,Y)"}, PlotStyle -> AbsoluteThickness[2],
PlotRange -> All, ImageSize -> 350]
Print["y0 = " , ans, " mg/L "]
Print["Concentration in Raff = " ,
funx[1]*(cxo - SuperStar[cx1]) + SuperStar[cx1], " mg/L "]
Print["Concentration in Loaded Solvent = " ,
funy[0]*(SuperStar - cy1) + cy1, " mg/L "]The 1st and 2nd boundary conditions are:y'[0] == 0, x[0] == 1
POSTED BY: Yong Wang
Posted 10 years ago
Many parts of this worry me.
 In[1]:= Clear[fx, fy, cxo, cy1, h, pey, noofc, hox, ans];
 (*Input data for the calculation of the outlet concentrations*)
 (*Flowrate of Cu Raff (m3/h)*)fx = 5;
 (*Flowrate of barren solvent (m3/hr)*)fy = 1;
 (*Subscript[U,3] Subscript[O,8] concentration in Cu Raff (mg/L)*)cxo = 267;
 (*Subscript[U,3] Subscript[O,8] concentration in barren solvent(mg/L)*)cy1 = 130;
 (*Effective column height (m)*)h = 9;
 (*Peclet number of organic phase*)pey = .7;
 (*no of compartments*)noofc = 13;
(*Height of transfer unit (m)*)hox = 1.65855537601783531;
(*Equilibrium constants represented by Langmuir model,cx^*=(a*cy)/(b*cy+1)*)
a = 0.00982089325801922;
b = -0.0002962283775423278;
(*Press once you complete the input section*)
nox = h/hox;
uxuy = fx/fy; dc = h/noofc; bd = h/dc; SuperStar = (a/cxo - b)^-1; SuperStarcx1 = (a*cy1)/(b*cy1 + 1);
b1 = cxo - SuperStarcx1; b2 = SuperStar - cy1; a1 = a*(SuperStar - cy1); a2 = a*cy1; a3 = b*(SuperStar - cy1);
a4 = (b*cy1) + 1;
(*yo=.3;*)
parta[yo_?NumericQ] := Module[{},
  ini = NDSolve[{-x'[z] - (nox/b1)*(b1*x[z] + SuperStarcx1 - (a1*y[z] + a2)/(a3*y[z] + a4)) == 0,
     y''[z]+pey*bd*y'[z]+((nox*pey*bd*uxuy)/b2)*(b1*x[z]+SuperStarcx1-(a1*y[z]+a2)/(a3*y[z]+a4)) == 0,
     y'[0] == 0, x[0] == 1, y[0] == yo}, {x[z], y[z]}, {z, 0, 1}];
  funx[z_] = ini[[1, 1, 2]];
  funy[z_] = ini[[1, 2, 2]];
  funy'[1] + pey*bd*funy[1]
  ]

In[17]:= lowz = 0.; highz = 1.;
Do[midz = (lowz + highz)/2.0;
  If[Sign[parta[lowz]] == Sign[parta[midz]], lowz = midz, highz = midz], {30}
  ];
ans = midz;

(*try=FindRoot[parta\[Equal]0,{ans,0,1},MaxIterations\[Rule]350,AccuracyGoal\[Rule]6];*)

In[20]:= Plot[{funx[z], funy[z]}, {z, 0, 1}, Frame -> True,
FrameLabel -> {"Dimensionless length (Z) ", "Dimensionless conc (X,Y)"},
PlotStyle -> AbsoluteThickness[2], PlotRange -> All, ImageSize -> 350]

Out[20]= ...PlotSnipped...

In[21]:= Print["y0 = ", ans, " mg/L "]

During evaluation of In[21]:= y0 = 0.458824 mg/L

In[22]:= Print["Concentration in Raff = ", funx[1]*(cxo - SuperStarcx1) + SuperStarcx1, " mg/L "]

During evaluation of In[22]:= Concentration in Raff = 3.36841 mg/L

In[23]:= Print["Concentration in Loaded Solvent = ", funy[0]*(SuperStar - cy1) + cy1, " mg/L "]

During evaluation of In[23]:= Concentration in Loaded Solvent = 1448.16 mg/L

The 1 st and 2 nd boundary conditions are : y'[0] == 0, x[0] == 1

Please test this carefully to find and correct more errors and misunderstandings
POSTED BY: Bill Simpson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract