Introduction to Computational Fluid Dynamics - Numerics - 4 - Classic Solver Algorithms
摘要
TLDRIn this lecture, Professor Steve Miller discusses advanced solver algorithms in computational fluid dynamics (CFD), particularly focusing on incompressible and compressible flow algorithms. Key topics include the artificial compressibility method to address oscillatory pressure patterns, the SIMPLE algorithm for pressure-linked equations, and the PISO method for implicit pressure correction. The lecture also covers preconditioning techniques for low-speed flows and the FDB method for high-speed flows, highlighting their significance in commercial CFD solvers. The importance of eigenvalues in determining information travel speed and the role of parameters in the FDB method are also discussed, along with the checkerboard problem in numerical solutions.
心得
- 📊 Understanding advanced CFD algorithms is crucial for accurate simulations.
- 🔄 The artificial compressibility method helps eliminate pressure oscillations.
- 🛠️ SIMPLE and PISO are key algorithms for solving incompressible and compressible flows.
- ⚙️ Preconditioning techniques improve convergence in low-speed flows.
- 🚀 FDB is effective for high-speed and hypersonic flow simulations.
时间轴
- 00:00:00 - 00:05:00
Professor Steve Miller introduces classic solver algorithms in computational fluid dynamics (CFD), focusing on incompressible and compressible flow algorithms. Previous discussions included steady and unsteady simulations, time integration, and numerical methods.
- 00:05:00 - 00:10:00
The artificial compressibility method (ACM) is explained as a technique to solve incompressible flow fields, addressing oscillatory pressure patterns in numerical solutions. The ACM modifies the continuity equation to include artificial compressibility terms that vanish at steady-state solutions.
- 00:10:00 - 00:15:00
The governing equations for incompressible flows are presented in non-dimensional form, with a focus on the continuity equation modified by artificial compressibility. The introduction of an artificial density term aims to eliminate fictitious oscillations in pressure calculations.
- 00:15:00 - 00:20:00
The eigenvalue analysis of the ACM is discussed, highlighting how disturbances in flow travel at modified velocities due to the introduction of artificial compressibility. This adjustment helps in solving stiff systems of equations more effectively.
- 00:20:00 - 00:25:00
The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm is introduced as a method to eliminate oscillations in pressure and velocity solutions. The staggered grid approach is explained, along with the predictor-corrector procedure used in the algorithm.
- 00:25:00 - 00:30:00
The SIMPLE algorithm's steps are outlined, including initial pressure guesses, momentum equation solutions, and pressure corrections. The iterative nature of the algorithm is emphasized, along with its application in commercial CFD codes.
- 00:30:00 - 00:35:00
The PISO (Pressure Implicit with Splitting of Operators) method is presented as an alternative to SIMPLE, allowing for larger time steps and improved stability without iterative procedures. The governing equations for PISO are discussed, including predictor and corrector steps.
- 00:35:00 - 00:40:00
Preconditioning techniques are introduced for low-speed and high-speed flows, addressing the stiffness of equations in these regimes. The importance of modifying eigenvalues to improve convergence in numerical solutions is highlighted.
- 00:40:00 - 00:45:00
The FD (Finite Difference) method is discussed for compressible flows, emphasizing the need for flow field-dependent variations to handle high-speed and low-speed regions effectively. The role of parameters in controlling numerical stability and accuracy is explained.
- 00:45:00 - 00:54:15
The lecture concludes with a summary of the discussed algorithms, emphasizing their applications in commercial CFD solvers and the importance of understanding the underlying principles for effective simulation of fluid dynamics.
思维导图
视频问答
What is the purpose of the artificial compressibility method?
The artificial compressibility method is used to solve incompressible flow fields and eliminate oscillatory pressure patterns in numerical solutions.
What does SIMPLE stand for in CFD?
SIMPLE stands for Semi-Implicit Method for Pressure-Linked Equations.
What is the main advantage of the PISO algorithm?
The PISO algorithm allows for pressure corrections without iterative procedures, enabling larger time steps and improved stability.
What are preconditioning techniques used for in CFD?
Preconditioning techniques are used to improve convergence in low-speed or creeping flows by modifying the stiffness of the system of equations.
What is the FDB method used for?
The FDB method is used for handling high-speed supersonic and hypersonic flows, particularly in scenarios with shockwave interactions.
How does the SIMPLE algorithm work?
The SIMPLE algorithm involves a predictor-corrector approach to iteratively solve for pressure and velocity corrections in incompressible flows.
What is the significance of eigenvalues in CFD algorithms?
Eigenvalues determine the speed at which information travels through the flow, affecting the stability and convergence of numerical solutions.
What is the role of the parameters s1, s2, s3, and s4 in the FDB method?
These parameters control the influence of convection and diffusion terms in the governing equations, adapting the scheme to the flow characteristics.
What is the checkerboard problem in CFD?
The checkerboard problem refers to oscillatory patterns in pressure solutions that arise from numerical artifacts in incompressible flow simulations.
What is the focus of the next lecture in this series?
The next lecture will focus on stability analysis for simpler CFD schemes and how convergence is measured.
查看更多视频摘要
3 Type of People You Should Never Trust | Denzel Washington Motivation
Belajar Colab Part : 2 Machine Learning Library
Understanding 3D Reconstruction with COLMAP
Microsoft Build 2025 Keynote: Everything Revealed, in 14 Minutes
Long Distance Relationships Tested! |
Pelatihan SMILE Logistik Provinsi Kaltim, Kalsel & Kaltara 19-22 Mei 2025 - Sesi Panel
- 00:00:00welcome back to introduction to
- 00:00:01computational fluid dynamics I'm
- 00:00:03Professor Steve Miller today we'll be
- 00:00:05talking about the classic solver
- 00:00:06algorithms which reside in many
- 00:00:08commercial CFD solvers in the previous
- 00:00:13class we discussed the idea of unsteady
- 00:00:15or steady simulations we then looked at
- 00:00:17time integration and implicit explicit
- 00:00:19schemes and we of course looked at the
- 00:00:21numerical methods that operate to create
- 00:00:24these types of solutions today we'll be
- 00:00:26talking a little bit more advanced
- 00:00:27solvers for particular types of being
- 00:00:30compressible and then compressible flows
- 00:00:31in particular we will be talking about
- 00:00:34the so called incompressible flow
- 00:00:36algorithms ACM simple and piezo and then
- 00:00:40we'll look at the predominantly used for
- 00:00:42compressible flow algorithms piezo which
- 00:00:45is a modification of the incompressible
- 00:00:47piezo algorithm preconditioning
- 00:00:48techniques and of course F dB
- 00:00:51let's first look at the artificial
- 00:00:52compressibility method this is basically
- 00:00:55used to solve incompressible flow fields
- 00:00:57which of course comes from the knee
- 00:01:00which is a little bit misleading is of
- 00:01:02course it's compressible these are
- 00:01:04basically pressure based formulations
- 00:01:06and are used for incompressible flows to
- 00:01:08keep pressure fields from becoming
- 00:01:10oscillatory what that means is that
- 00:01:12you'll see in a particular flow a
- 00:01:14checkerboard like pattern in the
- 00:01:16pressure data from the solver when you
- 00:01:19look at the solution which appears to be
- 00:01:21converged
- 00:01:22so these checkerboard patterns appear
- 00:01:24for certain reasons which you'll see in
- 00:01:26a few minutes the ACM method is used to
- 00:01:29of course try to eliminate these
- 00:01:30checkerboard patterns in the pressure
- 00:01:32solution this formulation is based on
- 00:01:36field variables often involving the
- 00:01:38pressure velocity and temperature where
- 00:01:40compressible flows who are of course
- 00:01:41look at and examine densities Momentum's
- 00:01:44and energies so the main difficulty in
- 00:01:48these incompressible solutions is of
- 00:01:49course the calculation of the pressure
- 00:01:51and that's apparent in simpler methods
- 00:01:54where you find like an oscillatory
- 00:01:56pattern and the pressure which is
- 00:01:57fictitious because of the numerical
- 00:01:59method let's return to our
- 00:02:02incompressible governing equations we've
- 00:02:05written these in index or Einstein
- 00:02:06notation and in non-dimensional form we
- 00:02:09have continuity for steady flows
- 00:02:13and and here the variables with their
- 00:02:18non dimensionalization save the Infinity
- 00:02:20al-saleem scale and of course new some
- 00:02:24infinities at ambient viscosity now in
- 00:02:28the artificial compressibility method
- 00:02:30the continuity equation will be modified
- 00:02:32to include artificial compressibility
- 00:02:34terms these will vanish when
- 00:02:36steady-state solutions are reached for
- 00:02:38example partial Rho tilde partial T
- 00:02:40tilde plus VII is zero so you notice we
- 00:02:43not much line C's equations here we
- 00:02:45write with a tilde as an artificial
- 00:02:48density so this artificial density we
- 00:02:51want to drive to zero for a particular
- 00:02:53or incompressible flow you notice that
- 00:02:57term is missing on the continuity
- 00:02:58operator nonetheless we are trying to
- 00:03:02equate this term to the product of the
- 00:03:05artificial compressibility factor beta
- 00:03:07so we've introduced a new factor beta in
- 00:03:10our equations of motion which we call
- 00:03:12artificial and it's just a coefficient
- 00:03:14which we set and we try and drive this
- 00:03:17term of course it's a small positive
- 00:03:19coefficient but we try and drive these
- 00:03:22terms so that this equation is going to
- 00:03:24be true so let's write out our equations
- 00:03:27in vector notation so it can be more
- 00:03:29compact and we'll try and find an eigen
- 00:03:31value of it in the first equation we
- 00:03:33have of course our vector form or if the
- 00:03:36definition is in W ad and B here w a WB
- 00:03:41w this is all put and of course a vector
- 00:03:44form for pressure and velocity so for
- 00:03:46example in W we had a partial P partial
- 00:03:48T and parcel BJ personalty right there
- 00:03:51the coefficients of a and B are written
- 00:03:54down here along with their particular
- 00:03:56derivatives so you notice that these
- 00:03:57coefficients are generally known in one
- 00:04:01particular space and time and CFD
- 00:04:03solution let's for an example look at
- 00:04:06the eigenvalues of a rain here to find
- 00:04:12eigenvalues we would take of course the
- 00:04:15determinant of AI minus lambda I with
- 00:04:17died any matrix or lambda I are the
- 00:04:20eigenvalues for each particular order of
- 00:04:24the matrix so if a is of course four by
- 00:04:26four we should
- 00:04:27for eigenvalues if we do the hard work
- 00:04:30we can find that the eigenvalues are
- 00:04:31going to be u u and u plus the square
- 00:04:36root of U squared plus beta and u minus
- 00:04:38the square root and u minus the square
- 00:04:41root of U squared plus beta this is very
- 00:04:44much like the eigenvalue analysis which
- 00:04:45i talked about earlier in the class here
- 00:04:48we call beta is now viewed as an
- 00:04:51artificial speed of sound or artificial
- 00:04:53compressibility of course if get beta
- 00:04:55goes to zero then we return and find the
- 00:04:57original form the caressable equations
- 00:04:59the original eigenvalues this can
- 00:05:02interpret it very physically in that you
- 00:05:04can view that the disturbances in the
- 00:05:07flow are traveling at a velocity U
- 00:05:09velocity U and have law say u Plus u and
- 00:05:13a velocity of u minus u so you can see
- 00:05:17if u is near zero or perhaps the mean
- 00:05:20flow velocity for compressible flow then
- 00:05:22I'll have some eigenvalues which are
- 00:05:24very very close to zero this is a
- 00:05:26problem and it might maintain or create
- 00:05:28a system which is rather stiff and that
- 00:05:30is meaning it's slang in a sense that
- 00:05:33it's hard to solve with a linear algebra
- 00:05:35technique or a numerical technique in
- 00:05:39general nonetheless by adding beta you
- 00:05:41can see that we move the potential eigen
- 00:05:43values for some particular problems away
- 00:05:45from their zero components and this is
- 00:05:48excellent and a very wise thing to do is
- 00:05:51of course then the information will
- 00:05:52travel the numeral data through the
- 00:05:54solution added to domain quicker we
- 00:05:57would like to keep beta high enough such
- 00:05:58that the pressure waves will be allowed
- 00:06:00to travel enough to balance the viscous
- 00:06:02effects of the solutions which are
- 00:06:04usually found which we're going to show
- 00:06:06in this particular ACM method through
- 00:06:08the typical crate Nicholson method note
- 00:06:10that we showed in a slide in the
- 00:06:12previous lecture the crank Nicholson
- 00:06:14method now by doing this by having some
- 00:06:18positive coefficient beta will now have
- 00:06:20a well-defined coefficient or set of
- 00:06:22equations so the ACM method can be
- 00:06:25summarized in by simply adding in some
- 00:06:27artificial speed of sound or
- 00:06:28compressibility factor to the
- 00:06:30incompressible equations of motion and
- 00:06:32simply solving them we have a typical
- 00:06:34well-known technique like crank
- 00:06:36Nicholson by doing this we will have
- 00:06:37removed this so-called checkerboard
- 00:06:39pattern we can see in some
- 00:06:41Solutions you might come across this
- 00:06:43yourself one day in some commercial
- 00:06:44solvers and you'll have to change the
- 00:06:46solver method especially if you're an
- 00:06:49incompressible flow now if we do not use
- 00:06:52this particular scheme of artificial
- 00:06:54compressibility then two other most
- 00:06:56popular approaches are the so-called
- 00:06:58simple and piezo approaches neither
- 00:07:00these are really simple to understand or
- 00:07:02implement but their acronyms are aren't
- 00:07:07just kind of a nice one because of
- 00:07:08course it turns out to be simple but
- 00:07:10it's nothing but simple so you might be
- 00:07:13wondering what simple stands for it's an
- 00:07:14acronym and it's called which I'll call
- 00:07:16it simple for now on to just save time
- 00:07:19it'll be the simple implicit method for
- 00:07:21pressure linked equations what a
- 00:07:23mouthful this also can be used to try
- 00:07:27and eliminate the so-called oscillations
- 00:07:29of the solution which we call with like
- 00:07:31a slang word of checkerboard patterns of
- 00:07:33velocity R pressures in each direction
- 00:07:35at every other node in the solution so
- 00:07:37if you plot say a plane your solution
- 00:07:39you see a checkerboard prop pattern is
- 00:07:41it cannot be there it's actually a
- 00:07:43numerical artifact and if you look at
- 00:07:46the actual residual of the equations
- 00:07:48you'll find out that the mass is
- 00:07:49actually not being conserved this is a
- 00:07:51problem too so you're not really finding
- 00:07:53a numerical solution to your equations
- 00:07:55of motion for incompressible flows and
- 00:07:57you'll find that these particular
- 00:07:58difficulties can be bypassed through the
- 00:08:01use of so-called staggered grids within
- 00:08:03the outer room those is simple what that
- 00:08:05means is we might have a standard grid
- 00:08:07which we create with our grid generation
- 00:08:09package and then we're going to be using
- 00:08:12the solver to only look at it from an
- 00:08:14artificial sense as a stator grid for
- 00:08:17example one set of the solver will look
- 00:08:19at points 1 3 5 7 etc and the other part
- 00:08:23of the solver will look at nodes 2 4 6 8
- 00:08:27etcetera so one will look at even a
- 00:08:28normal look at odd and if you go up to
- 00:08:30the next row then we stagger it by an
- 00:08:32increment of one we've already talked
- 00:08:34about the predictor-corrector procedure
- 00:08:36in the previous class here we'll look at
- 00:08:38it again the I remember the idea in a
- 00:08:40particular character procedure is to the
- 00:08:42that we make an initial prediction which
- 00:08:44is like a mid step and then we correct
- 00:08:46it to the next time step or iteration of
- 00:08:49the solver here we're going to do a
- 00:08:51predictor corrector struck with the
- 00:08:53so-called success
- 00:08:54pressure corrector step of course that's
- 00:08:57part of the name of the acronym and
- 00:08:59we'll say that now we're going to
- 00:09:00decompose the pressure which is the
- 00:09:02pressure we really care about as an
- 00:09:04estimated pressure with a bar and a
- 00:09:06prime these are not the so-called
- 00:09:08average pressures or fluctuating
- 00:09:10pressures which we've already talked
- 00:09:11about and which we'll talk about later
- 00:09:12in turbulence modeling in the class
- 00:09:15likewise we can do the same
- 00:09:17decomposition that there's an estimated
- 00:09:19velocities and corrected velocities
- 00:09:22written with the bars and prime
- 00:09:24respectively so we now have a type of
- 00:09:26decomposition
- 00:09:26nonetheless we'll link this in with our
- 00:09:29particular solver into something like a
- 00:09:31predictor correct or type method in the
- 00:09:35simple method the semi implicit method
- 00:09:37for pressure linked equations will try
- 00:09:39and relate the pressure Corrections to
- 00:09:41the velocity Corrections by
- 00:09:42approximating momentum equations so now
- 00:09:45that we have this decomposition we will
- 00:09:46approximate them among equations as you
- 00:09:49see here once again you can see now we
- 00:09:51have a right-hand side of a correction
- 00:09:53and I excuse me a yes correction and the
- 00:09:58left-hand side will be the original
- 00:09:59left-hand side momentum equations both a
- 00:10:01nonlinear term you can also solve this
- 00:10:03for u Prime like here on the right with
- 00:10:05delta T terms now you'll notice
- 00:10:07something interesting there's no viscous
- 00:10:09terms in here that's of course because
- 00:10:11it's not part of this part of the
- 00:10:13approximation we've dropped the
- 00:10:14nonlinear terms I'm excuse me the
- 00:10:16viscous terms no we'll try and combine
- 00:10:20these particular sets of equations with
- 00:10:22the continuity equation so we decompose
- 00:10:24the momentum equation with the corrector
- 00:10:25step and now we're going to insert it
- 00:10:27and combine it with of course the
- 00:10:28continuity equation this is the
- 00:10:30so-called and when we find the rather
- 00:10:32popular pressure correction Poisson
- 00:10:35equation named after of course the
- 00:10:37mathematician but he didn't find this
- 00:10:39equation it just has a similar form to
- 00:10:41his of course famous PDE which also
- 00:10:43carries his name so now we'll have a
- 00:10:46correction of pressure with index I I
- 00:10:49will go as the native density over the
- 00:10:51delta T the time step times the change
- 00:10:54of the axial divergence of velocity
- 00:10:57minus the predicted divergence and
- 00:11:00they'll be equal to this right hand side
- 00:11:02upon simplification for I equals 1 and 2
- 00:11:04this is a 2d example in formulation it
- 00:11:07can easily be extended as
- 00:11:08reading now we still have to enforce
- 00:11:11explicitly mass conservation at the
- 00:11:13current iteration before I go to the
- 00:11:14next time step of course that is to make
- 00:11:17the corrector step so we solve these
- 00:11:20particular equations through the simple
- 00:11:22decomposition of the predictor equation
- 00:11:24method through the simple algorithm now
- 00:11:26let's look at the actual form of the
- 00:11:28simple algorithm I've made a photocopy
- 00:11:32from a single part of Chung's book on
- 00:11:35CFD which shows the method and of course
- 00:11:39he drives it in much greater detail than
- 00:11:40we have in this class but let's go
- 00:11:42through the method together so that we
- 00:11:43can understand it and keep them
- 00:11:45equations in mind of the simple method
- 00:11:47as we look at it step one or Part A will
- 00:11:53guess a pressure P bar to each grid
- 00:11:56point for the initial iteration we can
- 00:11:57simply set it to an ambient value or
- 00:12:00anything we want it's like the initial
- 00:12:01condition then we'll solve the momentum
- 00:12:03equations to find V bar
- 00:12:05that's the predictor step at the
- 00:12:07staggered grade I plus 1/2 I might as
- 00:12:09one have j plus 1/2 and shave my as 1/2
- 00:12:11so here they're not going to every other
- 00:12:13grid point like we were showing or
- 00:12:15discussing earlier we're actually going
- 00:12:17to half grid points or points in between
- 00:12:19this is just like a midpoint or method
- 00:12:21which advances in time but now we're
- 00:12:23doing the same operation to space so you
- 00:12:26see this is the motivation for the
- 00:12:27method it's exactly equivalent then I'll
- 00:12:29solve the pressure correction so we've
- 00:12:31made the prediction now we make the
- 00:12:32pressure correction which we showed of
- 00:12:34course for P Prime in the previous page
- 00:12:36and since the corner grid points will be
- 00:12:38avoided is we can't perform this
- 00:12:40operation on the corners then the scheme
- 00:12:43will be so-called semi implicit and not
- 00:12:45fully implicit as we previously shown so
- 00:12:48it's a semi implicit method the
- 00:12:50correction of the pressure and the
- 00:12:52velocity will then be performed with the
- 00:12:54equations we've shown and here they are
- 00:12:56of course again we'll find the pressure
- 00:12:58at the new step as P U and be based on
- 00:13:01of course the predicted values and of
- 00:13:06course the corrected values so that's
- 00:13:08step D a is shown here as a substitution
- 00:13:12and notice the predictor step of course
- 00:13:14contains the viscous terms that's very
- 00:13:16important
- 00:13:17now now that we've gone to the next time
- 00:13:19level we can simply replace the previous
- 00:13:22intermediate values of pressure and
- 00:13:25velocity with the new corrector values
- 00:13:26of P and B and we return to step B we
- 00:13:29then iterate from b to e over and over
- 00:13:33through every iteration so for example
- 00:13:35if we've done 1,000 iterations step
- 00:13:38iteration one will perform a through e
- 00:13:40and iteration two through a thousand
- 00:13:42perform e through b of course and repeat
- 00:13:44and then every iteration can do a
- 00:13:46convergence check this is a simple
- 00:13:49implementation of the 2d incompressible
- 00:13:51simple algorithm which can of course
- 00:13:53excuse lis be extended to three
- 00:13:54dimensions and there's many papers
- 00:13:56written about this so when you select a
- 00:13:58simple algorithm in the CFD code which
- 00:14:00are usually commercially based to
- 00:14:02solving incompressible flows this is the
- 00:14:05algorithm that's actually solving the
- 00:14:06equations we're looking at note this
- 00:14:08does not contain a turbulence model and
- 00:14:10these equations can become much more
- 00:14:12complicated if of course we use this
- 00:14:13simple algorithm with a particular
- 00:14:15turbulence closure so make sure that you
- 00:14:18understand that
- 00:14:19very important point let's look at some
- 00:14:23examples to illustrate graphically
- 00:14:25what's happening with the simple method
- 00:14:26the correction the pressure will
- 00:14:29actually accelerate convergence if we
- 00:14:31used as a standard integration technique
- 00:14:33then integration that is the number of
- 00:14:36iterations and computer wall time and
- 00:14:38power will takes much longer compared to
- 00:14:40a simple method so we can actually
- 00:14:43modify the method slightly even more to
- 00:14:45find even better values but of course
- 00:14:48introducing some parameter alpha it's an
- 00:14:51empirical factor so you see now the
- 00:14:53pressure will be decomposed into its
- 00:14:56predictor step in this corrector step
- 00:14:58but the corrector step has an extra
- 00:15:00coefficient alpha we just numerically
- 00:15:02experiment and find a value of point 8
- 00:15:04which seems to be optimal we don't have
- 00:15:07of course the same two other terms or
- 00:15:09nabla P and nabla P now really the
- 00:15:15previous discussion is for
- 00:15:17time-dependent flows where we have time
- 00:15:19to moon in terms but this can become
- 00:15:21even more efficient if we use and assume
- 00:15:24that the flow is steady so we might use
- 00:15:27the following finite volume
- 00:15:29discretization
- 00:15:30and you might recall that we just talked
- 00:15:32about how we're looking at the grits
- 00:15:34which are staggered and 1/2 steps and
- 00:15:37you might recall that we are looking at
- 00:15:39the particular values at half steps for
- 00:15:41example I've 1/2 so for the control
- 00:15:45volume view we might have a cell
- 00:15:47centered value where my cursor is and we
- 00:15:50might go into I 1/2 direction so the
- 00:15:52control volume for you and of course P
- 00:15:55are actually staggered relative to the
- 00:15:57central control volume for P this is of
- 00:16:00course part of the reason why they call
- 00:16:02it a semi implicit method for pressure
- 00:16:04linked equations it's not explicit it's
- 00:16:06not implicit but it is semi implicit in
- 00:16:08that it's a mixed scheme for study flows
- 00:16:11the simple out room becomes even simpler
- 00:16:14if you will and we can substitute a
- 00:16:16conservation variable fee and alpha will
- 00:16:19be called the under relaxation parameter
- 00:16:21for subscripts P and n B so these P and
- 00:16:25n B's will denote the node under
- 00:16:27consideration within the two dimensional
- 00:16:29finite volume grid now we can improve
- 00:16:32convergence even more through a simple
- 00:16:35little method we can write a sub e minus
- 00:16:38sum of a n B of UE prime goes as AE of P
- 00:16:43prime pima's PE prime then e will be
- 00:16:46easy use of e prime plus te P prime P
- 00:16:49minus PE prime this last point is simple
- 00:16:52is not that important the big picture
- 00:16:55what is important is the overall
- 00:16:56algorithm you'll see there's certain
- 00:16:59advantages this symbol but you'll also
- 00:17:01note that it's an iterative procedure
- 00:17:03and it's semi implicit what if we want
- 00:17:06to try and improve the method a bit and
- 00:17:07look at a truly implicit method which
- 00:17:11might of course have better stability
- 00:17:12characteristics during our Brussels
- 00:17:14schemes
- 00:17:14so another incompressible solver is
- 00:17:16called piezo which is the pressure
- 00:17:18implicit with splitting of operations
- 00:17:20you'll see and have noted that simple
- 00:17:25requires an iterative procedure so it
- 00:17:28will obtain solutions that is PISA
- 00:17:30without interest procedures and with
- 00:17:33very large time steps without much
- 00:17:34computational effort relatives a simple
- 00:17:37scheme and other schemes this is why
- 00:17:39it's one very popular scheme to select
- 00:17:41in computational fluid dynamic
- 00:17:44the conservation of mass will be
- 00:17:46designed to satisfied with
- 00:17:47predictor-corrector steps so is still of
- 00:17:50course uses the same predictor-corrector
- 00:17:51idea the governing equations will now be
- 00:17:55a momentum equation and pressure
- 00:17:57correction equation let's look at them
- 00:17:59now so we have a left-hand side which is
- 00:18:01discretized which is marches in time
- 00:18:03with a velocity descritization between
- 00:18:06the next time step and the current one
- 00:18:08where we know the data and the right
- 00:18:10hand side will be of course the stress
- 00:18:12shear stress tensors and strain with the
- 00:18:16pressure term say partial P partial X
- 00:18:18and partial P partial Y at steps n plus
- 00:18:211 now note this is truly an implicit
- 00:18:24scheme as we discussed before because of
- 00:18:26course we have time values of n plus 1
- 00:18:27remember n is the current time level and
- 00:18:30n plus 1 is the time level which we're
- 00:18:32trying to find and minus 1 would be a
- 00:18:34time level in the past there'll be a
- 00:18:36pressure correction associated with this
- 00:18:38let's look at this the new pressure
- 00:18:40which is discretized will go as a
- 00:18:43negative 4 over delta T of the
- 00:18:45differences in velocities minus this
- 00:18:48viscous term so here s IJ of course is
- 00:18:52the derivatives of some one convection
- 00:18:54in viscous diffusion terms we've lump
- 00:18:55them all together as a right-hand side
- 00:18:58term and we can calculate that usually
- 00:19:00because of course is the previous and
- 00:19:01next time stuff it's a closed form
- 00:19:03equation for of course didn't known and
- 00:19:06unknown field variables at n n n plus 1
- 00:19:08respectively here we listed the rest of
- 00:19:12the equations and we define s IJ and of
- 00:19:16course the shear stress no there's a
- 00:19:18predictor step and then there's a to
- 00:19:20corrector steps so the simple method
- 00:19:22have one predictor step and two
- 00:19:24corrector steps this one simply has to
- 00:19:25you can take a second to look at the
- 00:19:28form of equations in the current time
- 00:19:29levels the first predictor step uses a
- 00:19:32time level denoted by star or asterisk
- 00:19:36superstar and the second one uses a
- 00:19:38double supe star so through this simple
- 00:19:41predictor and corrector algorithm we can
- 00:19:43implement the Pizzo method for the
- 00:19:46incompressible 2d system of equations
- 00:19:49this is seems like a lot of work and
- 00:19:52there are certain characteristics of
- 00:19:53these schemes which are beneficial for
- 00:19:55example their implicit which takes them
- 00:19:57more power
- 00:19:58from the next time stuff but the time
- 00:19:59steps we can make our humongous and
- 00:20:02maintain stability and we eliminate
- 00:20:04those so-called oscillations of the grid
- 00:20:06so we can greatly increase the stability
- 00:20:08of the scheme by splitting s IJ term
- 00:20:12which we defined here the nonlinear term
- 00:20:15in shear stress term viscous term and
- 00:20:18will split into his Dinoland 9 now terms
- 00:20:20for example we would split s IJ with I
- 00:20:25that is this tensor s IJ with
- 00:20:28derivatives I as partial partial X of
- 00:20:31Rho U feet minus K V partial X as an
- 00:20:35example now we can expand this term
- 00:20:39through simple discritization so here's
- 00:20:42the continuum form one can discretize it
- 00:20:43using of course our standard
- 00:20:45differencing finite differencing
- 00:20:47techniques which we talked about in our
- 00:20:48numeric class and we can apply up
- 00:20:51winding so up whining is going to take
- 00:20:53these parts of the stencil at midpoints
- 00:20:56and bias them towards the upwind
- 00:20:58direction for example if the current
- 00:21:01velocity is positive we would use this
- 00:21:03so-called upwind direction and if it's
- 00:21:05negative one so look at the differences
- 00:21:08between the indices the positive
- 00:21:11velocity you use I and I minus one and
- 00:21:13the negative velocity we use indexes I
- 00:21:14plus 1 and I so this means the stencil
- 00:21:17is biased in the upwind direction that
- 00:21:19uses more grid points or excuse me
- 00:21:22nodes for the differencing approach in
- 00:21:25the direction of the wind
- 00:21:27so using these operations we can find
- 00:21:30the P so with a splitting operation
- 00:21:32which we just showed and we'll arrive at
- 00:21:35these particular equations for Rho UV at
- 00:21:37I plus 1 and I minus 1 with of course an
- 00:21:41up winning operator which is equivalent
- 00:21:43to of course positive and negative
- 00:21:45velocities over its high speed you know
- 00:21:49it might be a very very small negative
- 00:21:51velocity a very high speed positive
- 00:21:54velocity so we can also rewrite s IJ and
- 00:21:57these coefficients alpha beta and gamma
- 00:21:59and put them on the scheme once we do
- 00:22:02all that we can try to analyze the
- 00:22:04system for a structured grid and you'll
- 00:22:06find a system like this a tri diagonal
- 00:22:09dominant of beta gamma and alpha
- 00:22:12beta and a gamma and alpha this is the
- 00:22:16diamond off diagonal turn decomposition
- 00:22:18which we talked about of course on the
- 00:22:20first line of slide 15 now that's an
- 00:22:24incompressible flow and the whole piezo
- 00:22:27algorithm which you can of course
- 00:22:28implement when you run the piezo
- 00:22:31incompressible solver unsafe fluent or
- 00:22:33star or other commercial software codes
- 00:22:35you can imagine it's a lot of work to
- 00:22:37implement these just for the laminar
- 00:22:39equations but once again they can also
- 00:22:41be implemented and applied to flows with
- 00:22:44a turbulence model now what if we want
- 00:22:47to consider a flow which has
- 00:22:48compressibility which is perhaps in a
- 00:22:50general sense seen as a flow with mach
- 00:22:52number greater than 0.3 there's some
- 00:22:55modifications in the algorithm which we
- 00:22:56have to examine we would have to split
- 00:22:59SI je once again to Dino Dino parts as
- 00:23:02we discussed and we can use the piezo
- 00:23:04algorithm apply to a compressible flow
- 00:23:06we would once again have a predictor
- 00:23:08step and to correction steps between the
- 00:23:12correction steps we would solve for
- 00:23:14pea-soup star minus PN which is known PN
- 00:23:17is known this is the current time step
- 00:23:19in the known value and certain to obtain
- 00:23:21the new velocities V stubble soup star
- 00:23:24those values would of course be used to
- 00:23:27find and solve for Rho we can subtract
- 00:23:30these two particular equations and find
- 00:23:33this equation which is labeled from
- 00:23:36Chung's book we would look for the
- 00:23:39triple soup star which of course in this
- 00:23:42case with I I would be to enforce the
- 00:23:45conservation of mass we would then find
- 00:23:48solutions which lead to this form where
- 00:23:50VI
- 00:23:51triple soup star and P Double soup star
- 00:23:53is vo n plus 1 and and P n plus 1 so
- 00:23:56this would complete the splitting
- 00:23:58process where V triple star and piece
- 00:24:00double soup star will imply the exact
- 00:24:03solutions of V of n plus 1 P n plus 1
- 00:24:05which can be shown right here in these
- 00:24:07particular equations for these
- 00:24:09particular procedures you can see the
- 00:24:11full derivation and implementation in a
- 00:24:14paper by Asuma Gausman and Watkins from
- 00:24:171986 this is the so-called PSO
- 00:24:20compressible form
- 00:24:22let's now examine the modification
- 00:24:24explicitly for compressible flows we
- 00:24:27would add an additional corrector stage
- 00:24:29which would be incorporated because of
- 00:24:31course now there's coupling between
- 00:24:32momentum energy and pressure
- 00:24:34there's fluctuations of density and
- 00:24:36temperature with the effect of
- 00:24:37compressibility which we did not have
- 00:24:39before and they must be common form so
- 00:24:41we can discretize a piezo scheme in a
- 00:24:43particular time levels like I've written
- 00:24:45down in the middle of the page you can
- 00:24:47summarize a new predictor-corrector
- 00:24:48steps as now we have a momentum
- 00:24:52predictor step an energy predictor step
- 00:24:55an intermediate momentum predictor step
- 00:24:58and another momentum corrector step a
- 00:25:01second one then we'll have an energy
- 00:25:04corrector and finally a third momentum
- 00:25:07corrector we can then and have to use a
- 00:25:09continuity equation based on of course
- 00:25:11the third corrector steps and the
- 00:25:14original density and the pressure
- 00:25:16equation which we can use to find a new
- 00:25:20pressure corrector step equation these
- 00:25:23are all the equations they're used to
- 00:25:25find a two or three dimensional
- 00:25:28compressible flow of the navier-stokes
- 00:25:31equations with the pleasure implicit
- 00:25:33splitting operation method so take a
- 00:25:36second and just to review these
- 00:25:37equations because of time in the class
- 00:25:39you don't have to review these but if
- 00:25:41you're curious and you're running piezo
- 00:25:43methods you'll have to know that these
- 00:25:44are what are programmed and validated
- 00:25:47with most commercial solvers it's very
- 00:25:49rare today to find research-based
- 00:25:51solvers that use the PSO method we've
- 00:25:54looked at incompressible flows and some
- 00:25:57workhorse schemes implemented in most
- 00:25:59commercial solvers to find the solutions
- 00:26:02and we just looked at those equations
- 00:26:03and gave you a very broad overview of
- 00:26:06the algorithm to solve them now we'll
- 00:26:08look at preconditioning techniques now
- 00:26:10for very low speed flows which are
- 00:26:12incompressible or creeping or very high
- 00:26:16speed flows that are supersonic that
- 00:26:18contain regions that are very very very
- 00:26:20slow and locally incompressible then the
- 00:26:24density based formulations which we
- 00:26:25showed their chart traditional are very
- 00:26:27slow to converge this is where the idea
- 00:26:29of pre conditioning techniques came in
- 00:26:31and you might have seen how we formed
- 00:26:34coefficient matrices in the pre
- 00:26:36previous methods of say for example
- 00:26:38simple we might be able to use
- 00:26:39preconditioning or linear algebra
- 00:26:41techniques on those matrices to change
- 00:26:44them from a very very stiff system of
- 00:26:47equations to one which is less stiff or
- 00:26:50one that can be handled by a numerical
- 00:26:51solver with much better ease why are
- 00:26:54very very slow speed flows difficult and
- 00:26:57highly stiff to solve well that's
- 00:26:59because the acoustic speeds are so much
- 00:27:01higher than the flow velocity and the
- 00:27:03mean velocity might be very very very
- 00:27:05close to zero and this means that some
- 00:27:08information is travelling very fast
- 00:27:09through the fluid domain and someone's
- 00:27:11very trimmed very very very slow and so
- 00:27:14it would take a very long for the
- 00:27:16information the Travelodge domain and
- 00:27:17find say a steady solution for a flow
- 00:27:20which is very very slow speed or very
- 00:27:22high speed with a very slow speed region
- 00:27:24so we'll try and examine the transition
- 00:27:27from spatially varying compressible to
- 00:27:28incompressible regions and vice-versa
- 00:27:31this will use the density based
- 00:27:33formulation and we'll use
- 00:27:35preconditioning matrices on the
- 00:27:38particular time-dependent terms this
- 00:27:41will essentially and mathematically
- 00:27:43improve the convection eigenvalues for
- 00:27:45low Mach number or incompressible flow
- 00:27:47regimes let's take a look at these
- 00:27:49formulations now as before we'll rewrite
- 00:27:52our equations in vector form which
- 00:27:55Akopian so this is couldn't perhaps form
- 00:27:57to a curvilinear or transform coordinate
- 00:28:00system using a finite difference
- 00:28:01approach we have a vector U which is our
- 00:28:03solution vector of P the velocities and
- 00:28:05temperature F and G of course our flux
- 00:28:08vectors here will rewrite a as partial u
- 00:28:13partial Q like we've done in the second
- 00:28:15equation and Q will then be this
- 00:28:18particular matrix so we take the partial
- 00:28:20derivative of U with respect to Q and
- 00:28:22then we marae write the first equation
- 00:28:24as second with a equals partial u
- 00:28:26partial Q here we have a alpha sub P and
- 00:28:29beta sub T these coefficients are
- 00:28:32defined at the lower part of the page so
- 00:28:35it's a lot of work to do this type of
- 00:28:36decomposition and rewriting the equation
- 00:28:39in the form of the second equation on
- 00:28:42the page nonetheless you can see is now
- 00:28:44in a vector form which we've shown for
- 00:28:46you now let's try and look at the eigen
- 00:28:48values and the main
- 00:28:49because the eigenvalues of this matrix
- 00:28:51will dictate the types of speeds that
- 00:28:54the information travels through the flow
- 00:28:55as we're finding a steady solution so
- 00:28:59let's look at the eigenvalues of a
- 00:29:00through the same equation we did before
- 00:29:02we'll take the inverse of a times B sub
- 00:29:05I minus lambda I you can see here's a
- 00:29:08and B matrices in the coefficient
- 00:29:10equations you'll see that incompressible
- 00:29:13limits that is the densities going to
- 00:29:16some constant informally will find that
- 00:29:18the eigenvalues become rather stiff
- 00:29:20which is an algebraic set of equations
- 00:29:23which are ll condition an accoustic
- 00:29:25speeds will become infinite this means
- 00:29:27that as the Mach number the flow becomes
- 00:29:29very small and goes to zero and the
- 00:29:32density becomes a constant and that
- 00:29:33cannot fluctuate that is we're solving
- 00:29:35the incompressible set of equations the
- 00:29:37speed of sound will go to infinity this
- 00:29:39is rather troubling this of course it's
- 00:29:41the unphysical thing but is a property
- 00:29:42of making the incompressible assumptions
- 00:29:45that the acoustic waves travel at
- 00:29:46infinity this is not appropriate for
- 00:29:49finding an incompressible flow solutions
- 00:29:51of course in true and compressible flows
- 00:29:53acoustic waves do not in reality travel
- 00:29:56at infinite speed they travel at roughly
- 00:29:59the square root of DM a times R times T
- 00:30:02the ratio specific heats times the gas
- 00:30:04constant times the static temperature
- 00:30:06and so we want to modify these
- 00:30:09eigenvalues in some way and so that we
- 00:30:12can make the system less stiff so we can
- 00:30:14find it will do this through examination
- 00:30:17the first column of the time to Cobian
- 00:30:19which contains the derivative of the
- 00:30:20density versus pect pressure at some
- 00:30:22constant temperature so if we look at
- 00:30:24partial Rho partial P at T it goes as
- 00:30:26gamma over the speed of sound squared
- 00:30:28this is a well-known equation in
- 00:30:30compressible flow it relates the changes
- 00:30:33of density with pressure at some
- 00:30:35temperature and it'll go as gamma over
- 00:30:37this beautys on squared so this
- 00:30:40derivative will vanish for
- 00:30:41incompressible flows in the speed of
- 00:30:42sound go to infinity this is the core
- 00:30:44problem written down mathematically to
- 00:30:46get around this problem we're going to
- 00:30:48modify density essentially by
- 00:30:50multiplying it by beta
- 00:30:52we'll have Rho a beta so T a partial Rho
- 00:30:55partial P at some temperature and
- 00:30:57expanding out the right-hand side we can
- 00:30:59rewrite it as 1 over V R which is now a
- 00:31:02reference velocity my
- 00:31:04one over C sub P which is the
- 00:31:05coefficient a specific heat at constant
- 00:31:07pressure times density times partial Rho
- 00:31:10partial T at some constant pressure you
- 00:31:12can see the mathematical operations here
- 00:31:14from thermodynamics to write up these
- 00:31:16equations and expansions what might V
- 00:31:19sub R be through this introduction we
- 00:31:21might simply set it for the speed of
- 00:31:23sound for compressible flows or the
- 00:31:26square root of say the square of the
- 00:31:29velocity for incompressible flows for
- 00:31:32her house boundary layers as an example
- 00:31:35so basically we do have to set it to
- 00:31:37some reference velocity that shouldn't
- 00:31:40be too much of a problem a signal
- 00:31:42reference velocity for the flow we can
- 00:31:43try to adjust the eigenvalues of the
- 00:31:45pre-ignition system and their
- 00:31:47adjustments will be done through of
- 00:31:49course the eigenvalue definition with a
- 00:31:51and we'll right now that capital lambda
- 00:31:54goes to dyno of u u YY soup star plus
- 00:31:58the speed of sound
- 00:31:59- u superstar - the speed of sound soup
- 00:32:02star we've now defined two new variables
- 00:32:05Y sweep star YouTube star now capital
- 00:32:08lambda will introduce as the diagonal as
- 00:32:11which are the eigen values of the system
- 00:32:14as u u u the swash should be u YouTube
- 00:32:19star plus the speed of sound soup star a
- 00:32:21soups are and YouTube star minus u soup
- 00:32:25start so the three eigenvalues are you
- 00:32:28you and you and the last two are you
- 00:32:29plus the speed of sound and u minus B so
- 00:32:33YouTube star and the speed of sound soup
- 00:32:35star are now defined in this manner
- 00:32:37these are found from of course the eigen
- 00:32:39value operations with the substitution
- 00:32:42of our modified speed of sound relation
- 00:32:46you'll now find when the velocity is
- 00:32:48greater than speed of sound that the
- 00:32:50eigen values will become u plus the
- 00:32:53speed of sound and u minus V the speed
- 00:32:54of sound and if the reference velocity
- 00:32:57is approximately zero are very very
- 00:32:58small then all the eigen values will be
- 00:33:00of the same order so if a very very
- 00:33:02small creeping flow like here if the
- 00:33:05overall velocity of the flow is say 0.01
- 00:33:08meters per second then all the eigen
- 00:33:10values will be modified to become of the
- 00:33:12same order of U this is good this means
- 00:33:14our system is not very Stephanie
- 00:33:17easily solved with a numerical method by
- 00:33:19simply applying this preconditioning
- 00:33:20technique these preconditioning
- 00:33:23techniques are actually pioneered in
- 00:33:25linear algebra systems if we just have a
- 00:33:27linear system of ax equals B we can
- 00:33:29precondition the matrix a and B with the
- 00:33:32same multiplication and find a
- 00:33:35pre-condition system which modifies
- 00:33:36Augen values and then we can solve the
- 00:33:38new system of equations opposed to
- 00:33:41solving the original set which might be
- 00:33:42impossible this is exactly what we're
- 00:33:44doing with this CFD technique we're
- 00:33:46modifying the eigenvalues through a
- 00:33:49modification of the speed of sound
- 00:33:51relation with a reference value
- 00:33:53introduction we can further improve the
- 00:33:55efficiency in time accurate solutions we
- 00:33:57may try and utilize a dual time stepping
- 00:33:59method which we showed previously of
- 00:34:01course for a midpoint type method for
- 00:34:03simple and P so we'll introduce the
- 00:34:05pseudo time derivative so for finding a
- 00:34:07steady flow we will introduce a pseudo
- 00:34:09time derivative and we'll try and drive
- 00:34:11that to the zero in a linearized
- 00:34:14iteration step form so we can now take
- 00:34:16our vector equation which was shown on
- 00:34:19slide 23 the second equation which we
- 00:34:22discussed its formulation and we'll find
- 00:34:25a right-hand side age which of course
- 00:34:27will let be the compress the
- 00:34:29incompressible form the continuity
- 00:34:31equation so the pseudo time step here
- 00:34:34DITA will approach infinity as the
- 00:34:36pseudo time step vanishes and will
- 00:34:38recover of course the steady original
- 00:34:40governing equations that is this term
- 00:34:42will go away over time and we have the
- 00:34:44same terms from the modified vector
- 00:34:47equation that's really the transition
- 00:34:51from incompressible flow to compressible
- 00:34:53and then looking at incompressible or
- 00:34:55compressible flows of course with a
- 00:34:57preconditioning technique let's turn our
- 00:34:59attention to a completely different
- 00:35:01concept for compressible flows that
- 00:35:03contain both high speed flow and very
- 00:35:05low speed flow flow field dependent
- 00:35:08variation was developed to handle high
- 00:35:10speed supersonic flows and hypersonic
- 00:35:13flows where there's shockwave bound
- 00:35:16chill interactions a shock wave boundary
- 00:35:18layer interaction in a very small region
- 00:35:19of the flow we have compressible
- 00:35:21turbulence we have flow which is
- 00:35:23creeping like in the viscous part the
- 00:35:25turbulent boundary where near the wall
- 00:35:27where the flow is essentially zero the
- 00:35:29velocities and very
- 00:35:30very high-speed flow not far away from
- 00:35:32the law where the Mach number might be
- 00:35:3310 or 20 in this case you might have
- 00:35:36flow at the wall which is about Mach
- 00:35:37number zero and away from the wall might
- 00:35:39be about Mach number 20 and so we can
- 00:35:42return to our generalized form of the
- 00:35:44conservation Knab your Stokes equations
- 00:35:46which we shown it's a vector form in our
- 00:35:49fluid dynamics part of the class well a
- 00:35:51partial u partial T plus partial F
- 00:35:53partial X post partial G partial x
- 00:35:55equals zero so here you'll see something
- 00:35:58very simple we can see of course a
- 00:36:00vector form the equations which we'll
- 00:36:02try and operate on now will expand in an
- 00:36:05explicit form for U of n plus 1 which
- 00:36:08would be the solution of the time sub n
- 00:36:10plus 1 a form of the Taylor series upon
- 00:36:12this particular first set of equations
- 00:36:15by expanding the first term in the
- 00:36:17Taylor series you'll find what we I've
- 00:36:19written on the bottom of the page a 27
- 00:36:21well of U of n plus 1 will go as U of n
- 00:36:24plus delta-t
- 00:36:25of partial U of n plus s a partial T
- 00:36:28plus delta T squared over 2 parcel to U
- 00:36:31of M plus SB partial T over squared over
- 00:36:342 plus higher order terms
- 00:36:35so you see I've done something a little
- 00:36:37bit strange here I've kept the second
- 00:36:39order term the first order term the
- 00:36:41zeroth order term and I've thrown away
- 00:36:42the however terms I truncated them and
- 00:36:45I've done the expansion about N and U of
- 00:36:47M plus 1 with delta T now the new time
- 00:36:51levels on the right hand side that's my
- 00:36:53choice to write them that way I've
- 00:36:55essentially written with something like
- 00:36:57an implicit method is unwrapping in time
- 00:36:59level n plus si and + SB so instead of
- 00:37:03an it it's a well defined sort of like
- 00:37:06midpoint method I'm actually choosing
- 00:37:08two particular midpoints between n time
- 00:37:12level N and n plus 1 and I call these SA
- 00:37:14and SB as a operates the odd time
- 00:37:19derivative and SB operates on the even
- 00:37:21time derivative let's see why we've made
- 00:37:23that choice here now we can split and
- 00:37:27rewrite this first and second order time
- 00:37:31derivative terms the right-hand side the
- 00:37:33first one I've written as partial U of n
- 00:37:35plus si over partial T goes as partial U
- 00:37:38n partial T plus si of partial Delta u n
- 00:37:41plus one
- 00:37:42over partial T and the ranges of si must
- 00:37:46be bounded between 0 & 1 I can rewrite
- 00:37:49the second-order equation here in the
- 00:37:52same methodology also the coefficient s
- 00:37:55B is between 0 & 1 so I redefine these
- 00:37:58and I can take these and insert and
- 00:37:59recover my of course Taylor expansion
- 00:38:01here I've written Delta U of n plus 1
- 00:38:04will go as U of n plus 1 minus u n so
- 00:38:06it's just a shorthand notation between
- 00:38:08the two night-time level solutions
- 00:38:10remember these are vectors for row
- 00:38:12momentum and energy we can now
- 00:38:14substitute these equations in the
- 00:38:16previous one and we will find U of n
- 00:38:19plus 1 goes as unit and plus delta T
- 00:38:21plus the SI term plus delta T squared
- 00:38:23over 2 plus the SB term that's what this
- 00:38:25middle equation is we want to do this
- 00:38:28through jacobians and of course we want
- 00:38:31to solve
- 00:38:31perhaps in a finite difference approach
- 00:38:33to a computational domain from the
- 00:38:36physical space which has uniform spacing
- 00:38:38course so we'll introduce the jacobians
- 00:38:39once we do that we can introduce them in
- 00:38:42terms of the convection diffusion and
- 00:38:44diffusion gradient terms and write the
- 00:38:46first and second derivatives of the
- 00:38:47cancer variables in this form this is
- 00:38:50simply our original equation with the
- 00:38:51right hand side moved the left hand side
- 00:38:53moved to the right and then we'll also
- 00:38:56find a second-order term with our
- 00:38:58definitions above you can combine these
- 00:39:00two particular derivatives and find this
- 00:39:02intermediate equation we'll find partial
- 00:39:04T partial T squared goes as partial
- 00:39:07partial X I of AI plus bi times partial
- 00:39:10FJ partial XJ plus partial GJ partial XJ
- 00:39:13plus this mixed derivative we can then
- 00:39:17substitute these terms into our form of
- 00:39:20U plus M of 1 equals u n plus our higher
- 00:39:24order terms and we'll find this large
- 00:39:27equation it will be for Delta of U and
- 00:39:31plus 1 will go as delta T times all
- 00:39:33these terms which you previously
- 00:39:35substituted but you'll see in this
- 00:39:37equation we have the typical flux terms
- 00:39:39on the right hand side plus an sa times
- 00:39:41the partial derivative of the flux terms
- 00:39:43at n plus 1 plus a second delta T
- 00:39:47squared over 2 of the first order
- 00:39:50derivatives plus the second order time
- 00:39:54derivative terms
- 00:39:55coefficient SB with the flux terms plus
- 00:39:58higher order terms so these
- 00:40:01substitutions and changes the variations
- 00:40:06with Espeon SP are simply resubstitute
- 00:40:08in the previous equation and software
- 00:40:10for the time LaLanne plus one that's a
- 00:40:12lot of work and we've done this in a
- 00:40:14very careful way to find an implicit
- 00:40:18scheme of course for f DB why did we do
- 00:40:22this and I'll try to explain this
- 00:40:24physically now the parameters SA and SB
- 00:40:28will have certain physical analogies and
- 00:40:31we'll find them and set them based on
- 00:40:34the current flow solution so SB and si
- 00:40:37are not constant coefficients they will
- 00:40:39vary depending on the local flow
- 00:40:41properties this is a kind of model if
- 00:40:44you will now let's talk about si first
- 00:40:46which is of course operating on the
- 00:40:49first time derivative if si si shows you
- 00:40:52two temporal changes of convection it
- 00:40:54may be viewed from changes in Mach
- 00:40:56number between adjacent nodal points we
- 00:40:59will set si we'd apply that there's no
- 00:41:01changes in the convection fluctuations
- 00:41:03let's look at an example for the
- 00:41:05variation of si for a shock tube problem
- 00:41:08a shock tube problem of course is where
- 00:41:10we have a long pipe and we have a
- 00:41:11high-pressure gas separated by a
- 00:41:13low-pressure gas a diaphragm separates
- 00:41:16them and an explosion occurs which blows
- 00:41:18apart the metal diaphragm and the high
- 00:41:20pressure region is next to low pressure
- 00:41:22region this creates a shock wave which
- 00:41:24travels to the right in the system from
- 00:41:26high pressure to low pressure
- 00:41:27of course the shock will eventually
- 00:41:30reach the end of the tube and will
- 00:41:31reflect here we show particular
- 00:41:33solutions of the shock tube problem we
- 00:41:37also show relative velocities pressures
- 00:41:40and temperatures you can see in this
- 00:41:43particular case si will vary depending
- 00:41:47on the particular solution of where it
- 00:41:49is in the system
- 00:41:50let's look instead of a shock on a bound
- 00:41:52Euler problem for a boundary layer flow
- 00:41:54si would be one that is dependent on the
- 00:41:57fluctuations of diffusion such as the
- 00:41:59boundary layer flows therefore si would
- 00:42:01be dependent on something like the
- 00:42:03changes of the reynolds number of the
- 00:42:04Peck
- 00:42:05a great working model for sa could be
- 00:42:07formed and so fdv might actually have
- 00:42:10models that vary sa depending on of
- 00:42:14course the Flex vectors F and G it would
- 00:42:17be very beneficial to find general
- 00:42:18models for SA and SB we've only showed
- 00:42:21you two particular models si for a shock
- 00:42:23to problem and SB for of course a
- 00:42:26boundary layer problem but this is
- 00:42:28rather cumbersome in that we have to
- 00:42:29create entire models just to vary SA and
- 00:42:32SB through the flow field we'll try and
- 00:42:35seek out just general models and how
- 00:42:37might we do that well to provide these
- 00:42:40variation of these parameters and
- 00:42:42changes with convection diffusion we
- 00:42:44would need to know something about the
- 00:42:46current flow field so we would need to
- 00:42:48make our models dependent on current
- 00:42:50flow field which is kind of like this
- 00:42:51equation here at the top of the page
- 00:42:53we note that SA and SB to find accurate
- 00:42:56solutions for different kinds of flow
- 00:42:58fields have certain physical properties
- 00:42:59you might know that si with Delta F will
- 00:43:03be called s 1 si with Delta G will be s
- 00:43:063 SB with Delta G will be s 2 and s B
- 00:43:09with Delta G will be s 4 so now we
- 00:43:12divided SA and SB into 4 different
- 00:43:15parameters depending on if it's
- 00:43:17operating on the fluxes f or g so s 1
- 00:43:20will look at his first-order convection
- 00:43:22fdv parameter s 3 will be the
- 00:43:24first-order diffusion ftv parameter and
- 00:43:262 & 4 will be the second-order terms for
- 00:43:29the convection diffusion parameters
- 00:43:30let's illustrate these particular
- 00:43:33breakdown of s1 through s4 for these
- 00:43:35parameters the convection parameters
- 00:43:38will call s 1 and s 2 and parameters
- 00:43:40that are altered by diffusion physics
- 00:43:42and mechanisms will call s 3 S 4 s 1 and
- 00:43:45s 2 are going to be more so dependent on
- 00:43:47Mach number and s 3 and s 4 are going to
- 00:43:49be more so dependent on Reynolds number
- 00:43:50and so will look at the first order
- 00:43:52parameters first we can define simple
- 00:43:56models based on particular Mach numbers
- 00:43:58for s 1 and s 2 here's one functional
- 00:44:01form for s 1 and 1 functional form for s
- 00:44:03to the first order F TV parameters can
- 00:44:07be shown on the right for s 1 and s 2
- 00:44:09with Y axis for F TV parameters s 3 and
- 00:44:12s for diffusion based so for diffusion
- 00:44:15based parameters these simple models are
- 00:44:18chosen
- 00:44:18based on say minimum zero or one
- 00:44:20depending on the Reynolds numbers so you
- 00:44:22can see it s1 has two are depending on
- 00:44:25Mach numbers in a very very similar form
- 00:44:27is given for say s3 s4 which are
- 00:44:30dependent on rounds number and the
- 00:44:32peclet number let's look at one
- 00:44:34particular variation for a shock in a
- 00:44:37corner here the flow moves from left to
- 00:44:40right an initial in LOC Inlet shock
- 00:44:43forms and then in the corner another
- 00:44:45oblique shock form so we have two
- 00:44:47oblique shocks which are being turned of
- 00:44:49course by the flow so the flows pretty
- 00:44:52much uniform in before the shock after
- 00:44:54the shock and after the corner shock and
- 00:44:57the corner shocks and the inland shock
- 00:44:59have extremely large gradients so we
- 00:45:02would expect s1 and s2 to be zero in the
- 00:45:05regions between the shocks at the shocks
- 00:45:09themselves who would probably have very
- 00:45:10very high values that's one most two for
- 00:45:13the diffusion parameters we would expect
- 00:45:15s3 and s4 to be zero and we might have a
- 00:45:19rotational flow near the shock where we
- 00:45:21would of course have s3 and s4 as one so
- 00:45:24by looking at this simple problem we can
- 00:45:26estimate what the values are and indeed
- 00:45:28s1 s2 s3 and s4 are found as a matter of
- 00:45:31the solution as a function of space for
- 00:45:34steady solutions filed through F DV what
- 00:45:38s1 through s4 essentially do are
- 00:45:40parameters which scale and weight
- 00:45:42certain terms in the equations of motion
- 00:45:44if say the diffusion terms are dominant
- 00:45:47or the convection terms are dominant and
- 00:45:51those go of course as the first and
- 00:45:53second order derivatives in the Taylor
- 00:45:56expansion of U respectively let's
- 00:45:58summarize our full form of the f DB
- 00:46:01equations and so if you have a research
- 00:46:02solver or commercial solver and you
- 00:46:05select f DB for these very high speed
- 00:46:07flows then these an equation you'll be
- 00:46:09solving the change in the vector U which
- 00:46:13is the field variables at time plus one
- 00:46:16plus delta T times s 1 plus s3 minus
- 00:46:19delta T squared of the s 2 s 4
- 00:46:21parameters plus delta T of the original
- 00:46:25fluxes minus delta T squared over two of
- 00:46:28the modified flexions plus hard our
- 00:46:30terms will go to zero
- 00:46:32so we have a single vector equation for
- 00:46:34the continuity momentum energy equations
- 00:46:37with weighting factors s1 through s4
- 00:46:39which are dependent on the flow field
- 00:46:41solution you can always ask them main s1
- 00:46:44through s4 in advance as part of the
- 00:46:45initial condition but it's honestly best
- 00:46:47to let them be calculated by the initial
- 00:46:50condition itself of densities velocities
- 00:46:52pressures and temperatures so you can
- 00:46:53see the physical phenomena is in these
- 00:46:55flows are actually dictated by the F DV
- 00:46:57parameters themselves and our inherent
- 00:46:58models now what we haven't done is
- 00:47:02discretize this set of equations
- 00:47:04necessarily with say a finite difference
- 00:47:06method of finite element method or
- 00:47:08finite volume method so we would have to
- 00:47:10take this vector equation and discretize
- 00:47:12it with one of these three methods once
- 00:47:15we do that we would also have to apply
- 00:47:17to balance model if we really care about
- 00:47:20modeling turbulence and high-speed flows
- 00:47:22which we'll discuss later in the class
- 00:47:23so this is in the end of the story of F
- 00:47:26DV it's only the beginning and you can
- 00:47:28imagine how a researcher would have to
- 00:47:30take this equation discretize it perhaps
- 00:47:33of a finite element method for high
- 00:47:35order discretization and then apply it
- 00:47:38and even modify these equations further
- 00:47:40to have additional equations in the
- 00:47:43vector form with a turbulence model this
- 00:47:45could actually take years of someone's
- 00:47:47life and even additional years to
- 00:47:50program it in a contemporary CFD solver
- 00:47:52so it's impossible to learn and
- 00:47:54understand any of these massive methods
- 00:47:56especially flow field dependent
- 00:47:58variation in a single class you're
- 00:48:00probably feeling completely overwhelmed
- 00:48:01and that's okay
- 00:48:03it's just a tip of the iceberg of the
- 00:48:05field of CFD which of course is perhaps
- 00:48:07seen as a mild why entity methodology
- 00:48:11and a field let's take away some
- 00:48:13physical points of the fdv method the
- 00:48:15first-order ftv parameters s1 and s3
- 00:48:18will control all the high gradient flow
- 00:48:20phenomena such as shocks and turbulence
- 00:48:21and so where we have high gradients of
- 00:48:23shock waves or very intense turbulence
- 00:48:25s1 and s3 Multani these parameters will
- 00:48:28be calculated through changes of local
- 00:48:30Mach numbers and Reynolds numbers within
- 00:48:32each element of the flow field there
- 00:48:35could be viewed as actual local elements
- 00:48:38of the flow the contours of s1 and s3
- 00:48:40will resemble the axial flow field
- 00:48:42features themselves
- 00:48:44frig's
- 00:48:45and the shock problem contour of s1 and
- 00:48:48s3 would very much like go along these
- 00:48:51types of rotational flows the fact in
- 00:48:54that contours of s1 and 3 s 3 will
- 00:48:56resemble say Mach number density
- 00:48:58contours has been done demonstrating
- 00:49:00this previous example the role in s1 s3
- 00:49:03of course is to provide computational
- 00:49:05accuracy the second order STV parameters
- 00:49:08are a little bit different they will
- 00:49:09provide exponentially proportional terms
- 00:49:13to the first order SVD parameters
- 00:49:15however their role is more to adequately
- 00:49:18provide computational stability and
- 00:49:20artificial viscosity they were
- 00:49:22introduced as you recall through the
- 00:49:23Taylor series of expansion on a second
- 00:49:25order terms if s1 is 0 these terms would
- 00:49:29represent convection terms this applies
- 00:49:32that if s1 is approximately 0 lumen the
- 00:49:34convection is small this might be a
- 00:49:37rather stationary type flow with a low
- 00:49:39velocity the computational scheme will
- 00:49:41automatically be altered to take this
- 00:49:43effect in account when the governing
- 00:49:45equations are predominantly parabolic
- 00:49:47that is very low speed or transonic if
- 00:49:51s3 is approximately zero in our
- 00:49:53solutions it'll likely be a hyperbolic
- 00:49:55type system of equations and the flow
- 00:49:57would be very very high-speed the scheme
- 00:50:00will automatically switch to something
- 00:50:01like an Euler equations where the
- 00:50:03equations are generally hyperbolic and
- 00:50:05are a very high speed hypersonic flow
- 00:50:07without many terms of turbulence or
- 00:50:10shock 6 era if we look at our solution
- 00:50:12and s1 and s3 are not 0 and they're
- 00:50:16rather dominant then we'll see that
- 00:50:19these indicates that the simulation is
- 00:50:21rather a mixed hyperbolic parabolic and
- 00:50:23elliptic in nature which of course is
- 00:50:26like an ave Stokes system of equations
- 00:50:28and the convection diffusion terms will
- 00:50:30be balanced so if s1 s2 and s3 s4 are
- 00:50:36not 0 and they're equal will of course
- 00:50:38have a rather good balance between
- 00:50:40convection and diffusion terms these all
- 00:50:44these cases can actually be seen
- 00:50:46everywhere in a type of flow and this
- 00:50:48will always be the case for
- 00:50:50incompressible flows at low speeds you
- 00:50:52would not typically choose F DB for a
- 00:50:54purely incompressible flow that would be
- 00:50:56a waste of computational resources and
- 00:50:58you should look at something like the P
- 00:50:59scheme or simple or some simpler scheme
- 00:51:02which we looked at in a previous class
- 00:51:04unique properties of the FTC scheme will
- 00:51:08be capable to control the presser
- 00:51:10oscillations adequately without
- 00:51:11resorting to separately hyperbolic terms
- 00:51:13like the Poisson equation to pressure
- 00:51:15Corrections
- 00:51:16that's one wonderful thing about the F
- 00:51:18DV scheme is we won't see pressure
- 00:51:19oscillations like we had in the other
- 00:51:22schemes and you don't have to resort to
- 00:51:24some sort of like intermediate grid
- 00:51:26overlaps or for example like I plus 1/2
- 00:51:29or I minus 1/2 as we saw in for example
- 00:51:33the simple scheme now the F DV scheme
- 00:51:35could be applied to incompressible flows
- 00:51:37if we delicately balance s1 and s3 and
- 00:51:41of course it happens automatically and
- 00:51:43say the near wall region of that
- 00:51:45turbulent boundary there where this flow
- 00:51:47is high-speed the flow is completely
- 00:51:49incompressible that is the global Mach
- 00:51:50number is zero but there still flow
- 00:51:52fluctuations then of course s1 will be 1
- 00:51:54then will set s1 equal to one explicitly
- 00:51:57and we'll let the variation the Primmer
- 00:51:58s3 be determined through the flow solver
- 00:52:01at the end of the day all these floats
- 00:52:03terms and interactions between
- 00:52:05convection diffusion will happen
- 00:52:06automatically if we use a general model
- 00:52:08for the for fdv parameters in this class
- 00:52:12we summarized some major algorithms for
- 00:52:15numerical solutions of incompressible
- 00:52:16and compressible flows these are
- 00:52:18obviously a little bit beyond an
- 00:52:20introduction to CFD class and are often
- 00:52:22taught more so in depth in Graduate CFD
- 00:52:24classes nonetheless we have looked at
- 00:52:27some workhorse algorithms so to speak
- 00:52:30for many advanced commercial and
- 00:52:32industrial solvers you will see a CM
- 00:52:35simple piezo and variations like
- 00:52:38compressible piezo and maybe
- 00:52:40FDB in many of these types of solvers
- 00:52:43incompressible flows are certainly
- 00:52:45handled excellently by the ACM simple
- 00:52:48and piezo algorithm if you have a
- 00:52:50compressibility in your flow your solver
- 00:52:52should handle of course these
- 00:52:54fluctuations for example through the
- 00:52:56modification of piezo preconditioning
- 00:52:59techniques for creeping flows which are
- 00:53:00incompressible or compressible and of
- 00:53:03course the very famous FDB technique for
- 00:53:06very high speed flows and hypersonics
- 00:53:08and high speed SuperSonics now there's
- 00:53:10many other algorithms
- 00:53:12in fact we could have a whole class just
- 00:53:13talking about all the different types of
- 00:53:15CFD algorithms I was careful and just
- 00:53:17choose three or four of the most popular
- 00:53:20ones which you'll see as choices in the
- 00:53:21solvers from this class try and take
- 00:53:25away the overall approach and realize
- 00:53:27what's happening the solvers and why
- 00:53:29some of these choices are made for
- 00:53:30example to eliminate the so called
- 00:53:32checkerboard problem pressure and
- 00:53:34compressible flows or to handle
- 00:53:36seamlessly the different types of
- 00:53:38dominant terms are convection and
- 00:53:40diffusion in a high speed turbulent
- 00:53:42boundary layer next time don't try and
- 00:53:44ground ourselves in stability analysis
- 00:53:47for simpler CFD schemes the same
- 00:53:49stability analysis can be applied to
- 00:53:51these complicated schemes but it's
- 00:53:53usually not performed and only examine
- 00:53:55numerically then we'll talked about
- 00:53:56residual and convergence and look at the
- 00:53:58so-called l1 and l2 norms if we run
- 00:54:00these schemes or schemes we previously
- 00:54:02talked about then we will look at how
- 00:54:04convergence is achieved and how we
- 00:54:06measure air well then formally define
- 00:54:08and look at convergence and give
- 00:54:10examples thank you very much for your
- 00:54:12time today I'm professor Steve Miller
- CFD
- computational fluid dynamics
- solver algorithms
- incompressible flow
- compressible flow
- artificial compressibility
- SIMPLE
- PISO
- preconditioning
- FDB