Unstable 6 equation ode (2024)

47 views (last 30 days)

Show older comments

Joseph Improta about 17 hours ago

  • Link

    Direct link to this question

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

  • Link

    Direct link to this question

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

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');

Unstable 6 equation ode (2)

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

Torsten about 13 hours ago

Direct link to this comment

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

  • Link

    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

  • Link

    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

Unstable 6 equation ode (5)

Sam Chak about 12 hours ago

Direct link to this comment

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

  • Link

    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

Unstable 6 equation ode (7).

The 6th state equation is probably unnecessary because

Unstable 6 equation ode (8) can be directly derived

andUnstable 6 equation ode (9), will also be a function of k.

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

Unstable 6 equation ode (11)

I understand why you coded Unstable 6 equation ode (12).

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.

Unstable 6 equation ode (14)

Note that Unstable 6 equation ode (15) depends on the solution of Unstable 6 equation ode (16), and Unstable 6 equation ode (17) depends on both solutions of Unstable 6 equation ode (18) and Unstable 6 equation ode (19), and Unstable 6 equation ode (20) in Eq. (11) requires the solution of Unstable 6 equation ode (21).

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

Unstable 6 equation ode (22)

Since the authors didn't make Unstable 6 equation ode (23) as the 6th ODE, most likely they used the approximation of Unstable 6 equation ode (24).

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

  • Link

    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

  • Link

    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 Unstable 6 equation ode (27) because Unstable 6 equation ode (28).

Sam Chak about 11 hours ago

Direct link to this comment

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

  • Link

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

Unstable 6 equation ode (30)

Joseph Improta about 3 hours ago

Direct link to this comment

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

  • Link

    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

  • Link

    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.


Unstable 6 equation ode (33)

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

Unstable 6 equation ode (2024)
Top Articles
Latest Posts
Article information

Author: Cheryll Lueilwitz

Last Updated:

Views: 6642

Rating: 4.3 / 5 (74 voted)

Reviews: 89% of readers found this page helpful

Author information

Name: Cheryll Lueilwitz

Birthday: 1997-12-23

Address: 4653 O'Kon Hill, Lake Juanstad, AR 65469

Phone: +494124489301

Job: Marketing Representative

Hobby: Reading, Ice skating, Foraging, BASE jumping, Hiking, Skateboarding, Kayaking

Introduction: My name is Cheryll Lueilwitz, I am a sparkling, clean, super, lucky, joyous, outstanding, lucky person who loves writing and wants to share my knowledge and understanding with you.