Talk Abstract:

Seminar on Industrial Problems

ISDP: Inverse Shape Design Procedure

May 28, 1999

Presented by:

Gary S. Strumolo
Scientific Research Laboratory
Ford Motor Company

The talk will be held at 570 Vincent Hall
10:10 am RealAudio (80 kbps)    RealAudio (28.8 kbps)

Most computational fluid dynamics (CFD) analyses today follow the direct method: specify the shape of an object and CFD will compute the velocity or pressure along its surface. This report will show how to perform the inverse operation; namely, how to specify a velocity or pressure distribution and compute the shape of a body needed to achieve this distribution. The procedure, called ISDP, is fairly robust and can obtain the desired shape even when the initial guess is far from the final one. It works by migrating a set of vortex elements, whose strengths are related to the desired distributions, until they reach spatial positions consistent with them. Often, the final shape is reached with only a few dozen iterations, each of which is computed directly without any matrix inversion. Two types of boundaries are allowed: free and fixed. Free boundaries are ones where the velocity/pressure is known but the shape isn't. Fixed boundaries are ones where the shape is prescribed but the velocity/pressure is unknown.

The method is currently limited to potential flow in two dimensions, although the procedure can be extended to three. Sample problems will illustrate its capability, accuracy, and ease of use.

Introduction

Most computational fluid dynamics (CFD) analyses today follow the direct method: specify the shape of an object and CFD will compute the velocity or pressure along its surface. This report will show how to perform the inverse operation; namely, how to specify a velocity or pressure distribution and compute the shape of a body needed to achieve this distribution. The procedure, called ISDP, is fairly robust and can obtain the desired shape even when the initial guess is far from the final one. It works by migrating a set of vortex elements, whose strengths are related to the desired distributions, until they reach spatial positions consistent with them. Often the final shape is reached with only a few dozen iterations, each of which is computed directly without any matrix inversion. Two types of boundaries are allowed: free and fixed. Free boundaries are ones where the velocity/pressure is known but the shape isn't. Fixed boundaries are ones where the shape is prescribed but the velocity/pressure is unknown.

Broadly speaking, there are two types of inverse methods:

• Iterative use of a direct method
• True inverse methods

The first deals with successive guesses at the final shape. At each iteration, the latest shape is analyzed using the direct method and its actual surface velocity or pressure distribution is compared to the desired one. The error is used to construct a correction to the geometry for the next iterate. This approach is one step above trial-and-error, however, with the "corrections" hopefully providing some rationale for making shape changes. In contrast to this, the ISDP is a true inverse method, where the final shape is determined directly through a series of iterations.

We'll begin by describing the procedure for a purely free boundary case. Three shapes of increasing complexity will be sought: a circular cylinder, a spoiler-type wing, and the centerline profile of the Taurus. We will then describe how this technique can be modified to handle mixed boundary (free and fixed) conditions. This procedure has been awarded a US patent. We conclude with a discussion on stability considerations and the extension of the procedure to three dimensions.

ISDP: free surfaces

Let's consider the case of free boundaries alone. The problem can be stated as follows: determine the shape of an object such that the flow over it results in a prescribed distribution of pressure or velocity. Since the flow is assumed to be potential-like, thus pressure (here represented by the dimensionless pressure coefficient C_p) and velocity are related through Bernoulli's equation:

.

Thus, C_p equals one at stagnation points (V = 0) and zero when the local velocity equals that of the free stream. Higher local velocities yield negative C_p values. Because of this relationship, we can assume that we always begin an analysis with a desired velocity distribution.

The procedure employed is best illustrated using the following diagram:

Figure 1. Steps in inverse shape design for the free boundary case [1].

The surface velocity v_s is chosen as a function of distance along the proposed body, say, from the leading to the trailing edge. We begin by making an initial guess for the shape and dividing it into a series of M element chains of length Ds_m (steps a,b). We further subdivide the surface into an "upper" and "lower" half. At the center of each element, we permanently assign a surface vorticity g (s_m) equal to the surface velocity v_s there. Velocity components u_m, v_m at each element m can then be calculated which include the influence of all other elements plus the external free stream flow (step c). Thus, for M elements with vorticity at the centers we have

where is the external free stream flow and Dq_m is the self-induced convection velocity of element m, given in terms of the profile slope b_m by

The profile slope is measured clockwise from the horizontal, as shown in the following diagram:

Figure 2. Illustration of profile slopes [1].

With these velocity components in hand, we can realign the upper and lower surface of our object by treating them as flexible chains with rigid, straight-line links. Each link is aligned to be parallel to the velocity direction at that point (step d). This realignment proceeds along both the upper and lower surfaces from the leftmost point of the profile, which is held fixed throughout the iterations. Numerical round-offs in the calculations may cause the trailing edge computed by marching along the upper surface to be different from that obtained by the lower surface. Since the shape must be closed for a converged solution to exist, a two-part correction is applied at this point. First, the upper and lower surfaces are rotated so that the lines connecting the fixed point to the (possibly) two trailing edge points coincide. Then the surfaces are stretched/shrunk so that the distance from leading to trailing edge is maintained (steps d,e). This stretching/shrinking of the profiles means that the initial guess need not be a small perturbation from the final shape. Indeed, we will show that a solution can be obtained even with the two far apart.

To assist in the rapid convergence of the algorithm, a relaxation factor f is used. This quantity, between zero and one, regulates how much of the newly-computed shape is to be used after an iteration, i.e.,

.

If f is one, then the newly-computed shape is used in its entirety as the new guess. If f is less than one, the new shape is a combination of the previous shape and the newly-computed iterate. At the end of an iteration the velocity is then computed to be

For the iterations to stop this should agree with the desired velocity distribution to within a prescribed tolerance. Ideally, the final shape should have it's velocity/pressure distribution computed directly and compared to the desired one as a last check.

We will illustrate this procedure with a few examples. The first two have analytical solutions to provide the velocity/pressure distributions while the last is computed numerically. This saves us the need to back check our results via a direct solver.

Circular Cylinder (unit radius)

While this may seem to be a simple problem, it is actually quite a challenge because of the shape near the stagnation points. It is well known that the velocity along the surface of the cylinder under potential flow is given by:

where q is the polar angle measured from the front stagnation point. The initial shape chosen was a thin ellipse with a major-to-minor axis ratio of five. Both the upper and lower surfaces were divided into 20 segments by simply marking the ellipse every nine degrees. At the center of these segments we assigned a velocity according to the above formula. The relaxation factor f was set to 0.5. Figure 3 shows the results of our computation; keep in mind that the desired shape is a unit radius circle. The blue line shows the predicted shape after 30 iterations (square markers). The agreement with a circle is extremely good. Given the initial guess, this example shows that the initial and final shapes can be quite different; the algorithm will still seek out the correct shape. It also illustrates an expansion solution, where the initial guess expands outward to reach the final shape.

In the previous calculation we used 20 segments to define each of the upper and lower surfaces. A reasonable question to ask at this point is how sensitive is the result to the "mesh" used. To answer this we re-ran the calculation using only eight segments per upper/lower surface (a reduction of 2.5 times). The number of iterations and relaxation factor was kept the same as before. The results are combined in figure 3.

Figure 3. 30^th iterate using 20 divisions (solid line, squares) and 8 divisions (dashed line, triangles).

For both meshes, the markers still lie on the correct final shape! This clearly illustrates the robustness of the algorithm.

Spoiler-type Wing

The next level of complexity is provided by wing designs. Aerofoils of this type can be used for rear-deck spoilers. We chose a special set of wings for which the shape and velocity/pressure distributions are known analytically. To provide the mathematical basis for our analysis, we will begin by illustrating the Joukowski transformation, which converts a circle in the z plane to an aerofoil in the z plane:

Figure 4. Joukowski's transformation.

The equation for the aerofoil in the z plane is given by

where the polar coordinates (r,q ) of the circle in the z plane can be expressed as

The angles are defined by

If curve C is arranged to pass through the singular point B at (a,0) then the aerofoil will have a cusped trailing edge, which is not representative of practical aerofoils. The Joukowski transformation can, however, be used to produce aerofoils with a rounded trailing edge. To achieve this we introduce an offset e3 and set

By picking values for r₀, e1, e2, and e3 (hereafter referred to as set A) we can create different wing shapes. Putting this all together leads to the following expression for the surface velocity on the aerofoil:

where the bound vortex strength G is given by

The wing shape derived by picking set A = (0.25, 0.02, 0.03, 0.01) is:

Figure 5. Sample wing used in our analysis (inverted to resemble a rear end spoiler orientation).

We can use the same wing profile but change the free stream angle of attack slightly to give quite different surface velocity (or pressure) distributions over the same object shape.

Let's begin with the 5° case. As with the circular cylinder, we start with an initial guess. However, this time our guess will be larger than the final shape, illustrating the algorithm's ability to also construct a contraction solution. The ellipse chosen has a major axis equal to the final chord length (which is assumed given) and a major-to-minor axis ratio of two. Ten segments make up the upper and lower surfaces with velocity values assigned as before.

Figure 6, shows the state of the computation after the 20th iteration. The desired curve is given by the blue solid line. The predicted shape is given by the red square markers. Except for an overshoot at the trailing edge (which is difficult to remove completely) the agreement with the desired distributions is very good:

Figure 6. Comparison of the computed shape (squares) to the desired one (line).

As mentioned earlier, by changing the angle of attack for the free-stream flow, we can obtain a different velocity/pressure distribution for the same object shape. This provides us with another test of the algorithm robustness: apply two distinctly different velocity distributions and see if it can generate the same profile shape. Using the analytical pressure curves , we ran both cases and compared the two computed profiles. These are shown in figure 7 along with the exact solution; the match is quite good.

Figure 7. Squares - 5° case; triangles - 0° case; solid line - exact solution.

Taurus Profile

As a final test of the algorithm, we developed a vehicle shape based on the centerline profile of the Taurus and computed the C_p distribution along it's upper and lower surfaces. Keep in mind that the distribution used is for a two-dimensional body and should not be confused with the centerline C_p from an actual vehicle. This latter distribution, though two dimensional, is taken from a fundamentally three-dimensional flow and the difference can be considerable. For example, consider the C_p distributions for centerline cuts through a circular cylinder and a sphere. The cross sections are the same, namely, a circle. However, the pressure coefficients along the cross section are not:

For our vehicle test case, the upper and lower surfaces were divided into 36 segments each and a circle was chosen as the initial guess. After 20 iterations, and using a relaxation factor of 0.9, we obtained the shape shown in figure 8. While there is some inaccuracy in capturing the C_p distribution near the front stagnation point, the agreement for both the upper and lower surfaces is quite good. The predicted shape also shows remarkable agreement with the desired one, particularly given the coarse mesh used for this type of body.

Figure 8. Plot showing Taurus shape recovered after 100 iterations using a circle as the initial guess.

In essence, the Taurus was "carved" out of the circular mold using the prescribed pressure distributions.

ISDP: Mixed Surfaces

The previous algorithm assumes that all surface shapes are to be computed, under the constraint of having a certain velocity/pressure distribution. But consider the design problem for a rear-end spoiler. We might want to know what the shape of the spoiler should be when it lies next to a rear deck whose shape and position is fixed. Thus, we would like to be able to solve a mixed problem where the boundaries can be divided into two classes. The first contains free surfaces, where the velocity/pressure distribution is given but the shape is to be determined. The second contains fixed surfaces, where the shape is known but the velocity/pressure distribution is unknown. Consider the following problem:

