Shooting method - Matlab ODE

Question

I'm trying to solve these equations of hypersonic adiabatic flow over a flat plate. I did all the simplifications and got these equations for the stagnation point flow. $$\left(Cf''\right)' + f f'' = \left(f'\right)^2-g$$ $$\left(\frac{C}{p_r} g'\right) + fg'=0$$

For a calorically perfect gas, $g=h/he$, and $he = \rm stagnation ~enthalpy$

$p_r$ is a constant and $C=g^{-1/3}$

$f(0)=0$ (over the streamline of the body)

$f'(0)=0$ (this is $u$, the velocity over the body which is null)

$f''(0)=?$ white suggests to try as a first value $0.664/\sqrt{C_w}$

I'm doing this for the adiabatic case so

$g(0)=g_{wall}$

$g'(0)=0$ (no thermic flux)

UPDATE

Control conditions:
$f'(\delta)=1$ and $g(\delta)=1$ I want to integrate from $\delta=0$ to $\delta=6$, in the end, if $f'(\delta)$ is different than $1$ I'll have to change the value assumed for $f''(0)$.

How can I solve these equations, someone suggested me using the shooting method but how could I code this in matlab? I would like just a help on how to implement the equations on matlab, not a final solution.

Answer 1

To change this into a system of first order ODEs, you can define

$$Y(x;s) := \left( \begin{split} f''\\ f'\\ f\\ g'\\ g \end{split} \right) = \left(\begin{split} Y_1\\ Y_2\\ Y_3\\ Y_4\\ Y_5 \end{split} \right)\,.$$ Your first equation then gives $$ CY_1' + C'Y_1 + Y_3Y_5= (Y_2)^2 - Y_5\,,$$ or, with $C=g^{-1/3}$, $C' = -\frac13g^{-4/3}g'$, $$ Y_5^{-1/3}Y_1' - \tfrac13Y_5^{-4/3} Y_4Y_1 + Y_3Y_5= (Y_2)^2 - Y_5\,,$$ and the second differentiated is $$\left(\frac{Y_5^{-1/3}}{p_r}+Y_3\right)Y_4' + \left(\frac{-\frac13Y_4Y_5^{-4/3}}{p_r}+Y_2\right)Y_4 = 0\,,$$ with the additional equations $$Y_2' = Y_1, \quad Y_3' = Y_2, \quad Y_5' = Y_4\,.$$

Bringing these together, gives $$ Y'(x;s) = \left(\begin{split} Y_5^{1/3}\left(\tfrac13Y_5^{-4/3} Y_4Y_1 - Y_3Y_5 + (Y_2)^2 - Y_5\right)\\ Y_1\\ Y_2\\ - \left(\frac{Y_5^{-1/3}}{p_r}+Y_3\right)^{-1} \left(\frac{-\frac13Y_4Y_5^{-4/3}}{p_r}+Y_2\right)Y_4\\ Y_4 \end{split}\right) $$ and initial conditions $$Y(0;s) := \left( \begin{split} s\\ 0\\ 0\\ 0\\ g_{\text{wall}} \end{split} \right) \,.$$

You can use an ODE solver to integrate this the find $Y(\delta)$ for some $\delta>0$, and change $s$ to match your additional condition $$F(s) := Y_2(\delta;s) - 1 = 0\,,$$ or $$F(s) := Y_5(\delta;s) - 1 = 0\,.$$

This could be by bisection (find 2 values $s_1$ and $s_2$ with $F(s_1)<0$, $F(s_2)>0$, and then home in by testing $\tilde{s}=0.5(s_1+s_2)$ and replace $s_1$ or $s_2$ as appropriate) or even Newton's method (more work is required to get $F_s(s)$, but it's possible).