[MUSIC]
In this video, I'll review the MATLAB learner template for
calculating an unsteady solution.
Let's get started.
I want to call the function flow around cylinder unsteady.
The output will be omega, which is the scale of vorticity and
the computation will be for a Reynolds number of 60.
First, we define the grid, will use 101 grid points in the xi
direction and 202 grid points in the theta direction.
The number of grid intervals in xi is n minus 1,
the number of grid intervals in theta is m minus 2 because
there are two ghost points, there will be a ghost point at minus h and
a ghost point at 2 pi, 2 pi is just a repeat of 0.
The grid spacing is based on the number of grid intervals in theta,
since theta goes from 0 to 2 pi, the grid spacing here is 2 pi divided by capital M.
We then set up the grid for xi starting at 0 and ending at big N times h.
And we set up the grid for theta starting at minus h and ending at 2 pi.
This is a time evolution code, so we need time stepping parameters.
For the greater I'll only ask you to integrate
the odes from t equals 0 to t equals 0.5.
The vortex street starts around t equals 1000.
Before you get to the street, you pass through what looks like a steady state
solution, but then the steady state solution becomes unstable.
So if your computer is not that fast,
you may have to run the code overnight to see the street.
We put the starting and ending times into the variable tspan.
We initialize the scale of vorticity field to 0.
The field is an n by m matrix.
Here I want you to construct the matrix A that shows up in the psi equation.
You'll have to construct the matrix A that takes into account both
the interior grid points and the boundary grid points.
The matrix A is independent of time, so you need to only construct it once.
Then after you construct a you can compute it's L U decomposition which you
can then use to solve the psi equation, that will speed up your computation.
In addition, in this place of the code you can compute and
store any time-independent constants like exponential factors
that show up in the equations and that don't depend explicitly on the time step.
So you only need to compute them once and can save some computational time.
Now we advance the solution using ode23.
The relative tolerance is set to 10 to the minus 3,
that's actually the default tolerance but I write it here explicitly.
We will only time advance the interior values of omega and
not the values on the boundaries.
So here I replaced the omega that contains all the grid values by
an omega that contains just the values at the interior points.
Basically we are removing the values on the four sides of our
computational rectangle.
We then turn omega into a column vector using the column operator.
The output from ode23 will be t and omega,
t is a vector containing the integration times and
omega is a matrix containing the values of omega at each integration time.
The call to ode23 contains the function omega_eq.
This function specifies the right hand sides of the differential equations and
we're going to need to write it.
The function omega_eq has as its arguments omega,
L and U from the L U decomposition of the matrix A and
whatever additional arguments you want,
including the factors that don't depend on time.
And then the other arguments of ode23 include tspan
which we've already assigned, omega containing
the initial conditions of the vorticity and options.
After ode23 finishes running, I want to put the boundary
values back into omega, so I wrote a short code to do this.
I first defined a temporary variable called temp to be a zero matrix
of size n by m, and then I put omega into the interior values of temp.
We need the omega at the last time step and
I make use of the n key word to do this and
we also need to reshape the column vector omega back into a matrix.
Finally, I reassign omega to equal temp and
remove the temp variable.
To finish computing omega we still need its boundary values but
the boundary values on the cylinder depend on the stream function.
So in this bit of code you need to re-compute the stream function,
once you have the stream function you can set the boundary conditions.
Finally, I wrote a code that will plot the vorticity for you.
This is the function plot_Re60 and
that's the end of the main body of the code.
The remaining task is to write the function for
the system of differential equations, you'll do that here.
The function returns the vector d_omega_dt,
that's the time derivatives in the omega equation.
The arguments of your function are omega, capital L,
capital U and whatever additional arguments you add.
So what you're passing to this function is just the interior points of omega.
Here, I added a piece of code that expands omega to provide space to add
the boundary values similar to what you already did in the main body of the code.
So, in this function you need to first compute the stream function here,
add the boundary values to omega here and compute the derivatives of omega here.
The derivatives of the interior values should be put in
d_omega_dt which should be an n minus 2 by m minus 2 column vector,
and after you write this function, your code should be complete.
Good luck in writing this code.
I think this particular solution method pulls together
a lot of the material you learned in this specialization, and
after you get a working code that passes the greater you can go ahead and
compute the Karman street one of the most iconic images in fluid dynamics.
We'll see how to do that next.
I'm Jeff Chasnov, thanks for watching and I'll see you in the next video.