shooting method solving compressible boundary layer equations

28 Ansichten (letzte 30 Tage)
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.
  2 Kommentare
darova
darova am 13 Okt. 2019
Equations for image and code don't look similar
First equation can be simplified to second order?
zhou liu
zhou liu am 15 Okt. 2019
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.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

darova
darova am 15 Okt. 2019
bvp4c should work
see attached script
  7 Kommentare
darova
darova am 16 Okt. 2019
Indeed, you differentiated C function in a correct way
The shape of results looks similar
img1.png
What about Y axis? Maybe it's about scaling?
123.png
Mouli Bhaskar Duddupudi
Mouli Bhaskar Duddupudi am 28 Feb. 2022
Bearbeitet: Mouli Bhaskar Duddupudi am 28 Feb. 2022
I have been looking into this problem. What could be reason for Y-axis? Is it because of scaling?

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Mohamad
Mohamad am 11 Dez. 2019

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by