Prev: skip bad files in batch processing loop?
Next: Interpolating a point
From: Jay on 19 Jul 2010 13:57
I am trying to model a converging/diverging sonic nozzle, but having problems solving for roots of an equation. Initially I set up the code to solve for Mach number step by step in a 'for' loop, trying to use 'fzero' to find the root (Mach number) and then calculate dM/dx using a backwards difference based on the Mach number that MATLAB solved for. However, solution becomes difficult for MATLAB's solver to handle as the Mach number approaches 1.

I have calculated area (A), dA/dx, and diameter (Diam) for each point in the nozzle, am assuming that k (gamma of gas) and f (friction) = constant.

The equation that I am working with is (where m = Mach number that I'm trying to solve for and n is a counter to specify position down the nozzle):

f = @(m) (((1-m^2)/m)/(1+.5*(k-1)*m^2))*dMdx(n) + dA(n)/(A(n)*dx) - (1/2)*k*m^2*(fric/Diam(n));

When the above equation is solved for zero, this should give the instantaneous mach number of gas flowing through the nozzle.

Since I have been having so much trouble with the iterative solution, now I'm trying to use one of MATLAB's ordinary differential equation solvers, but having trouble with syntax.

Here are my inputs:
A = 40x1 matrix of area at each point
dA = 40x1 matrix of change in area at each point
Diam = 40x1 matrix of diameter of nozzle at each point
dx = constant; step size as I move down the nozzle
k = constant; gamma value for gas
f = constant; friction factor inside nozzle
n = counter number to keep track of position in nozzle

I have created a separate .m file called dMdx to solve for Mach number and change in Mach number when called with one of MATLAB's ode solvers which has the equation re-arranged in the form dM/dx = ... (in the following, f is dM/dx)

f = (1+(1/2)*(k-1)*m^2)*((1/2)*k*m^2*fric/Diam(n) - (1/A(n))*(dA(n)/dx))*m / (1-m^2);

I'd like to get this to output a matrix of Mach numbers for each point as specified by the counter "n" (or one each step (dx) down the nozzle).

Any help for imputing these matrices into the correct syntax for the ODE solver would be greatly appreciated - I've been stuck on this for a while.
From: Torsten Hennig on 20 Jul 2010 02:28
> I am trying to model a converging/diverging sonic
> nozzle, but having problems solving for roots of an
> equation. Initially I set up the code to solve for
> Mach number step by step in a 'for' loop, trying to
> use 'fzero' to find the root (Mach number) and then
> calculate dM/dx using a backwards difference based on
> the Mach number that MATLAB solved for. However,
> solution becomes difficult for MATLAB's solver to
> handle as the Mach number approaches 1.
>
> I have calculated area (A), dA/dx, and diameter
> (Diam) for each point in the nozzle, am assuming that
> k (gamma of gas) and f (friction) = constant.
>
> The equation that I am working with is (where m =
> Mach number that I'm trying to solve for and n is a
> counter to specify position down the nozzle):
>
> f = @(m) (((1-m^2)/m)/(1+.5*(k-1)*m^2))*dMdx(n) +
> dA(n)/(A(n)*dx) - (1/2)*k*m^2*(fric/Diam(n));
>
> When the above equation is solved for zero, this
> should give the instantaneous mach number of gas
> flowing through the nozzle.
>
> Since I have been having so much trouble with the
> iterative solution, now I'm trying to use one of
> MATLAB's ordinary differential equation solvers, but
> having trouble with syntax.
>
> Here are my inputs:
> A = 40x1 matrix of area at each point
> dA = 40x1 matrix of change in area at each point
> Diam = 40x1 matrix of diameter of nozzle at each
> point
> dx = constant; step size as I move down the nozzle
> k = constant; gamma value for gas
> f = constant; friction factor inside nozzle
> n = counter number to keep track of position in
> nozzle
>
> I have created a separate .m file called dMdx to
> solve for Mach number and change in Mach number when
> called with one of MATLAB's ode solvers which has the
> equation re-arranged in the form dM/dx = ... (in the
> following, f is dM/dx)
>
> f = (1+(1/2)*(k-1)*m^2)*((1/2)*k*m^2*fric/Diam(n) -
> (1/A(n))*(dA(n)/dx))*m / (1-m^2);
>
> I'd like to get this to output a matrix of Mach
> numbers for each point as specified by the counter
> "n" (or one each step (dx) down the nozzle).
>
> Any help for imputing these matrices into the correct
> syntax for the ODE solver would be greatly
> appreciated - I've been stuck on this for a while.

Why did you already prescribe the output points in
x-direction ?
If you have an ordinary differential equation
dM/dx = f(M,x) with initial conditon M(x=0) = M_0,
the ODE integrator will solve it for you over [0;L]
where L is the nozzle length.
To get output at specific x-values, you can use
a special MATLAB routine for interpolation after
the integration has finished.
But maybe I misunderstood your question ...

Best wishes
Torsten.
 | 
Pages: 1
Prev: skip bad files in batch processing loop?
Next: Interpolating a point