shooting method solving compressible boundary layer equations

Hi, i hope someone can help me. I want to find the solution to the compressible boundary layer equations, this problem is part of my thesis project, but I'm running into some problems.

here is the ODEs

here is the boundary conditions

here is the matlab code

function [x,y] = shooting
% Use fsolve to ensure the boundary function is zero. The result is the
% unknown initial condition.
opt = optimset('Display','off','TolFun',1E-20);
F = fsolve(@(F) eval_boundary(F),[0.0826,0,0,0,25.8],opt);
% Solve the ODE-IVP with the converged initial condition
[x,y] = solve_ode(F);
end
function [x,y] = solve_ode(F)
%basic value of outer flow
Te = 218.6;   mach = 8.0;   Pr = 0.75; 
he = 1004.68 * Te;    ue = mach * sqrt(1.4 * 8.29586 * Te / 0.02885); 
ratio1 = ue * ue / he;    cm = 110.4 / Te;  ck = 110.4 / Te;

% C = sqrt(y(5))*(1.0+ratio2)/(y(5)+ratio2);
% Solve the ODE-IVP with initial condition F on [0 100] (arbitrary upper bound)
[x,y] = ode45(@(x,y) [  -y(3) * y(1) / ( sqrt(y(5)) * (1.0 + cm) / ( y(5) + cm ) );
                        y(1) /  ( sqrt(y(5)) * (1.0 + cm) / ( y(5) + cm ) );
                        y(2);
                        -Pr * ( y(3) * y(4) /  ( sqrt(y(5)) * (1.0 + ck) / ( y(5) + ck ) ) + ratio1 * y(1) * y(1) /  ( sqrt(y(5)) * (1.0 + cm) / ( y(5) + cm ) ) ) ;
                        y(4) /  ( sqrt(y(5)) * (1.0 + ck) / ( y(5) + ck ) )    ] , [0 100] , F); %solve BVP                
end
function [g] = eval_boundary(F)
% Get the solution to the ODE with inital condition F
[x,y] = solve_ode(F);
% Get the function values (for BCs) at the starting/end points
y2_start = y(1,2); %f(0) = 0
y3_start = y(1,3); %f'(0) = 0
y4_start = y(1,4); %g'(0) = 0
y1_end = y(end,1); %f''(end) = 0
y2_end = y(end,2); %f'(end) = 1
y4_end = y(end,4); %g'(end) = 0
y5_end = y(end,5); %g(0) = 1
% Evaluate the boundary function
g = [
        y2_start
        y3_start
        y4_start
%         y1_end
        y2_end-1
%         y4_end
        y5_end-1
        ];
end

this code is come from Mohammad A Alkhadra (https://ww2.mathworks.cn/matlabcentral/fileexchange/69310-solving-blasius-equation-with-the-shooting-method) for solving blasius boundary layer equations. I changed the equations and BCs but it can not get correct answers (the code return complex number which must be wrong).

I have checked many times for the code and equations deduction, but did not find the wrong place.

I don't know if it is because these ODEs are differential equation with variable coefficients ？ this is one different point from the original code.

[darova]

Equations for image and code don't look similar

First equation can be simplified to second order?

[zhou liu]

Thanks for your quick comment!

the equations deduction is come from https://flylib.com/books/en/4.201.1.212/1/

The author also offered a JAVA code using his own ODE solver.

package TechJava.MathLib;
import TechJava.Gas.*;
public class CompressODE extends ODE
{
	double he, ue, Te, mach, Pr, C;
	double ratio1, ratio2;
	
	// The CompressODE constructor calls the ODE
	// constructor passing it some compressible
	// boundary layer specific values. There are five
	// first order ODEs and two free variables.
	
	public CompressODE(double Te, double mach) {
		super(5,2);
		this.Te = Te;
		this.mach = mach;
		he = 3.5*AbstractGas.R*Te/0.02885;
		ue = mach*Math.sqrt(1.4*AbstractGas.R*Te/0.02885);
		ratio1 = ue*ue/he;
		ratio2 = 110.4/Te;
		Pr = 0.75;
	}
	
	// The getFunction() method returns the right-hand
	// sides of the five first-order compressible
	// boundary layer ODEs
	// y[0] = delta(C*f'') = delta(n)*(-f*f'')
	// y[1] = delta(f') = delta(n)*(f'')
	// y[2] = delta(f) = delta(n)*(f')
	// y[3] = delta(Cg') = delta(n)*(-Pr*f*g' –
	// Pr*C*(ue*ue/he)*f''*f'')
	// y[4] = delta(g) = delta(n)*(g')
	
	public void getFunction(double x, double dy[],
	                        double ytmp[]) {
		C = Math.sqrt(ytmp[4])*(1.0+ratio2)/
		(ytmp[4]+ratio2);
		dy[0] = -ytmp[2]*ytmp[0]/C;
		dy[1] = ytmp[0]/C;
		dy[2] = ytmp[1];
		dy[3] = -Pr*(ytmp[2]*ytmp[3]/C +
		ratio1*ytmp[0]*ytmp[0]/C);
		dy[4] = ytmp[3]/C;
	}
	
	// The getE() method returns the error, E[], in
	// the free variables at the end of the range
	// that was integrated.
	
	public void getError(double E[], double endY[]) {
		E[0] = endY[1] - 1.0;
		E[1] = endY[4] - 1.0;
	}
	
	// This method initializes the dependent variables
	// at the start of the integration range. The V[]
	// contains the current guess of the free variable
	// values
	
	public void setInitialConditions(double V[]) {
		setOneY(0, 0, V[0]);
		setOneY(0, 1, 0.0);
		setOneY(0, 2, 0.0);
		setOneY(0, 3, 0.0);
		setOneY(0, 4, V[1]) ;
		setOneX(0, 0.0);
	}
}

import TechJava.MathLib.*;
public class ShootingCompress
{
	public static void main(String args[]) {
		
		// Create a CompressODE object
		
		CompressODE ode = new CompressODE(218.6, 8.0);
		
		// Solve the ODE over the desired range with the
		// specified tolerance and initial step size.
		// The V[] array holds the initial conditions of
		// the free variables.
		
		double dx = 0.1;
		double range = 5.0;
		double tolerance = 1.0e-6;
		double V[] = {0.0826, 25.8};
		
		// Solve the ODE over the desired range
		
		int numSteps =
			ODESolver.ODEshooter(ode, V, range, dx, tolerance);
			
		// Print out the results
		
		System.out.println("i eta Cf'' f' f");
		for(int i=0; i<numSteps; ++i) {
			System.out.println(
				""+i+" "+ode.getOneX(i)+" "+ode.getOneY(i,0)+
				" "+ode.getOneY(i,1)+" "+ode.getOneY(i,2));
		}
		System.out.println();
		System.out.println("i eta Cg' g");
		for(int i=0; i<numSteps; ++i) {
			System.out.println(
				""+i+" "+ode.getOneX(i)+
				" "+ode.getOneY(i,3)+" "+ode.getOneY(i,4));
		}
	}
}
 
 

The equations and BS settings are same as that in his code, but I cann't get the correct answer using MATLAB.

Answer 1

[darova]

bvp4c should work

see attached script

[zhou liu]

Thank you soooooo much for your generous help！

I saw your MATLAB code and found that I omitted a key condition so mislead you. So sorry!

In fact, the parameter C is a function of g.

According to the reference,

So the discription of the function should be different.

I put my function into your code as follow (equation order is different):

function bvp4c_mathworks
nspan = [0 1];
init = zeros(1,5);
solinit = bvpinit(nspan,init);
sol = bvp4c(@ode4,@bc4,solinit);
nx = sol.x;
% f = sol.y(1,:);
% df  = sol.y(2,:);
% d2f = sol.y(3,:);
% g = sol.y(4,:);
% dg= sol.y(5,:);
plot(sol.y(2,:),nx)
xlabel('\it{ f^{(1)}}');
ylabel('\eta');
% hold on
% plot(nx,f,'r')
% hold off
%legend('f','df','d2f','g','dg')
end

function du = ode4(n,u)
%basic value of outer flow
Te = 218.6;
mach = 8.0;   
% Pr = 0.75; 
he = 1004.68 * Te;    
ue = mach * sqrt(1.4 * 8.29586 * Te / 0.02885); 
ratio1 = ue * ue / he;    
cm = 110.4 / Te;  
% ck = 110.4 / Te;

du = zeros(5,1);

du(1) = -u(3) * u(1) / ( sqrt(u(5)) * (1.0 + 0.5050) / ( u(5) + 0.5050 ) );                       % d3f
du(2) =  u(1) /  ( sqrt(u(5)) * (1.0 + 0.5050) / ( u(5) + 0.5050 ) );                      % d2f
du(3) = u(2);             % df
du(4) =  -0.75 * ( u(3) * u(4) /  ( sqrt(u(5)) * (1.0 + 0.5050) / ( u(5) + 0.5050 ) ) + 25.6446 * u(1) * u(1) /  ( sqrt(u(5)) * (1.0 + 0.5050) / ( u(5) + 0.5050 ) ) ) ;                       % d2g
du(5) = u(4) /  ( sqrt(u(5)) * (1.0 + 0.5050) / ( u(5) + 0.5050 ) );  % dg
end

function res = bc4(u0, un)
res = [u0(2)        
       u0(3)        
       u0(4)        
       un(2)-1      
       un(5)-1];    
end

it shows:

>> bvp4c_mathworks2

Error using bvp4c (line 251)

Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in bvp4c_mathworks2 (line 5)

sol = bvp4c(@ode4,@bc4,solinit);

by the way, according to the reference, the result should be

So if the nspan should be large?

Thank you again!!

[darova]

I saw your MATLAB code and found that I omitted a key condition so mislead you. So sorry!

oh no! all my work is pointless!

[zhou liu]

I‘m sorry.....

It's so nice of you. Thanks very much~

You proposed a different useful frame for solving this problem at least, haha:）

I will continuously update my work on it untill the correct answer.

[darova]

But can i have one more try?

I differentiated C function and it looks like below

According to this origin equations can be re-written as:

Unable to solve the collocation equations -- a singular Jacobian encountered.

I just changed initial conditions: all ones. See attached script. This time it should be cool!

[zhou liu]

Thanks again!

You have different function structure from mine. I define Cf'' as a function of eta, but you define f'' as a function. So both of us are correct in math, right?

In addition, I tried your code, and indeed the five parameters(f, f', f'', g, g' ) are consistent with the boundary values. But it seems not consistent with the correct answer as follows:

PS：

I I differentiated C function and found:

which is different from you :(

However, both of our equations can not get the correct answer, and I found that the "nspan" is quite sensitive. When I increase the span number from [0 1] to [0 25], it will get "a singular Jacobian encountered". Sad...

[darova]

Indeed, you differentiated C function in a correct way

The shape of results looks similar

What about Y axis? Maybe it's about scaling?

[Mouli Bhaskar Duddupudi]

I have been looking into this problem. What could be reason for Y-axis? Is it because of scaling?

Answer 2

[Mohamad]

See below:

https://www.mathworks.com/matlabcentral/fileexchange/73587-compressible-similarity-solution-for-flat-plate-bl