47 views (last 30 days)

Show older comments

Joseph Improta about 17 hours ago

Edited: Torsten about 2 hours ago

- Castro-Orgaz and Chanson - 2022 - Free surface profiles of near-critical instabiliti.pdf

Open in MATLAB Online

Hello all, I am working on solving a 6 equation ODE to model the shape of a free surface in a channel. The paper that I am referencing only uses 5 euqations but i have no idea how they are solving for one of the variables so I added the variable and equation to the ODE. I have tried almost solver avalible in matlab and am very stuck on what to do next. Ive been in contact with the papers authors but they havent been a huge help. If anyone has a suggestion, it would be a huge help.

Here is the code:

clear; close all;

% Initial conditions and parameters

x0 = 0; %m

xfinal = 1; %m

h0 = 0.11; %m

hx0 = 0;

hxx0 = 0;

s0 = 1/163; % bottom slope

F1 = 1.07; % Froude number

N = 0.142857;

K = 2; % curvature

q = 0.12; % m^2 / s

g = 9.81; % gravity

k = 0;

theta = 0.00009;

beta0 = (1 + N)^2 / (1 + 2 * N);

ep0 = (h0 * hxx0) / (1 + hx0^2);

ep1 = hx0^2 / (1 + hx0^2);

S0 = (h0^2 / 2) + beta0 * (q^2 / (g * h0)) * (1 * ep0 * ((2 / (K + N + 2)) - ((1 - 2 * N) / (K + 2 * N + 2))) - ep1 * ((N + 1) / (N + 3)));

f0 = [h0; hx0; S0; theta; k;0];

% Solve the ODE

options = odeset('RelTol',1e-5,'AbsTol',1e-7);

[x, f] = ode23tb(@undularODE, [x0 xfinal], f0, options);

Warning: Failure at t=3.074452e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.092265e-17) at time t.

% Plot the results

figure;

plot(x, f);

legend('h', 'hx', 'S', 'theta', 'k','N');

xlabel('x');

ylabel('y');

title('Solution');

Here is the UdularODE function:

function dfdx = undularODE(x, f)

% Parameters

g = 9.81;

N = 0.142857;

q = 0.12; % m^2 / s

K = 2;

s0 = 1/163; % bottom slope

utheta = 1.1; % m/s velocity at momentum thickness

theta = f(4);

Nu = 1.12*10^-6; % kinematic viscosity

% Computed parameters

beta = (1 + N)^2 / (1 + 2 * N);

ep1 = f(2)^2 / (1 + f(2)^2); % epsilon 1

Re = 178458; % Reynolds number

Cf = 0.075 / (log(Re) - 2)^2; % skin friction coeff

ff = 4 * Cf * (1 + N)^2; % friction factor

sf = ff / 8 * (1 + 2 * (f(1) / 1)) * 1.47^2; % friction slope

Ue = (1 + N) * (q / f(1)); % max velocity at boundary layer edge

k = utheta / Ue;