Figure 9. Schematic of rear-end spoiler problem showing fixed, free, and base surfaces.

The fixed surface corresponds to that piece of the boundary that we know and are interested in. The base surface is added to the fixed surface so we have a closed surface over which potential flow can be computed.

How do we proceed with these mixed boundary conditions? The algorithm is as follows:

1. Discretize the initial-guess free surface G_free and the fixed surface G_fixed into M and N segments, respectively. Construct the base surface G_base and append it to G_fixed so that we have two closed surfaces.
2. Assign vortex strengths permanently to G_free.
3. Solve the matrix problem on using a direct solver to obtain which yields potential flow over this composite surface.
4. Compute the velocities on G_free by using the previous procedure in the "all free surface" case but this time including the effects of the additional vortices .
5. Move the segments on G_free according to the velocities obtained from (4). A few iterations at this point is recommended (i.e., looping over steps (4) and (5)) although it may not be necessary to iterate until convergence before proceeding to the next step.
6. Loop back to (3) and repeat until convergence is achieved.

The net result of this procedure is the determination of the shape of the free surface with a prescribed velocity/pressure distribution along it as well as the computation of this distribution along the fixed portion of the boundary!

This new procedure enables us to consider more relevant problems such as the determination of a rear spoiler shape in the vicinity of a prescribed rear deck shape (2-D problem) and a side mirror shroud shape in the neighborhood of the A-Pillar, side glass (3-D problem).

Future Applications and Improvements

Uniqueness of the solutions

The greatest danger in using this approach is prescribing a velocity/pressure distribution for which a solution does not exist, or one where the surfaces cross each other in an unphysical fashion. Assuming a reasonable solution exists, however, a general proof of uniqueness cannot be given, particularly when dealing with a numerical solution. But with additional constraints on the shape, it does appear that for bodies whose upper and lower surfaces are single valued and nonintersecting that a unique solution is obtained. The trick lies in how the slopes of the leading and trailing edges are defined.

For example, in the circular cylinder case we imposed that the two segments emanating from the leading edge did so at prescribed angles (nearly ± 90°). This is not an unreasonable or burdensome constraint since at a stagnation point the velocity is zero so the surface must be normal to the free stream flow just prior to it. A similar constraint could have been imposed at the trailing edge as well but we allowed it to remain free. Despite this, we can see from figures 3 and 4 that the correct profile was obtained. In both the spoiler-type wing and Taurus cases we imposed slope conditions at both the leading and trailing edges.

Stability considerations

Despite the apparent robustness of the algorithm, it is possible for the iterations to sometimes stall or even diverge if the relaxation factor used isn't optimal. One improvement might be to use optimization theory to determine at each iteration a value for the factor that produces the greatest reduction in our residual error estimate. Thus, we could have an iteration sequence where the relaxation factor changes throughout, producing the most rapid convergence to the final shape.

Extension to three dimensions

This is the greatest challenge, but one which clearly must be met if we are to extend the scope of the present work. The direct solutions in three dimensions can be easily done using boundary integral techniques. The question is how do we define the vortex distribution for a three-dimensional body? Consider the following diagram for a section of a body using generalized coordinates:

Figure 10. Surface vorticity components of a vortex sheet (from [1]).

Whereas before we used one vortex along a segment of the two-dimensional body, we now need to use a pair of vortices for a surface patch. These two vortices are not independent, however, since we must have

.

Thus, only one is truly independent. We could consider, say, g₁ as 'bound' vorticity giving rise to a 'shed' vorticity g₂. These bound vortices might then be moved in a fashion similar to that described earlier to reach the final desired shape. This remains a subject for possible future research.

Reference:

[1] Lewis, R. I., Vortex Element Methods for Fluid Dynamic Analysis of Engineering Systems, Cambridge University Press, 1991.

Back to top of page