Hy!

I'm an MSc student and I'm currently working with a 3-D Finite Difference Method for general curvilinear coordinate systems. Its application is mainly supersonic flow simulations about rocket-like configurations.In order to generate a 3-D mesh for the rockets, the approach is to generate a 2-D plane and rotate it about the symmetry axis. My problem is with the singularity of the coordinate trasnformation that appears in this upstream centreline with this rotation. Does anybody here happen to know how to deal with boundary conditions and, foremost, metrics treatment in this singularity line?

Thanks in advance.

The problem is trivial when you think of it in terms of finite volumes. No flux can travel across a line. In FV terms this singularity does not exist. The extension to FD (finite difference) is analogous. You simply don't take derivatives across that line. Period!

I discussed this topic at length on page 57 of my dissertation, which also makes the connection between FV terms and FD metrics. The topic is addressed in the paragraph, "One may wonder about the applicability of Equation 4.16 near the poles."

Again, my dissertation (and related papers) can be downloaded from the following location,

www.microcfd.com/papers.htm

Hy, Axel.

Excuse me for taking so long for replying. I'm busy working on my thesis. =)

Firstly, thanks for your help. Secondly, I've read your thesis and I could conclud that, actually, no special care should be taken at such singularity because the metric terms automatically goes to zero. Is that right?

However, I still have one doubt and I'll try to explain it by words. Say you have 3-D Finite Difference M. for a rocket configuration. The 3-D grid is generated by a rotation around the upstream centerline of a 2-D grid. Consider I the longitudinal direction, J, the normal one and K the azimuthal (as in your thesis). Moreover, axis 'x' goes longitudinally from the rocket nose to rear, 'z' is the vertical axis and 'y' goes to right hand side.

Consider now the pitching plane and another plane *normal* to it. Consider now, for an example, the metric term: I_y = d I / d y. If you are on the pitching plane, then: I_y = 0. If you are in the plane that is normal to this pitching plane, then: I_y < 0. This derivative increases its (absolute) value near the centerline because the I lines are along with the 'y' axis. IN the centerline, however, it is ZERO, for this plane. Nevertheless, this centerline is where it should have the largest value. This is an inconsistency that I cannot understand!

In other words, for the pitching plane, I_y is zero, which is correct. For the plane that is *normal* to it, it should have the largest value there! That's crazy!

Any idea about it? Best regards

I do not know what is the pitching plane. To construct a FD grid that has the same property as FV, one may put the first grid point at half-grid size away from the polar axis. In this case, the first grid point acts like a cell-center in FV grid . The half-grid size is the half of grid spacing in your "X" direction.

>Excuse me for taking so long for replying. I'm busy working on my thesis. =)

No problem! Actually I have been very busy myself recently, which is why I have retired a bit from this forum.

>Firstly, thanks for your help. Secondly, I've read your thesis and I could conclude that, actually, no special care should be taken at such singularity because the metric terms automatically goes to zero. Is that right?

My pleasure! - It depends. If you start out with a finite volume formulation and make the transition to finite difference via the formulas in my dissertation, then YES. Otherwise, in a finite difference framework, you should only take one-sided derivatives away from the polar axis; never differentiate across the pole.

I have read your questions, but the description of your coordinate system confuses me. Sorry! Could we stick to Figure 4-1 on page 55 of my dissertation. Let's say your rocket is on the launch pad or flying vertically 'up'. Then according to my figure y is up, x is to the right, and z is towards the viewer. In essence, the y-axis is the polar axis, or axis of rotation, along which you have the singularity. If you rephrase your questions according to this convention, I may be able to answer them.

Thanks once again!

>Otherwise, in a finite difference framework, you should only take one-sided derivatives away from the polar axis; never differentiate across the pole.

That's exactly what I do. No fancy stuff with FVM =)

>Could we stick to Figure 4-1 on page 55 of my dissertation. Let's say your rocket is on the launch pad or flying vertically 'up'. Then according to my figure y is up, x is to the right, and z is towards the viewer. In essence, the y-axis is the polar axis, or axis of rotation, along which you have the singularity. If you rephrase your questions according to this convention, I may be able to answer them.

Let me try again, based on this coordinate system you described. Coordinate system (x,y,z) is fixed and (I,J,K) follows the body. Say 'I' is the longitudinal direction, therefore, it follows the 'y' axis 'down'. 'J' spreads normal from 'I' axis and 'K' turns around 'I' axis. The rocket nose is placed in that origin of the figure. The singularity is the K planes collapsing into the polar axis, forward the rocket nose.

Say you want the metric term given by dJ/dz. If you're on the surface z=constant, than dJ/dz = 0. However, if you're on surface x=constant, |dJ/dz| > 0, including the polar axis. So for two differnet surfaces (x=constant, z=constant), two different values are obtained for the same metric term at the same region (polar axis).

That's my doubt. If you could not understand it again (it's difficult it explain it by words), may I send you an image of the problem?

Thanks, Best regards, Biga

>Let me try again, based on this coordinate system you described. Coordinate system (x,y,z) is fixed and (I,J,K) follows the body. Say 'I' is the longitudinal direction, therefore, it follows the 'y' axis 'down'. 'J' spreads normal from 'I' axis and 'K' turns around 'I' axis. The rocket nose is placed in that origin of the figure. The singularity is the K planes collapsing into the polar axis, forward the rocket nose.

Right there is your first problem: You said that (x,y,z) is fixed and (I,J,K) follows the body. So how can 'I' always point along the longitudinal direction? - If you look at Figure 4-1, 'I' is the LOCAL direction away from the body, and 'J' and 'K' are orthogonal to that local direction. Depending on where you are on the body or further out on the mesh, these three axes (I,J,K) can point ANYWHERE!

I think we agree that the rocket, as well as the rest of the 3-D mesh is generated from a 2-D mesh by rotation around the y-axis? Let's also say that the nose of the rocket points in the POSITIVE y-direction.

>Right there is your first problem: You said that (x,y,z) is fixed and (I,J,K) follows the body. So how can 'I' always point along the longitudinal direction?

For a rocket, in my case, 'I' is *almost* along the longitudinal direction! =D

>I think we agree that the rocket, as well as the rest of the 3-D mesh is generated from a 2-D mesh by rotation around the y-axis? Let's also say that the nose of the rocket points in the POSITIVE y-direction.

Agreed!

>For a rocket, in my case, 'I' is *almost* along the longitudinal direction!

This is what I can't visualize, but perhaps you can email me a GIF or JPG of your rocket and surrounding mesh. Are you only talking about the rocket exhaust, or the entire rocket (head to tail)? I am considering the entire rocket and surrounding flow field.

Is your mesh cylindrical on the outer boundary? Because mine (Fig. 4-1) is spherical. Perhaps this is our discrepancy!?!

This is cylindrical! Imagine the following grid plane:

_______outer______

/ |

/ | exit / mesh | z ^ / ___________| | |_________/nose wall .--> x centreline y

This plane is rotated about the centreline. 'I' goes from nose to exit. 'J', from wall to outer boundary. 'K', rotates around the symmetry line (centreline).

Got it?

I 'got it' now! Sorry about the misunderstanding! - Let me think about your case for a while before I get back, but the same concept applies.

>Say you want the metric term given by dJ/dz. If you're on the surface z=constant, than dJ/dz = 0. However, if you're on surface x=constant, |dJ/dz| > 0, including the polar axis. So for two different surfaces (x=constant, z=constant), two different values are obtained for the same metric term at the same region (polar axis).

According to Eq. 4.15 on p.56 of my dissertation, dJ/dz = Jz would be a scalar multiple of the z-component of the outward unit normal to the cell face along the J direction, which is normal to the I direction, and thus 'almost' normal to the y-axis, since your rocket is very slender.

It is difficult to talk about planes of constant z or x, because the intersection of those planes with the local grid yields different values for Jz depending on where you are. If your cell face is pointed along the z-axis, then Jz =1, if pointed along the x-axis, then Jz = 0 (which somewhat agrees with your statement).

On the pole itself (the y-axis) Jz varies as you sweep the 360 degrees around the y-axis. That may appear confusing physically, but computationally, this is what you want. If you unwrap the mesh, you will see that the points along the polar boundary are not really connected on the computational mesh. It just happens so physically that they coincide.

Picture this: Say you started out with a square piece of paper which is marked every 1cm along all of its four straight edges. This square paper would be the equivalent to your computational domain, your grid being perfectly Cartesian, made up of 1cm x 1cm cells. Now you wrap this square sheet around a globe, and through proper stretching of the paper you actually manage that all marked points along two parallel edges coincide with the north and south pole. You can do the same process in a cylindrical plane, it just requires more stretching.

I think when you start writing your code, it will become clear. Sorry, but I cannot explain it any better.

>On the pole itself (the y-axis) Jz varies as you sweep >the 360 degrees around the y-axis. That may appear >confusing physically, but computationally, this is what >you want. If you unwrap the mesh, you will see that the >points along the polar boundary are not really connected >on the computational mesh. It just happens so physically >that they coincide.

That's precisely my doubt. Physically they coincide, but computationally it's what I usually call a brick! For the face of this brick where the physical domain is just a line, I got different values for some derivatives. What you say is that: compute the derivatives as they should be and that's ok? That's what I do now... yet I still can't completely figure out why...

As long as you don't take any differences ACROSS the pole, you will be O.K. If you worked in FV, this whole thing would not even be an issue, which was my first recommendation. Try to understand those formulas I gave you in my dissertation, and eventually you will 'see' it.

I do FD, but I don't take any difference through the pole, so I think I'm safe! Nice!

Thanks for your help!

My pleasure!