Tip:
Highlight text to annotate it
X
We have seen the approach – the computational plane approach, where we solved the equations
not in the physical plane, that is xyz, but in a transformed plane, the computational
plane, where it is expressed in terms of the coordinates are psi, eta and zeta. So, this
transformation, this need for solving the equations in the computational plane also
requires us to transform the equations, which are described in xyz into the corresponding
derivatives appearing in the psi, eta, zeta directions. So, this transformation of equations
is done systematically like the following. We will illustrate it with the for a simple
two-dimensional flow case. We will take case of heat conduction.
Let us say that we are looking at dou square T by dou x square plus dou square T by dou
y square equal to Q x, y. This is the equation that we have, where Q is the source term and
T is the temperature. And, of course, as a thermal conductivity and all that is included
in this. Now, when we are trying to solve this in a
computational plane in two-dimensions in terms of psi and eta, then our physical domain may
be for example, like this. This is x and this is y. This is transformed into the computational
domain. And, here if you say that this is A, B, C and D, this is A prime, B prime, C
prime and D prime. Therefore, there is this C prime, D prime corresponds to this and this
AB corresponds to this; A prime D prime corresponds to this. And, lines, which are uniformly distributed
here may correspond to lines, which are curvilinear like that; 1 2 3 4 5; 1 2 3 4 5. So, this
is the physical plane; and, this is the computational plane. And, we do not want to solve the equations
in the physical plane, because that is not a constant y line. But, since this is mapped
on to the computational plane like this in which this is a constant eta line and this
is a constant psi line here, therefore, we can think of doing it like this.
Now, we have write dou T by dou x. We can write in terms of this as dou T by dou psi
into dou psi by dou x plus dou T by dou eta into dou eta by dou x, because x is a function
of psi and eta in two dimensions. And, this is dou T by dou psi into psi x plus dou T
by dou eta eta x. And, what we want is dou square T by dou x square; that is, dou by
dou x of dou T by dou x. And, dou T by dou x is given by this. So, this is dou by dou
x of psi x dou T by dou psi plus eta x.
Now, we can write this as dou T by dou psi into dou by dou x of psi x plus psi x times
dou by dou x of dou T by dou psi plus dou T by dou eta times dou by dou x of eta x plus
eta x times dou by dou x of dou T by dou eta. And, each of these derivatives here, dou by
dou x is like this and that will have two components; this will have two components
and like that. So, we have to evaluate this as dou T by dou psi times psi xx plus psi
x times dou by dou x of dou T by dou psi. So, this will be… Here we take this. So,
dou by dou x of something is equal to… So, we can write this as psi x dou by dou psi
plus eta x plus dou by dou eta of this T. So, this is the operator, which is acting
on T. And, now, this is the same operator which will be acting on dou T by dou psi.
So, we can write this as psi x dou by dou psi of dou T by dou psi plus eta x plus dou
by dou eta of dou T by dou psi. And, here we have plus eta xx – second derivative
dou T by dou eta plus here again we have plus eta x times dou by dou x of this operator
psi x dou by dou psi plus eta x dou by dou eta of dou T by dou eta. So, we can see that
this becomes dou square T by dou psi square.
So, we can write this as dou square T by dou x square as psi xx dou T by dou psi and here we have this eta xx
dou T by dou eta plus psi x whole square dou square T by dou psi square. And, similarly, we have eta x whole square dou square T by
dou eta square.
Then, we have psi x eta x dou square t by dou eta dou psi. And then, we also have psi
x dou x dou square T…
So, psi x eta x dou square T by dou psi dou eta here.
So, we will have plus 2 psi x eta x dou square T by dou psi dou eta. So, this is what we
have for dou square T by dou x square. And similarly, we will have dou square T by dou
y square.
We will have dou psi by dou y here and dou eta by dou y. So, the operator here will become
psi y dou by dou psi plus eta y dou by dou eta of this.
So, we can write dou T by dou y equal to psi y dou by dou psi plus eta y dou by dou eta
of T. And then, we can again differentiate this and finally, we can get dou square T
by dou y square equal to something similar – psi yy dou T by dou psi plus eta yy dou
T by dou eta plus psi y whole square dou square T by dou psi square plus eta y square dou
square T by dou eta square plus 2 psi y eta y dou square T by dou psi dou eta.
So, this equation Q, which is x and y here will now become Q corresponding to psi and
eta.
So, we can put up this whole thing here as psi xx plus psi yy dou T by dou psi plus eta
xx plus eta yy dou T by dou eta plus psi x square plus psi y square dou square T by dou
psi square plus eta x square plus eta y square plus 2 of psi x eta x plus psi y eta y dou
square T by dou psi dou eta equal to Q of psi, eta.
This is the transformed equation from here.
And, one can see that in the transformed equations, we have the metrics of the transformation;
we have psi x, psi y, eta x, eta y. And, double derivatives of psi and eta – these are all
coming.
And, one can immediately see that this equation is much more complicated. There are many more
terms that are coming here for the general case of transformation, where both x and y
depend on psi and eta. In this particular case, it looks like x direction; there is
no change, but y direction is obviously changed. So, some of the derivatives may cancel out.
But, in the general case, all these things are present.
Now, what is also important is that… Whereas, this equation – we have only normal derivatives.
Here we have cross derivative term, which is coming here. And, we have these terms – the
first derivatives, which are not there in this. But, what is especially important is
that we have the cross derivative term, which comes here.
And, we have seen that if were to use a central differencing scheme for… Let us say that
this point i comma j.
Then, if we consistently use second order approximation, then the second derivative
dou square T by dou eta dou psi will be represented as in terms of this point, this point, this
point and this point. And, it would not have any contribution coming from this. So, this
will be T i plus 1, j plus 1; let us put it as a plus b tau i plus 1, j minus 1 plus c
times T i minus 1, j plus 1 plus d times T i minus 1, j minus 1. This whole thing divided
by some delta eta. This sort of approximation will come, where a, b, c, d are some numerical
coefficients. And, especially, what is important is that this discretization around point i
j does not have any contribution from T i comma j. So, when we put this together in
the form of overall…
When we convert this into A T equal to b type of situation, the contribution to the central
term, that is, i j from this derivative is 0. And, there is contribution term from half
diagonal terms.
So, a direct discretization of this will lead to loss of diagonal dominance. This is an
important factor when we deal with this. And, we have deal with… We will have to therefore
use a method, which does not depend on diagonal dominance for the solution of this, for example,
a direct method like Gaussian elimination or an iterative method like conjugate gradient
method. These are all methods which do not depend on diagonal dominance. But, if you
wanted to use Gauss-Seidel method for the solution of this or even the TDMA method,
they depend to some extent on diagonal dominance for assurance of convergence. So, there is
a problem with this. And therefore, one way of solving this is to take this term on to
the right-hand side and evaluate all these things for the kth iteration in terms of k
minus 1 based on the previous values. And, since all these things – T k minus 1 are
known at all i j from the previous iteration, this will go to the right-hand side. In the
discretization here, we will have only terms, which will be contributing to i comma j.
And, here we have normal derivatives appearing.
And, here we have something like the advection term also coming in this. So, this is the
first derivative and this is the second derivative. So, if we were to look at it as some kind
of generic scalar transport equation with advection and all that thing, there is a first
derivative coming here. So, one has to do a proper stability analysis for a discretized
scheme for this derivative and this derivative put together. So, that is another complexity
that we have to deal with. And finally, what we also see is that there are derivatives
– first derivative or a second derivative of the grid coordinate transformation that
is coming into picture. And, for these things to be faithful to the original equation, we
should have a smooth grid – smoothness of the transformation.
If the grid lines are smooth like this, then the second derivatives will be small; otherwise,
these will tend to dominate the whole solution and you may get away from the visca, the diffusion
type of solution that is expected for this. So, we would like to have a…
This is an artificial thing that is being introduced into the equations because of the
transformation. One can show that because of the transformation, the mathematical nature
of this equation does not change; that is, if we are starting with an originally elliptic
equation like this, then even this transformed equation will also be elliptic. If we are
starting with the parabolic, then this will be parabolic if… So, that criterion is…
That condition is maintained; that assurance is there that we are still dealing with a
parabolic or elliptic equation as we started out with… Therefore, the wellposedness is
not affected.
But, there are transformation dependent terms, which are coming here – first derivative
square and then second derivatives. So, the influence of these things should be suppressed
wherever possible. And, that is possible when we have a smooth grid. So, if we have a discontinuous
grid, so that we have discontinuous things, where the second derivative becomes ill-defined
and all that kind of thing, then we have problems.
And, the most important problem is also with respect to the appearance of the cross derivative.
And, this cross derivative is multiplied by these metrics here. And, the contribution
of these metrics, these terms will be higher if you have a non-orthogonal grid. So, if
the grid lines psi and eta are orthogonal and x and y are orthogonal, for example, these
lines here are orthogonal with each other, then this term will be identically 0 and it
can be shown in that way. But, in the general case, we cannot expect orthogonality.
For example, this is coming like this and this is here. So, this is cleanly not 90 degrees
at this point. And, this point here is not 90 degrees. So, that means that for a general
case of transformation, when you go back to here to here, this transformation, the resulting
grids are not necessarily orthogonal.
If they are not orthogonal, then this term, this contribution here is not 0; we will have
to live with cross derivatives. And, these derivatives can be treated in this deferred
correction mode. So, we call this as deferred postponed correction. Because there is one
term here which is not evaluated at the current situation, but at the previous situation value,
the convergence of this scheme is compromised, is delayed. And, the convergence, the correction,
which is coming after some time is going to delay it more if the value of the magnitude
of the correction is large. So, in cases where this term is large, that is, where we have
highly skewed grid…
If we have a skewed grid, which is not orthogonal, then this term may dominate and that may reduce
the overall rate of convergence of an iterative solution of this in a deferred correction
mode. So, these are some of the things that we will have to take into consideration when
we are trying to solve the transformed equations in the computational plane.
Now, this particular form of the equation, where it is like this, where we have psi xx
and eta xx and psi x like this, these are not very useful, because we are trying to
write approximations for these things.
And typically, our computational grid is very simple like this by design. And, it is very
easy to discretize this into equal spacing here and equal spacing here.
So, if you were to rewrite this not in terms of psi xx, but in terms of x psi psi like
this, then x psi psi, that is, dou square x by dou psi square can be written for example,
as xi plus 1 minus 2 x i plus x i minus 1 by 2 delta psi whole square; whereas, psi
xx, which is dou square psi by dou x square will have to be written as psi j plus 1 minus
2 psi j plus psi j minus 1 by delta x square.
In the computational plane, in the physical plane here, the delta x’s and delta y’s
are changing all the time. But, in the computational plane, it is easy to take this rectangular
and divide this into constant delta x and constant delta eta.
The evaluation of this is much more simple than evaluation of this. And, here you have
uniform grids and here you have non uniform grids. So, the accuracies of the evaluation
of this is also compromised.
So, making use of the relation between metrics expressed in terms of psi x and x psi, which
we have derived earlier, we have shown that psi x psi y psi z; eta x eta y eta z; zeta x zeta y zeta z can be expressed
as x psi x eta x zeta; y eta y psi y zeta; z psi z eta z zeta, which we expressed in
terms of the Jacobian. Now, we can go from here to here or from here to here. And, we
are saying that it is easier to deal with dou x square by dou psi square rather than
dou square psi by dou x square. So, we make use of this relation.
To rewrite this in terms of psi like this, x eta square plus y eta square dou square
T by dou psi square plus x psi square plus y psi square dou square T by dou eta square
minus 2 y eta y psi plus x eta x psi dou square T by dou psi dou eta equal to Q psi comma
eta by. In this form, we see that the metrics of transformation in terms of dou x by dou
psi and dou x by dou eta and dou y by dou psi, these are appearing especially when we
have a smooth grid, so that these first derivative terms will go to 0.
I think here we have made the assumption that we have a smooth grid. And, for a smooth grid,
this is equal to 0 and this is equal to 0, so that we can cancel out these terms and
we have a simpler form of this.
So, the idea of the computational plane approach now is to solve this equation for a given
relation between these lines and these lines. So, we are trying to map the physical domain
like this into this. So, this is…
For example, if you say that this is eta 1, eta 2, eta 3, eta 4 and eta 0… So, this
corresponds to eta 0. On this line here, this is eta 1, is constant; this is eta 2, is constant;
and, eta 3 is constant like this. And, these lines here corresponds to psi 1, psi 2, psi
naught like this.
If you want to discretize this, we need to evaluate dou psi by dou eta, dou y by dou
eta and dou x by dou psi like this. And, that comes from the transformation. And, how do
we get the information?
We can say that we know this boundary; and, one can say that this boundary and this boundary
are mapped. And, we know that this varies from say 0 to 1 and 0 to 1 like this. So,
we know what this line is. That is a horizontal line at eta equal to 1. And, we know that
eta equal to 1 corresponds to this line. And, since this is part of the boundary, we can
potentially evaluate the metrics corresponding to this. But, we do not know these lines.
That is part of the grid generation. So, the objective of the grid generation is to find
out all the interior points in the physical plane given that the external boundaries here
is A B C D, correspond to the computational plane of A prime B prime C prime D prime and
subject to the condition that in the computational plane, we have uniform grid spacing. So, this
point here is topologically the same as this point here.
Topologically means that we have i equal to 0 let us say 1 2 3 4 5 and j equal to 1 2
3 4 5. So, this is i equal to 2, j equal to 2; and, this is also i equal to 2 and j equal
to 2. So, if you were to write this in terms of i and j and similarly, i and j here, this
point is 2 comma 2 and this point is 2 comma 2. So, the objective of the grid generation
is to identify all the interior points in the physical plane given the interior points
in the computational plane and given the boundaries of the physical plane to map with the boundaries
of the computational plane. So, that is… And, we have a number grid generation methods
by which we calculation find all the interior points. And, one can even use the condition
of smoothness in order to derive these interior points, so that we can cross out these terms
here.
And, based on that, we will be able to find all the interior points. So, once we know
the interior points, we can evaluate this x eta. For example, x eta for this point here
is dou x by dou eta. So, it means that this is x of this point minus x of this point divided
by this distance. So, that is, we can write this as at constant i here, this is i equal
to 1, 2, 3, 4, 5. And, this is j equal to 1, 2, 3, 4, 5. So, we are writing x eta at
2, 4. So, this is 2, 4. And, we can write this as x 2 comma 5 minus x 2 comma 3 divided
by 2 delta eta. And, delta eta is what we are getting from computational plane. And,
x 2 5 and x 3 5 are known from the grid generation. So, from that… because we know all the x
y of the interior points. So, once we know the interior points, we can evaluate the derivatives
in this way. And, we put those things here.
So, after the grid generation stage once we have found all the interior points and once
we have numerically evaluated all these derivatives for the general case, then we will have an
equation with coefficients based on the grid generation and all that. So, these coefficients
will be constants. And, we have a second order equation with constant coefficients, which
we can then discretize as per the usual second order accuracy or fourth order accuracy as
we wish. And then, we can convert this into an equation like A T equal to b. And, this
A T equal to b can then be solved for T at psi and eta. And, using those values, we can
then come back to T at this point. So, the approach in the computational plane is to
first identify the mapping between the physical plane and the computational plane. The computational
plane is such that it is rectangular in case of two-dimensions. And, in the case of third
dimension, the zeta dimension will also vary between… It is like a cube going from 0
to 1 and 0 to 1 and 0 to 1 in the three directions.
And, within this, we can have uniform spacing of coordinates of lines of constant psi, constant
eta, constant zeta, so that you have small cubes of delta psi and delta eta and delta
zeta depending on the number of grid points we have in this direction, in this direction.
And, what we then say is that this line A prime B prime maps onto the physical plane
– in the physical plane, the line or the curve AB. And, we have as many number of points
here as number of points here. Then, we try to identify where we want to have these 5
points here and then we fix the points on this side and this side. This is already fixed.
And, based on our consideration of where we want to have the grid points, we can fix these
points here. Similarly, the points on these things are fixed and points on these are fixed,
points on these are fixed. Now, once we identify the boundary points,
then we have to get the interior points such that we have a smooth grid. And then, these
lines here map on to the constant psi and constant eta lines here. Based on that, we
derive all the interior points; that is, we find at every i j, what is psi i eta j and
what is x i y j, so that for a given point i j in the computational plane and the physical
plane, we get psi i eta j in the computational plane and x i y j in the physical plane.
Based on that, we can derive the metrics either like this or in this way. And, once the metrics
are derived, we can substitute them into the transformed governing equations to come up
with governing equations with coefficients, which depend on the metrics of the transformation,
these kind of transformations. So, at that point, we have an equation of this kind and
j, the determinant of the transformation is also evaluated as we have already seen in
terms of the metrics, so that we now have an equation, which has derivatives of the
variable in terms of psi and eta and zeta. And, for each of these derivatives, we make
finite difference approximations and then this equation is converted into a discretized
form. And, for good measure, because we are getting
extra terms in this, we will have to do a stability analysis to see that their discretization
scheme that we have got is reasonable. And, we also have to derive the stability conditions
for this particular scheme with these coefficients that are coming here for a given transformation.
And then, based on that, we finally select the suitable discretization scheme. And then,
using the discretization scheme, we convert the transformed equations into A T equal to
b. And, depending on the structure of this A metrics, we may choose to use either a Gauss-Seidel
method. If for example, the diagonal dominance is preserved, if it is not there, then we
have to choose something else, some other method or we have to treat this term as a
deferred correction. We finally use a particular method and then
we solve this to get T as a function of at i, j. And, T at an i, j is also T at psi i,
eta j; and, it is also equal to T at x i, y j. So, this is the computational plane and
this is the physical plane. So, this is how we can get a solution in this approach, which
calls for a fairly complicated approach for a complicated geometry. So, that means that
we have to first generate the grid knowing only the boundary transformation. And, from
the generated grid points, the internal point and the boundary points, we evaluate the metrics
of the transformation. And, we substitute the values of the metrics of the transformation
corresponding to each grid point and then get a transformed equation in the computational
plane. And then, we choose a discretization and then come up with A T equal to b and then
we solve this and then finally get the solution. This kind of solution can be applied to any
arbitrarily complicated geometry. If it is completely a two-dimensional problem,
then it may be possible to generate a grid, which is orthogonal. For example, if one uses
a conformable map techniques, then one can come up with an orthogonal transformation,
which will mean that the cross derivatives will not appear if the cross derivatives are
not there in the original equation. But, for a general three-dimensional case, it is not
possible to ensure an orthogonal transformation. And so, in such a case, the cross derivative
terms will appear and we have to deal with them. And, it would be the grid generation
and transformation would be very easy if we had an analytical expression for the transformation;
that is, analytical expression for psi in terms of x, y, z or x in terms of psi, eta,
zeta like that. In the general case, it is not possible.
For example, if you are considering a case where even a rectangular domain with three tubes
and you have flow going in the outer things like this and if you have another thing, which
is smaller. So, in the general case like this, may not possible to have an orthogonal thing,
but it is still possible to come up with a corresponding non-orthogonal body-fitted grid. And, this non-orthogonal
body-fitted grid is such that this is also a structured grid. A structured grid in the
sense that the structured grid has the grid points here or at intersections of constant
coordinate lines. For example, this may be a constant eta line, this may be a constant
psi line like this. So, all the grid points are located along intersections of constant
psi and eta lines. That is what is known as structured grid.
And for example, every point here is at an intersection of constant psi line and constant
eta line. This is also constant… So, we are looking at this transformed line here.
So, the corresponding point between this and this is along this constant eta line and this
constant psi line. So, the grid points here are intersections of the family of lines of
constant coordinate directions. And, in such a case, a cell here will have four phases
in two dimensions. And, in the case of three dimensions, it will have six phases. And,
not only that, if you know that this is i comma j, you know that immediate left neighbor,
that is, i minus 1 j and right neighbor – i plus 1 j and also the top one, which is i
j plus 1 and the bottom, which is i j minus 1 and also the front and back. So, in this
particular case in a structured grid, you know where you are with respect to the neighboring
points. And, that kind of structure is inherent in this. And, this structured grid will preserve,
for example, the diagonal structure of the A metrics in that we get for elliptic equations
and so on. So, the resulting A T equal to b based on the transformed equation will still remain diagonal although we have
expected terms coming here.
So, structured grid is useful from the point of view of knowing the neighbors, which is
important when we want to deal with higher order accurate approximations for the derivatives.
If you want to have a third order accurate one sided difference, you need to have four
neighboring points. So, that means that you must have i minus 1, i minus 2, i minus 3.
And, that kind of information is readily available in the case of structured grid. But, in unstructured
grid, that information is not available, because these points are not at the intersection of
these coordinate lines. So, in terms of knowing the neighbors, that is an important aspect.
And, this also has possibility of structured coefficient metrics A. And, we know that if
we have a structured coefficient metrics, then the application of the specialized methods
like Gauss-Seidel method or the TDMA and the strong LPC procedure. These kind of methods
are much easier with a structured coefficient metrics and there is an advantage with this.
And, there are also modeling advantages in the case of a structured grid.
For example, if you are looking at a flow domain like this and the flow is developing,
then the velocity profile may be like this here and then it may be finally going towards
this. And, this dimension may be typically some 5 centimeters and this dimension may
be something like 5 meters. So, you have a very long and thin pipe. In a structured metrics,
it is possible to make it into grids such that the aspect ratio – it is delta x by
delta y can be large without affecting the overall accuracy of the solution too much.
But, if you have an unstructured grid and if you have a large aspect ratio, then the
accuracy of the solution is going to be compromised. So, there are certain advantages for a structured
grid. And, that structured grid is actually carried
forward, is used in this case of non-orthogonal body-fitted grid approach to the solution
of the governing equations in which all the computations are done in the computational
plane and then the information is sent back to the physical plane. But, this approach
requires us to have a smooth grid; and then, the fairly complicated transformations of
the equations happens; and then, addition of new terms of the cross derivatives is a
possibility. So, there are certain things – advantages and disadvantages. But, this
approach is a generic three-dimensional approach, which has proved to be very successful in
dealing with fairly complicated geometries. So, in the next lecture, we will look at an
alternative view point, alternative approach to the solution of equations again on a physically
complicated geometry.