H = (1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(2/3); % shape factor

Rtheta = (Ue * theta) / Nu; % momentum thickness Reynolds number

N2 = (H-1)/2; % calculated value of N

dUedx = -1*Ue*f(6)*(1+N2)*(Ue/f(1))*f(2);%dUedx

gamma = (theta / Ue) * dUedx * Rtheta^(1/4); % Buri Shape Factor

dHdx = (3*k^2 - 5.5*k + 3.68)^(2/3) * ( (1.46*(1-k^2))/( (1.5+k)*theta*Rtheta^(1/4)* (gamma +0.118*(0.67-k)) ) );

%dNdx(end+1)= 0.5*dHdx;

% Define the system of ODEs

hx = f(2);

eq2 = (((g * f(1)) / (beta * q^2)) * (f(3) - (f(1)^2 / 2)) - 1 + ep1 * ((N + 1) / (N + 3))) / ...

((f(1) / (1 + f(2)^2)) * (2 / (K + N + 2)) - ((1 - 2 * N) / (K + 2 * N + 2)));

eq3 = f(1) * (s0 - sf);

eq4 = -1 * f(4)*((1 / (2 + 2 * N)) * (dHdx) - (1 / f(1)) * hx) * (H + 2) + Cf / 2;

eq5 = 1.46 * ((1 - f(5)^2) / ((1.5 + f(5)) * theta * Rtheta^(1/4))) * (gamma + 0.118 * (0.67 - f(5)));

eq6 = 0.5*dHdx;

% Return the derivatives

dfdx = [hx;

eq2;

eq3;

eq4;

eq5;

eq6];

%keyboard

end

##### 9 Comments Show 7 older commentsHide 7 older comments

Show 7 older commentsHide 7 older comments

Torsten about 13 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3205936

Edited: Torsten about 13 hours ago

Open in MATLAB Online

Why do you use

k = utheta / Ue;

and not

k = f(5)

? That's inconsistent.

Further if

H = (1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(2/3)

, then

dHdx = -2/3*(1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(-1/3) * (1.3 + 3*2*(0.7 - k)) * dk/dx

Does this equal your dHdx ?

Sam Chak about 13 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3205966

Open in MATLAB Online

On top of that, there is another inconsistency:

dUedx = -1*Ue*f(6)*(1 + N2)*(Ue/f(1))*f(2);

but in the paper

Sam Chak about 12 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3205981

Also, since H is a function of k in Eq. (15), then

.

The 6th state equation is probably unnecessary because

can be directly derived

and, will also be a function of k.

Because both depend on the state k, only is necessary, which is the 5th state equation in Eq. (13):

I understand why you coded .

Sam Chak about 12 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3206001

⋮

Edited: Sam Chak about 11 hours ago

Open in MATLAB Online

Also, below Eq. (10), you can see how U is defined. This U is needed in Eq. (16), but you coded it differently.

Note that depends on the solution of , and depends on both solutions of and , and in Eq. (11) requires the solution of .

(Edited to correct mistakes pointed out by @Torsten).

Since the authors didn't make as the 6th ODE, most likely they used the approximation of .

Try fixing the code below:

function dfdx = undularODE(x, f)

% definitions

h = f(1);

hx = f(2);

S = f(3);

theta = f(4);

k = f(5);

% Independent parameters

g = 9.81;

% N = 0.142857; % Probably not needed?

q = 0.12; % m^2 / s

K = 2;

s0 = 1/163; % bottom slope

utheta = 1.1; % m/s velocity at momentum thickness

Nu = 1.12*10^-6; % kinematic viscosity

Re = 178458; % Reynolds number

b = 1;

F = 1.47;

% Dependent parameters

ep1 = hx^2/(1 + hx^2); % epsilon 1

Cf = 0.075/(log(Re) - 2)^2; % skin friction coeff

ff = 4*Cf*(1 + N)^2; % friction factor

sf = (ff/8)*(1 + 2*h/b)*F^2; % friction slope

% k = utheta / Ue;

H = (1.3 + 1.3*(0.7 - k) + 3*(0.7 - k)^2)^(2/3); % shape factor

N = (H - 1)/2; % calculated value of N

beta = (1 + N)^2/(1 + 2*N);

Ue = (1 + N)*(q/h); % max velocity at boundary layer edge

Rtheta = (Ue * theta) / Nu; % momentum thickness Reynolds number

dHdx = ...; % depends on the solution of k(x)

dUedx = ...; % depends on the solutions of k(x) and h(x)

gamma = (theta / Ue) * dUedx * Rtheta^(1/4); % Buri Shape Factor

% dNdx(end+1)= 0.5*dHdx;

num2 = (g*h)/(beta*q^2)*(S - h^2/2) - 1 + ep1*(N + 1)/(N + 3);

den2 = h/(1 + hx^2)*(2/(K + N + 2) - (1 - 2*N)/(K + 2*N + 2));

% Define the system of ODEs

dhdx = hx;

dhxdx = num2/den2;

dSdx = h*(s0 - sf);

dtheta = - theta*(dHdx/(2 + 2*N) - hx/h)*(H + 2) + Cf/2;

dkdx = 1.46*(1 - k^2)/((1.5 + k)*theta*Rtheta^(1/4))*(gamma + 0.118*(0.67 - k));

% eq6 = 0.5*dHdx;

% Return the derivatives

dfdx = [dhdx

dhxdx

dSdx

dtheta

dkdx];

end

Torsten about 12 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3206011

Open in MATLAB Online

N = (H-1)/2

instead of

N = 0.142857;

And I don't think that derivatives depend on initial values:

dHdx = ...; % depends on the initial value of k

dUedx = ...; % depends on the initial values of k and h

Sam Chak about 12 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3206026

@Torsten is right. I incorrectly phrased the comments. When typing those comments, my mind was thinking whether there is an algebraic loop dependency between k and because .

Sam Chak about 11 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3206056

Please review the computation of Cf in the code against the formula in Eq. (14):

Joseph Improta about 3 hours ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3206246

@Torsten @Sam Chak thank you for your suggestions. Im gonna go though them and work on myt code. I will get back to you with what I find. I will say this is how my mind interpreted it initally: to calculate dH/dx. I am doing the chain rule with dH/dk and dk/dx to obtain dH/dx. To calculate dk/dx though, I need gamma which needs dUe/dx which then needs dN/dx. The only way I have found to calculate dN/dx is to do the chain rule with dN/dH and dH/dx.

Torsten 2 minutes ago

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2135551-unstable-6-equation-ode#comment_3206526

Edited: Torsten less than a minute ago

Open in MATLAB Online

Insert the expression for gamma in equation (13). Then use

dN/dx = 0.5*dH/dx = 0.5*(-2/3*(1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(-1/3) * (1.3 + 3*2*(0.7 - k)) * dk/dx)

and insert it in the expression for dUe/dx.

You get an equation of the form

dk/dx = U*0.5*(-2/3*(1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(-1/3) * (1.3 + 3*2*(0.7 - k)) * dk/dx) - (1+N)*U/h * dh/dx

Insert f(2) for dh/dx and solve for dk/dx.

Sign in to comment.

Sign in to answer this question.

### Answers (0)

Sign in to answer this question.

### See Also

### Categories

SciencesPhysics

Find more on **Particle & Nuclear Physics** in Help Center and File Exchange

### Tags

- fluid
- ode
- ode45
- unstable

### Products

- MATLAB

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français

- United Kingdom(English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)

Contact your local office