function hydrodynamics1

clear all
clc

% Modelling hydrodynamics
global ReL GaL eps epsL epsL0 A B h1 h2 SL0 g FL sigmaL rohL muL uL Temp de Eostar
global ReG GaG epsG FG rohG muG uG DR L Ptot VpointL mpointG 
% Temperatur [K]
Temp=293.15;
% Pressure [MPa]
Ptot=0.2;
% Fitting parameters A B h1 h2 with 40 ppi foam (from stemmet page 24)
A=8800; 
B=5.04; 
h1=4.25; 
h2=1.88;
% strut length [m]
de = 0.38e-3;
% dynamic viscosity of nitrogen (approximation for air) at 293.15 K and 0.2 MPa [uPa*s]
% (Source: http://webbook.nist.gov/chemistry/fluid/)
muG=17.596;
% dynamic viscostiy of water at 293.15 Kand 0.2 MPa [uPa*s] &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% (Source: http://webbook.nist.gov/chemistry/fluid/)
muL=1001.6;
% acceleration due to gravity [m^2/s]
g= 9.8;
% Bed porosity; (Source: maximum of I. Mohammed presentation, 2012: 
%Hydrodynamic investigation of a flow reactor with foam packings)
eps=0.94;  %??
% static holdup of 40 ppi foam [-] (Source: Stemmet 2008, page 21
%epsL0=0.115; not used
% Density of nitrogen (approximation for air) in 293.15 K 0.2 MPa[kg/m3]
% (Source: http://webbook.nist.gov/chemistry/fluid/)
rohG=2.2997;
% Density of Water in 293.15 K 0.2 MPa[kg/m3]&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&$$
rohL=998.25; 
% Reactor Diameter [m]
DR=0.5;%??
% Reactor length [m]
L=1.58;%??
% Liquid Flowrate [l/min]
%?? VpointL= 5;
% Liquid superficial velocity [m/s] 
uL=0.001415  %??(VpointL/1000/60)/(3.14*DR^2/4);
% Reynolds number of the Liquid ??
ReL=62.22731 %?? rohL*uL*de/(muL*(1-eps));
% Galileo number of the Liquid ??
GaL=rohL^2*g*de^3*eps^3/(muL^2*(1-eps)^3);
% ?? Hydrogen Flowrate [kg/m^2*s]
%?? mpointG=1;
% Hydrogen superficial velocity [m/s]
uG=0.1 %?? mpointG/rohG;
% Reynolds number of Hydrogen ??
ReG=224.1606 %?? rohG*uG*de/(muG*(1-eps));
% Galileo number of Hydrogen  ??
GaG=rohG^2*g*de^3*eps^3/(muG^2*(1-eps)^3);
% Liquid surface tension at 293.15 K
% [N/m](Source:http://de.wikipedia.org/wiki/Oberflaechenspannung)
sigmaL=0.07275;
% E?tv?s number 
Eostar=rohL*g*de^2*eps/(sigmaL*(1-eps)^2);
% saturation of the Liquid at the static holdup SL0=epsL/eps
SL0=1/(eps*(20+0.9*Eostar));

x0=[1.56608e+15 0.894955 7.88252e+13 0.045045]';
opts=optimset('MaxFunEvals',4000);
%opts=optimset('Display','iter','MaxFunEvals',4000);
x1 = fsolve(@NonlinEqs, x0, opts)
  function f=NonlinEqs(x)
  % x(1) = FG; x(2) = epsG; x(3) = FL; x(4) = epsL;
  f(1) = (1/(x(2)/eps)^h1)*(A*ReG/GaG+B*ReG^2/GaG)*rohG*g-x(1)/x(2);
  f(2) = (1/(((x(4)/eps)-SL0)/(1-SL0))^h2)*(A*ReL/GaL+B*ReL^2/GaL)*rohL*g-x(3)/x(4);
  f(3) = x(1)/x(2)-x(3)/x(4)+(rohG-rohL)*g;
  f(4) = eps-x(2)-x(4);

  end
end