Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab iteration issue (I think :P)

  1. Feb 8, 2008 #1
    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:

    % 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)



    Thank you very much for your time.

    Yours Sincerely,

    André Røsbak
    Last edited: Feb 8, 2008
  2. jcsd
  3. Feb 8, 2008 #2
    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

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


    for i=1:1:3

    would print:

    for i=1:0.5:2

    would print,
  4. Feb 8, 2008 #3
    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é
  5. Feb 10, 2008 #4
    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook