MATLAB Matlab iteration issue (I think :P)

[roesbak]

Hi there,

I am working on a project in a subject called "Multiphase Flow in Pipes." We are to make a "simulator", describing the pressure (step-wise) in a pipe inserted vertically into water. Air is injected at the bottom of the pipe with a given mass rate, m, and for simplicity, ideal gas, no slip conditions, and a given QL (liquid flow rate) are assumed. The pressure-losses taken into account are hydrostatic and friction only. (I can explain further on request, but my problem is, I hope, code-related, and not physics :P)

I started writing what I thought was a fairly simple piece of "code", and decided to test my program with one given m and QL. However, based on my graph, it looks as though the pressure-drop upwards in the pipe is constant, and does not change, no matter what m and QL I try.

I think I may have done something wrong in my "for" loop, but since I am very inexperienced with Matlab, I thought I would ask if anyone might give me a small pointer..

Code follows:
------------------------------

clf
% First, we set up the main constants
g = 9.81;
h = 50;
D = 0.15;
rhoL = 1000;
P(1) = 1*10^5 + (rhoL*h*g);
Ptop = 1.013*10^5;
Pref = 1.013*10^5;
nu = 1.5*10^-5;
viscL = 0.001;
rhoref = 1.2;


% Then we start entering the equations necessary to iterate.

dH = 0.5;
H = 0:dH:h;
NH = length(H);


m = 1;

ql = 0.001;

for i = 2:NH % Loop, upwards in the well.
% Calculate gas density and friction pressure drop

rhoG = rhoref * (P(i-1)/Pref); % density of gas

qg = (m/rhoG); % volumetric flow rate of gas

EG = (qg)/(qg + ql); % gas fraction

EL = (ql)/(qg + ql); % liquid fraction

rhom = (EL * rhoL) + (EG * rhoG); % density of gas/liquid mixture

viscm = (EL * viscL) + (EG * rhoG * nu); % viscosity of gas/liquid mixture

U = (ql+qg)/((pi/4)*D^2); % flow speed

phyd = rhom * g; % hydrostatic pressure loss

pfrict = (4/D) * 0.046 * ((rhom*U*D)/(viscm))^(-0.2); % frictional pressure loss

ptot = phyd + pfrict; % total pressure loss

P(i) = P(i-1) - (ptot*dH)

end


figure(1)
plot(H,P)
---------------------------------------

Thank you very much for your time.


Yours

André Røsbak

 

[FrogPad]

Can you ask a specific question, I don't really feel like looking through all your code. I'm sure I can help if you ask a question pertaining to MATLAB though.

Maybe someone else will look through it though.

However, since you were talking about for loops. The syntax is

for i=a:h:b
end

a=start value
b=end value
h="step size"

so...

for i=1:1:3
i
end

would print:
1
2
3

and,
for i=1:0.5:2
i
end

would print,
1
1.5
2

 

[roesbak]

Thank you for your reply.

I understand that you don't want to look at all of it, but I thought I should include it at least.

My question/problem is: My loop does not seem to be working. I think the problem is that my variable "pfrict" does not update itself for each step. Have I done something wrong when defining "prfict"? (i.e. does it "update" itself for each step, based on what I've written)

The graph I get is linear, i.e. the amount subtracted from the original pressure is the same, or at least veeeery similar, for each time step. This should not be the case. The results should at the very least vary for different m and QL inputs (defined in the lines just above where the for loop starts), but it does not.

Thank you.


- André

 

[LeBrad]

My question/problem is: My loop does not seem to be working. I think the problem is that my variable "pfrict" does not update itself for each step. Have I done something wrong when defining "prfict"? (i.e. does it "update" itself for each step, based on what I've written)

The graph I get is linear, i.e. the amount subtracted from the original pressure is the same, or at least veeeery similar, for each time step. This should not be the case. The results should at the very least vary for different m and QL inputs (defined in the lines just above where the for loop starts), but it does not.

For debugging purposes you should save all your variables at each time step. E.g.
pfrict(i) = (4/D) * 0.046 * ((rhom(i)*U(i)*D)/(viscm(i)))^(-0.2);
That makes it much easier to trace backwards and find where the problem is coming from.