A gridless technique for fluid/structural dynamic coupling on flexible membranes

A gridless method has been developed for the simulation of coupled fluid/structural interactions over arbitrary bodies. This method uses Eulerian-based points arbitrarily distributed over the computational domain with no formal connectivity as typically required for a traditional grid. Comparisons are made with known exact solutions for simple two-dimensional model problems. Methods of improving the accuracy of the current implementation by using higher order approximations have been implemented. Accuracy improvement by using point adaption has been investigated. Plane strain and axisymmetric shells have been added to the code structural code PRONTO2D for future fluid/structural calculations. To date, coupled fluid/structure calculations have not been made.

Portions of this document may be illegible in eledronic image products. Smagm are produced from the best available original document.             Higher-Order Taylor Table 4.4

Introduction
There are numerous important fluid and structural problems which require tracking a moving membrane responding to the forces induced by the motion of a fluid. Examples include the deployment and inflation process of a parachute and air bag; the aeroelastic phenomenon known as aileron buzz; and in the biomedical field, the opening and closing of pressure-driven heart valves and angioplasty balloon devices. Safety issues are also a concern in those areas where a rapid decompression may occur, such as a tire blowout at highway speeds. New techniques in mold processing using inflatable bladders may improve quality control by reducing anisotropic strain rates during manufacturing.

I '
Fluidstructure problems, such as the examples listed above, are difficult to solve using existing numerical techniques. The difficulty in simulating coupled problems persists because the solution algorithm, in addition to solving for the fluid dynamic field variables (e.g., density, pressure, velocity, etc.), must also provide for the motion of both the membrane surface and the surrounding mesh. Transient remeshing of unstructured grids is computationally expensive, and convergence and solution accuracy degrade appreciably when the motion of the grid results in a highly distorted mesh. Typically, such grid related errors occur where large-scale structural deformation and/or large flow gradients occur, and it is precisely in those regions where the greatest accuracy is required.
In order to avoid the problems associated with grid generation and grid distortion, we proposed to develop a new technique which provides for the solution of both the field variables and the structural motion utilizing a gridless technology based on a random point distribution. The strength of the methodology lies in its ability to model any geometry in two and three dimensions using "clouds" of discrete points with no requirement for the points to be connected to form a grid as in conventional computational fluid dynamics (CFD) and structural mechanics algorithms. The governing partial differential equations (PDE's) which determine the equations of motion are solved locally at every point in the domain.
The proposed method is not a finite-difference, finite-volume, or finite-element approach since differences, metrics, areas, or volumes, are not computed. The methodology is extremely versatile in that the governing equations of motion for the fluid dynamics are not limited to one class of problems. That is, the method permits one to solve the relevant equations based on the complexity of the problem at hand. Hence, the difficulty of simulating either a potential flow, Euler flow, or full Navier-Stokes flow, and the ability to model compressible or incompressible flows is dictated only by the set of PDE's solved and is easily amenable to a modular code structure. This modularity appears as a simple consequence of the fact that the equations of motion are independent of the point distribution and the solution algorithm is completely grid a free.
In this report we summarize research conducted during FY94 and FY95 on a gridless by Sandia National Laboratories' Laboratory Directed Research and Development (LDRD) program. There were a number of principal investigators for this project. These were divided into smaller teams, each responsible for one part of the project. In this report, each chapter 0 technique for fluidstructural dynamic coupling on flexible membranes. This project was funded describes one aspect of the project and was authored by those team members responsible for that aspect. If a chapter does not list specific authors, e.g., this chapter, then all members of the research team contributed to its content. Chapter 2 presents an overview of gridless methods along with a brief survey of recent work by others in this area. We also give an example of the gridless method applied to a simple one-dimensional problem and present the concept of coupled linear operators for solving systems of differential equations. In Chapter 3 we give a description of the specific methods currently implemented in our gridless code and present some numerical results for our twodimensional model problem: fluid flow described by Laplace's equation. Chapters 4 and 5 describe the results of efforts to improve the solution accuracy. Chapter 4 gives results of research into higher-order representation of the variables over the computational domain, and Chapter 5 presents the results of our investigations into point adaption. Chapter 6 contains results of research into shell elements, which are needed for coupled fluidstructural calculations. Chapter 7 contains a s~m m a r y of the work and recommendations for the direction of continued research in gridless methods. Because of this modular report structure, each chapter is a semiautonomous unit with its own equations and nomenclature. Symbols are defined in each chapter as they are used. We have, therefore, decided not to include a nomenclature section for the entire report.
From its inception to completion, this project had several personnel turnovers. One of the individuals responsible for the original concept of this project and the LDRD proposal was Anthony Thornton, currently Director of Diversity Leadership and Education Outreach. Due to his promotion in early FY94, he was unable to participate in the actual research. Christine Hailey left Sandia at the beginning of FY95 to take a position at Utah State University. She was able to continue her involvement during the summer of 1995 through a university research contract. Walter Wolfe joined the project in January of 1995.
Because of unexpected difficulties with some of the tasks, LDRD funding cutbacks in FY95, and personnel turnovers, we were not able to complete all of the tasks originally planned. Extension of the method to the Navier-Stokes equations and the actual coupling of the fluid and structural calculations were not completed.

Overview of Gridless Methods
Roy S. Baty 8z Walter P. Wolfe

Introduction and General Gridless Philosophy
In this chapter we present a summary of gridless methods, a brief review of recent work by others in this area, and an outline of our current gridless method applied to a simple, onedimensional problem. We also present an optimization of our calculations using the concept of coupled linear operators. The extension of the current gridless method to two dimensions-is presented in Section 2.5 and Chapter 3.
The gridless method we present here uses Eulerian-based points arbitrarily distributed over the computational domain with no formal connectivity as typically required for a traditional grid. Each point in the domain has associated with it a cloud of neighboring points. These neighbors are chosen by a scheme that considers their position relative to the point of interest. Derivatives in the governing PDE's are numerically approximated by assuming a functional dependence of each variable within the local cloud of points, where the number of points in the cloud is greater than or equal to the number of unknown derivatives. The number of points is usually greater than the number of unknown derivatives in order to improve the accuracy of the derivative estimates. We have assumed here that the form of the functional dependence is a Taylor series. The governing partial differential equations are then solved locally at each point in the computational domain. This gives a solution to the governing equations without any a priori assumptions concerning grid or element topology, and allows for addition, deletion, or movement of points with no change in solution algorithms.

Previous Gridless Investigations
In recent years, a few individuals have begun to investigate gridless techniques. Liska and Orkisz' and Batina2 have reported results of initial investigations into the use of gridless methods for two-dimensional applied mechanics and fluid flow applications, respectively. Liska and Orkisz calculated derivatives by assuming a second-order Taylor series expansion about the central point in a nine point cloud. This formulation gives eight equations for the calculation of five derivatives. To solve this overdetermined set of equations, they used a weighted least 3 squares formulation. They chose a weighting factor of l/Aj , where A j is the distance between the cloud's central point and the j'" neighboring point. The justification for the use of this weighting factor was not given. Boundary conditions were imposed through the use of ghost points.
One of the significant points to come from the work of Liska and Orkisz is the importance of point geometry when forming a local cloud. The most straightforward means of selecting points for a cloud is to pick the N points that are closest to the point of interest. However, this criterion often gives an ill-conditioned or singular set of equations due to an irregular distribution of points. They found that an additional criterion is often needed to insure that the points have a more uniform angular distribution about the point of interest. We confirmed this finding in our research.
Batina investigated the use of gridless methods for calculation of flow over 2-D airfoils, using both Euler and Navier-Stokes equations. He assumed a first-order Taylor series expansion for the cloud of points and used a least squares method to determine the derivatives. Details of the least squares method are not given. He also imposed boundary conditions by using ghost points. In the examples presented in reference 2, the points in the computational domain were determined by using an unstructured-grid generator and picking the centers of the resulting triangles as the points. The cloud of neighboring points was determined by choosing the centers of the triangles that shared sides with triangle of interest. This gives a four-point cloud where the points are reasonably evenly distributed around the point of interest. Any problems arising from irregular point distribution were, therefore, not addressed.

The Gridiess Numerical Method
The gridless numerical method developed here consists of two parts: (1) a scheme that approximates a function and its derivatives and, (2) a scheme that combines the function/ derivative approximation with iterative techniques to satisfy the differential equation and the boundary conditions. The solution proceeds in a interactive manner between these two parts, or steps, until a specified convergence criteria is met. To start the iterative process, an initial distribution for the function is assumed. The gridless scheme separates the approximation of the function and its derivatives from the solution of the differential equation so that arbitrary distributions of field points may be applied to compute the solution. In this section, a gridless method is derived for a 1-D model boundary value problem, illustrating the basic approach for linear differential equations.
To develop the scheme for a linear boundary value problem, consider the simple 1-D model problem defined on the closed interval [0,1]: a*$"(x) + a,$'(x) + a*$(x) = 0 and $ ( O ) = a and $(1) = p where a, p , and ai for i = 0, 1,2 are constants. Now, let { x n } be any finite partition of [0,1] such that xo = 0, xN = 1, and x j < xk for all j < k . Note that { x,} is not assumed to be uniformly spaced. Next, recall that the Taylor series for the function + expanded about any partition point xi in [O,l] is given by: where hi = x ix j . Then truncating the Taylor series after the second derivative and rearranging yields: Applying equation (2.4) at E ( E 2 2 ) points neighboring xi produces the linear system: Here the superscript T indicates the transpose of the vector. The boundary conditions given by equation (2.2) are satisfied by fixing the value of the solution @ at the points 0 and 1. For more general boundary value problems, such as potential flow, the derivative matrix D is reformulated to satisfy the boundary conditions on the boundary of the computational domain. The details of this reformulation are contained in subsection 3.2.3 on page 27.
As an example of how equations (2.6), (2.7) and (2.8) are evaluated, a typical calculation might use four grid points to evaluate the matrix equation (2.5): the point where the differential equation is solved and three neighboring points, generating a 3 x 2 matrix for D. Only determined and overdetermined linear systems are used to approximate the function and its derivatives. The determined case (D is a nonsingular square matrix) produces a unique solution for the derivatives, while the overdetermined case ( D is a rectangular matrix with nullity of zero) yields a unique least squares solution. The underdetermined case where there are more unknowns than points has an infinite number of solutions and was not consider in this study.
The differential equation (2.1) defining the boundary value problem can also be written in matrix form as: Equation (2.10) representing the differential operator is written as a matrix with the same row and column dimensions as the derivative matrix D by adding zeros as needed. This representation for L is used to simplify the analysis for the coupled linear operators presented in Sections 2.4 and 2.5.
To see how the basic gridless scheme works for the model problem, the 'local approximation defined by equations (2.5) and (2.9) must be applied at each grid point in (0,l). This requires defining the computational cells used to approximate the solution of the boundary value problem.  Once the local approximation of the boundary value problem has been determined, the gridless scheme computes the solution using an iterative process. The following five steps are performed on the computational domain: 1. Guess a value for the solution @ ( x i ) at each grid point in (0,l). Fix @ at xo = 0 and x N = 1 so that the boundary conditions are satisfied.
2. On each computational cell, compute the approximate derivatives E,(@) from equation (2.5). Since, in general, equation (2.5) represents an overdetermined system of equations, a least squares scheme must be used to compute 6. 3. Substitute 6 and @ into matrix equation (2.9) and define a local residual as a function of 0: 4. Compute @ such that -E < RL(@) < E for some small value E > 0. Any convergeng root-finding scheme can be used. Steps 2 through 4 are repeated at each point in the computational domain. 5. After steps 2 through 4 are completed for each cell in the domain, a global residual is computed and checked for convergence. The global residual, R,, is defined as - (2.14) where N is the number of field points in the computational domain. If the global convergence criteria has not been met, the solution procedure repeats starting with step 2.
To measure the convergence of this scheme, a global convergence criteria is defined as:

Least Squares Solutions for Coupled Linear Systems
In the preceding section, a general numerical scheme was derived to approximate solutions of differential equations on irregular meshes. For arbitrary linear differential equations, this numerical technique yields coupled systems of linear equations on each local computational cell. This section outlines an algebraic method to find the solution of an arbitrary coupled linear system on an arbitrary computational cell. The solution obtained for the coupled system on each computational cell at each solution iteration will be a least squares solution which will minimize the local residual computed in step 3 above for a fixed value of v . Recall that least squares solutions minimize the error in terms of Euclidean length.
Hence, a solution to equations (2.18) and (2.19) does not exist. However, a least squares solution may be found. The least squares solution is the best approximation to equations (2.18) and (2.19) in the sense that it minimizes the error defined in terms of the Euclidean norm.
The problem defined by the coupled linear system (2.16) and (2.17) has been solved in a general algebraic setting for abstract linear operators by Bryan, Kinyon and T~c k e r .~ Their work is applicable to any coupled linear system. Following the development of reference 3, a new matrix G is introduced: The differential operator L is written as a matrix with the same row and column dimensions as the derivative matrix D by adding zeros as needed. This representation for L is used to simplify the analysis for the coupled linear operators. In this example, the matrix D is just a simple matrix with integer entries which has no immediate physical or numerical interpretation. However, the matrix L can be associated with a linear differential operator and a matrix similar to L will be applied to approximate Laplace's equation in Section 2.5.
To compute the least squares solution of (2.21), first notice that the Moore-Penrose inverses of D and L are given by:

Least Squares Solutions for the Gridless Method
This section applies the least squares solutions for coupled linear systems to an example of the gridless method. The specific gridless scheme developed here is used to simulate steady, incompressible, potential flow in two-dimensions. Potential flow is described by Laplace's equation:  The linear partial differential equation (2.40) requires the numerical approximation of the second derivatives of $ . Therefore, the Taylor series used to approximate the derivatives of D will be truncated after the second order terms:

Overview
In this chapter we present a description of the procedures used in our current two dimensional gridless code and then present numerical results for our 2-D model problem.

Cloud Formation
The only requirement for the distribution of points in a gridless method is that it allows an accurate solution of the governing equations over the entire computational domain. Each point has associated with it a cloud of neighboring points that are used to determine the derivatives needed to solve the governing equations. In the following discussion, we have chosen the word focus to designate this central point and the word neighbors to designate the other members of a cloud.
We have assumed that the points are read into the code in a completely random order, i.e., there is no underlying structure to the point array. Therefore, the entire array of points must be searched for each focus to find the points that form its cloud of nearest neighbors. Although this is a time consuming process, it need only be done once for each point during the initial problem setup and, therefore, is not a large fraction of the total run time. Once found, the neighboring points' indices are stored in each focus' data structure. A point's index number is its position in the original point array. We made no attempt to develop efficient functions for this search.
The search procedure proceeds as follows for each point in the interior of the domain: 1. For each point or focus, the distances to all other points are computed.

2.
from reference 10. 3. The closest point's index number is stored in the focus' cloud array.

The next closest point is then checked to see if it is suitable for inclusion within the
cloud. After the first point, the angular position of each point relative to the focus is checked against the angular position of each of the previously selected neighbors. If

A .
This array of distances is sorted from smallest to largest using a Quicksort function 25 the point lies within a wedge of +180/W, centered on the focus, from any of the previously selected cloud members, then that point is rejected. N is the number of neighboring points which form the cloud. 5. The routine repeats step 4 until the set of N neighboring points has been selected.
We found that step 4 was necessary for a stable and accurate solution procedure. Without this step, the points of the cloud could be biased toward one side of the focus, resulting in an illconditioned matrix.
It is possible for the above search procedure to pick points that are on the other side of the body from the cloud's focus. In practice, this never happened with the interior points, but it did occasionally occur when selecting points for clouds on the body's boundary. For boundary points on the body's surface, the following steps were performed before the above procedure was initiated: 1. The two closest boundary points, one on either side of the focus, were first added to the cloud. 2. The plane tangent to the body at the focus was computed.

3.
Only those interior (non-boundary) points that lay "above" this tangent plane, where "above" means in the direction of the outward normal, were included in the above search.

Cloud Equations
The equations presented in this subsection are a repeat of the equations presented in Section 2.5. They are repeated here to facilitate the following discussion. Within each cloud of points, we modeled the functional variation of the unknowns in the governing equations with a 2nd-order Taylor series expansion about the cloud focus, i.e., where (xoy yo) are the coordinates of the focus, h and k are given by h = xxo and k = yyo, and the subscripts x and y denote partial derivatives in the usual notation. matrix is ill-conditioned or singular, SVD can detect the problem and, in some cases, can provide a solution. Within a cloud, the [ A ] matrix is constant unless the cloud's points move or change. Therefore, this matrix is formed and its decomposition calculated and stored at the beginning of the solution procedure. Derivatives are then easily recomputed whenever the [ f ] matrix changes. In the example problems described in Section 3.3, nine points were used for the interior clouds and six points for the boundary clouds.

Boundary Conditions
Rather than use ghost points to impose boundary conditions, we decided to incorporate them directly into the cloud equations. We have examined only two types of boundaries, far field and solid surface. Far field boundary conditions were imposed by fixing the values of the unknown field variables. At wall boundaries, we examined two different methods for imposing boundary conditions. The normal velocity condition at a solid surface is expressed as where A is the unit normal to the surface, V is the fluid velocis vector, and V , is the surface velocity. For the remainder of this discussion, we assume that V i = 0. For a two dimensional problem, this reduces to This equation can be rearranged to give If the function f in equation (3.1) is the velocity potential, 9 , then superimposing the above boundary condition on equation (3.1) gives (3.7) or The first two columns in the h2 k2 Neither method seems to offer an advantage in either speed or accuracy over the other. The second method has the advantage of being cleaner and easier to program than the first.

Solution Procedure
The solution procedure used for the example calculations given below follows the procedure described in Section 2.3 on page 14. The solution procedure consists of two nested iteration loops. The outer, or global, loop spans the computational domain and uses a Gauss-Seidel scheme to reduce the global residual of the governing equation. We found that a Gauss-Seidel scheme converges two to three times faster than a scheme using Jacobi iteration. At the beginning of each global loop, the values of the velocity potential at each point in the domain are updated, the set of derivatives computed, and the local residuals from the governing equation sorted from largest to smallest. The inner, or local, loop is then begun, starting at the point with the largest local residual and proceeding -in sorted order to the smallest. The inner loop calculations are performed for each cloud of points to find the value of the velocity potential at the cloud focus that temporarily brings the local residual of the governing equation to 0 & E. For the test cases discussed below, E = 1 x~O -~. Smaller values of E were tried, but did not improve the overall solution accuracy. These local calculations were performed using Brent's method. This is a modified secant search routine that is guaranteed not to diverge. If the routine detects any divergence in the search, it switches to a bisection routine. Timings showed that this method converged two to three times faster than pure bisection. The actual C function for Brent's method was taken from reference 10. Attempts to use a Newton-Raphson scheme were not successful. After the calculations for the last cloud are completed, the local and global residuals are recomputed. The calculations are terminated when the global residual is reduced by five orders of magnitude from the global residual of the starting solution. If the solution is not converged, another global iteration is performed.

Numerical Examples
In this section, the gridless numerical method is applied to simulate steady, incompressible, potential flow around simple two-dimensional cylinders. The inviscid flow field is computed for two geometries: a circular cylinder and an elliptic cylinder of aspect ratio 5: 1. We first describe these model problems and their exact analytical solutions. We then give the numerical results of the gridless method for the model problems and discuss the types and magnitudes of numerical errors discovered. The numerical results presented here were generated with and without application of the coupled linear operator described in Sections 2.4 and 2.5. Comparison of the gridless method with and without the coupled linear solutions shows that the basic gridless scheme closely approximates the optimal coupled solution for a given spatial resolution. Overall, the numerical results of the gridless method are in very good agreement with the known exact solutions for the circular and elliptic cylinders.

Model Potential Flow Problems
This section describes the exact analytical solutions to two simple model flow problems. To define the model problems, consider steady, incompressible, potential flow around an arbitrary two-dimensional cylinder with no circulation. The first model problem is to compute this ideal flow field for the circular cylinder of radius one, while the second model problem is to compute this flow for an elliptic cylinder of aspect ratio 5: 1. Recall that these idealized flows are described by Laplace's equation subject to simple boundary conditions on the surface of the cylinders and at infinity, i.e., V2@ = 0 V$ -A = 0 on stationary solid bodies V$ approaches the freestream velocity as r + 00 For both model problems, the gridless method is used to compute the velocity potential $ around the cylinders as specified by these equations. The results were used to compute the pressure coefficient on the surface of the circular and elliptic cylinders. The computed velocity model problems.
circular cylinder with no circulation is given by -U potentials and pressure coefficients are then compared with the exact analytical solutions of the As is well known from classical aerodynamics, the complex velocity potential around a (3.10) 2 where 6 = x + iy , U is the freestream velocity, and A = a U . Here a is the radius of the circular cylinder. All the numerical results in this report assume a = 1 and U = 1 . Recall, that the complex velocity potential is defined by where @ is the velocity potential, and w is the stream function. Differentiating the complex velocity potential with respect to 5 then yields the complex velocity. Therefore, equation (3.10) implies that the exact solution for the complex velocity around a circular cylinder with no circulation is given by I.
The exact analytical solution for potential flow around an elliptic cylinder can be obtained in a similar manner by transforming the complex potential obtained for the circular cylinder via conformal mapping. Let z(6) denote a conformal mapping which carries the unit circular cylinder onto an elliptic cylinder. It can be shown that the Joukowski map defined by z(6) = ( + 3 5 (3.13) with q = 2/3 will map the unit cylinder onto an elliptic cylinder of aspect ratio 5: 1. Now, recall that z maps the velocity potential of the circular cylinder to the velocity potential of the elliptic cylinder. So, given the velocity potential for the circular cylinder, F ( c) , the velocity potential for the elliptic cylinder becomes (3.14) where 6 = z-'(() is the inverse of equation (3.13). Then differentiating equation (3.14) with respect to 5 and applying equation (3.13) yields the complex velocity for the elliptic cylinder (3.15) where W ( 6 ) is the complex velocity around the circular cylinder and (3.16) .
Once the complex velocity has been determined, the pressure coefficient C, on the surface of a given cylinder may be computed from

Numerical Results
To test the gridless numerical method, several computational experiments were conducted. This section presents numerical results generated by these experiments for three versions of the gridless scheme for both model flow problems. The first gridless scheme considered is referred to as the basic scheme and applies the algorithm defined by steps 1 through 5 on page 16 of Section 2.3 without modification. The second and third gridless. schemes are referred to as coupled-schemes one and two and modify how the algorithm approximates the derivatives and the differential equation on each local computational cell. Coupled-scheme one computes the least squares solution for the For all of these calculations, the problem was initialized by assuming that the velocity potential was equal to the freestream velocity potential at each field point. The potentials at the outer boundary points were set equal to those of the exact solution potential to satisfy the boundary condition at infinity. This was done so that any errors in the solution would be due to the gridless numerical method and not attributable to boundary conditions. The calculations proceeded until the global residual was reduced by more than five orders of magnitude from the global residual of the starting solution, i.e., (3.18) The global residual is defined as the root-sum-square of the local residuals at each field point. The local residual is the deviation of the governing equation from zero, i.e., Since we know the exact solution for these test cases, we have also defined a local error, EL, as the difference between Q exact and @ calculated The global error, EG, is defined as the root-sum-square of the local errors I N (3.21) where N is the number of points in the computational domain.
.   Here, the zero circular cylinder angle is aligned with the X-axis. Since the gridless scheme yields a symmetric solution about the X-axis, the angle can be measured either clockwise or counter-clockwise. Figure 3.5 shows that the gridless scheme generated a pressure coefficient in very good agreement with the exact analytical result.

Circular Cvlinder Calculations
The overshoot in the global velocity potential error before approaching a steady-state value (Figure 3.3) is believed to be a numerical oddity with no real significance. During the code development, other solution techniques were examined. For some of them, the solution error approached the same asymptotic value from above, i.e., with no overshoot. We believe that given this point distribution and second-order derivative approximation, this is the best solution that the current method can achieve. In order to improve the solution accuracy, point adaption and/or a higher-order derivative approximation are needed, at least in the vicinity of the stagnation points.
A second test calculation was made for flow over the cylinder with the number of field points doubled in both directions, i.e., 160 by 82, for a total of 13,120 points. The convergence characteristics are shown in Figures 3.6 and 3.7. The distribution of the local error around the cylinder is shown in Figure 3.8. The increased spatial resolution reduced the global error by a factor of approximately 4.5. However, the number of global iterations required to achieve the same level of residual convergence increased by almost a factor of approximately 3.75. The computer time required to achieve this solution increased by a factor of approximately 16 over the previous solution, indicating that the current implementation gives an N2 method. We have made no effort to optimize the efficiency of the current code.
We also compare the basic gridless scheme with the two coupled schemes for potential flow around the circular cylinder. Figure 3.9 compares the global residual of the basic scheme with the coupled gridless schemes one and two. Both coupled schemes yield the same global residual as a function of iteration number. Moreover, the coupled schemes produce a residual which was slightly smaller than the basic scheme. For a given approximation of the derivatives and the differential equation, the least squares solution generates the approximation with the smallest residual. This suggests that the basic gridless scheme closely approximates the optimal least squares solution for a given spatial resolution. Figure 3.10 shows the global errors computed with all three schemes. They all produced the same global error as a function of iteration number, since all three schemes performed a local iteration on each computational cell by calculating the root of an equation. Recall that the local iteration process is defined in step 4 of the basic gridless method. The differences between the coupled and uncoupled gridless schemes were so small that the local iteration yielded the same global error. One main conclusion of Figures 3.9 and 3.10 is that the basic scheme produces numerical approximations which are essentially identical to the optimal coupled solutions.

Elliptic Cvlinder Calculations
Calculations were also performed for potential flow around an elliptic cylinder of aspect ratio 5: 1. Figures 3.1 1 and 3.12 show the field point distributions used in these calculations. For this example, the field point locations were obtained through the conformal mapping of the points for the circular cylinder (Figure 3.1). As for the first model problem, these calculations used 80 points circumferentially and 41 points radially for a total of 3280 field points. Figure  3.13 shows the convergence properties for this calculation using the basic gridless scheme. As with the circular cylinder, the global error approaches an asymptotic value while the global residual continues to decrease. Figure 3.14 shows the spatial distribution of the local error. The largest error observed is less then six-tenths of a percent of the exact solution and occurs at the flow stagnation points. Figure 3.15 shows the pressure coefficient on the surface of the cylinder. The zero elliptic cylinder angle is aligned with the X-axis. As in the circular case, the gridless scheme yields a solution symmetric about the X-axis, hence, the angle of Figure 3.15 can be measured either clockwise or counter-clockwise. The plot of the pressure coefficient demonstrates that the basic gridless scheme generated numerical results in excellent agreement with the exact analytical solution.

Overview of Tasks
In order to assist with the development of gridless CFD techniques, we were asked to analyze several aspects of the problem. The ideas for our tasks were based on work previous finite-element researchers performed using h and p techniques for improving the s o i u t i~n .~~* In general, h versions of the finite-element method fix the degree of the element and convergence is achieved by mesh refinement. The p version fixes the mesh and convergence is obtained by increasing the degree of the element. The gridless technique also appears to lend itself to a form of h-p adaption. The tasks summarized in this chapter focus on investigating the effect that the order of the Taylor series has upon the solution convergence and accuracy, analogous to the p finite-element techniques. Subsequent tasks, to be reported in phase two of our project, focus on h adaption techniques.
The first task required the development of a generalized recursion relationship for an nthorder, two-dimensional Taylor series. This was then incorporated into a C-language program. The routine requires the user to define the number of terms in the Taylor series. The second task was to assess the impact of additional terms in the local Taylor series approximation on convergence and accuracy of the solution of a two-dimensional equation, such as Laplace's.

Two-Dimensional Taylor Series
An expression for the Taylor expansion about a point (xo, y o ) is given by9 where N is the order of the Taylor series and the remainder term, R, + , is given by The point (5, q ) is somewhere on the line segment joining the points (xo, y o ) and ( x , y ) .
In this expression h and k are given by h = xxo and k = yy o . For example, a third-order Taylor series expansion is written as: where the subscripts x and y denote partial derivatives in the usual notation. In order to write a recursive relation for the two-variable Taylor  where (a + b)" is given by The binomial coefficient [ ] is given by (nL n! Both the binomial coefficient and the factorial term can (4.5) be easily calculated using canned computational routines, such as those found in Numerical Recipes in C." A C-language program was written in order to check the higher-order Taylor series algorithm. A function whose derivatives could be determined analytically, given by (4.7) was used to test the C-code. A surface plot of the function is presented in Figure 4 The numerical evaluation of the derivatives depended on two different parameters. The first parameter of interest was the number of points necessary to form a "cloud" surrounding the evaluation point. Clearly, in order to obtain higher-order derivatives, more points were needed. However, we were also interested in the effect the number of points in the cloud might have on the solution. For example, a third-order solution requires a minimum of 9 points, however it is possible to form a cloud with more than 9 points. The second parameter of interest was the cloud geometry. Circular, rectangular, and irregular clouds were formed.
In the first investigation, all points were equally spaced on a circle or circles, centered at (0.5,0.5) or (0.9,0.8). For example, for a cloud of 20 points, 10 points were equally spaced on a circle of radius 0.01 and another 10 points were equally spaced on a circle of radius 0.02 (see The results of the first investigation for the function are shown in Table 4.1 for the evaluation point (0.5, 0.5) and in Table 4.4 for the evaluation point (0.9, 0.8). Second, third, fourth, and fifth-order terms in the Taylor series were examined. Clouds contained either 5, 9, 14, or 20 points. The cloud geometry consisted of a single circle centered about the evaluation point with radius of 0.01 (denoted R = 0.01 in the tables) or families of circles centered about the evaluation point with increasing radii, as indicated in the tables.
We offer several comments about the results. We observe that the selection of cloud geometry greatly influences the solution, especially for the evaluation point (0.9, 0.8) which is best results for points selected to lie on either two or three circles. Spreading the points out over too many circles did not give as good of results for the fifth-order solution. Clearly, there is an optimum cloud geometry which is dependent on the function and evaluation point. In general, the fifth-order results are poor for a fifth-order Taylor series, the fourth-order results are poor for a fourth-order Taylor series, and so on. We noted that one must select a Taylor series order that is at least one order higher, preferably two orders higher, than the order of the derivative of interest. Since most governing equations of interest contain only first and second order derivative, it appears one could select a fourth-order Taylor series and obtain good estimates for the first and second derivatives.
The next investigation was to form several "clouds" which better resemble those that might be formed using a rectangular-based mesh. Again, the evaluation points were taken as (0.5, 0.5) and (0.9, 0.8). A rectangle is formed which surrounds the evaluation points. The midpoints of the sides of the rectangles are located a distance of 0.01 from the evaluation point.    I  '  I  '  I~I  '  I  '  I  '  I  '  I  '  I  By selecting a rectangular clustering of points, we ensured that not all points lie on a circle with a fixed radius. Hence, all results give fairly good approximation to the lower-order terms in the Taylor series. In fact, it is interesting that the 8-point rectangle which provides only a 2ndorder approximation to the Taylor series, gives reasonable results. Clearly, if one is interested in only first and second derivatives, a fourth-order solution using a double rectangle should give very good results.
Finally, we considered a more pathological case of cloud point distribution. In particular, suppose grid adaption or mesh refinementkoarsening has occurred and the points have clustered to one side of the evaluation point. Again, two evaluation points were used: (0.5,0.5) and (0.9, 0.8). In this case, seventeen neighboring points were selected to evaluate the derivatives and are shown in Figure 4  Order of Taylor Series Taylor series gives better results than the second-order Taylor series. These results are very consistent with the results found by finite-element researchers for p-adaption. Higher-order terms can be used to improve solution accuracy. However, these results suggest there is an optimum higher order. Clearly the fifth-order results are no better than the fourth-order results. The fifth-order Taylor series approximation may have provided too many degrees of freedom for estimation of the function. We also note that we would be hard pressed to answer the question: "What cloud geometry gives the best results?' The cloud geometry and number of points in the cloud affect the solution, but these results do not indicate which combination is best or worst. The irregular cloud geometry gave reasonable results, which is an important result as we proceed in looking at h-adaption where the points will be moved/added/deleted and irregular clouds may be formed. Finally, the function had a less steep gradient at the (0.5, 0.5) evaluation point and the numerical simulation did a better job of predicting the value of f , . As the function steepened, the scheme did not predict the f xx term as well.

Higher-Order Solutions of Laplace's Equation
The 2-D Laplace's code described in Chapters 2 and 3 had to be modified in two places in order to permit higher-order terms: in the formation of cloud geometry and in the formation of the [A] matrix used to evaluate derivatives (see equation (3.2) on page 26). First we noted that in order to include third-order Taylor series terms, at least 9 points must be used to form a cloud, similarly for fourth-order Taylor series at least 14 points must be used, and for fifth-order Taylor series at least 20 points must be used to form the cloud. Since we were interested in comparing higher-order results with the second-order results described in Section 3.3, we utilized the same point distribution, which results in cloud geometries similar to the sketch shown in Figure 4.

4.
In order to develop a robust code which would allow a user-defined number of points for a cloud, we began with the existing cloud formation algorithm. The best way to summarize the this algorithm is to imagine the evaluation point as the origin of a local x-y coordinate system. The existing algorithm selects eight points which are fairly uniformly distributed in all four quadrants of this x-y coordinate system. This is accomplished by calculating the distance all points are from the evaluation point, sorting the distances and then selecting the point which is closest to the evaluation point. The closest point is the first member of the cloud. The point which is second closest to the evaluation point will be allowed in the cloud formation provided it is not located within a pie-shaped wedge of 45" (k22.5") which emanates from the origin and contains the first-closest point. Similarly, the third point in the cloud formation will be third closest to the evaluation point provided it does not lie in the pie-shaped wedges which surround the previous two points. Should the third point lie within a pie-shaped wedge, it is rejected and the next closest point is selected as a candidate. This algorithm effectively finds eight cloud members which are fairly uniformly distributed in all four quadrants which surround the evaluation point.
Since we were interested in higher-order terms, we always needed to form clouds with more than eight points. We used the existing algorithm to select the first eight points. We then return to the list of sorted distances and select the next closest to the evaluation point which has not already been used in the cloud geometry. We form a 22.5' (k11.25") pie-shaped wedge centered around this point, return to the sorted list of distances and select the next closest value. If this value is not within the 22.5' pie region, it becomes a member of the cloud geometry, we form another 22.5" pie region surrounding this point and then we select the next closest point which does not lie in the previous two points' pie-shaped wedges. This process is repeated until the required number of points have been selected. Again, we are trying to ensure the cloud contains points which are more uniformly distributed in all four quadrants surrounding the evaluation point.
The second modification to the 2-D Laplacian code was to calculate an appropriate [A] matrix which would contain the higher-order terms. We were able to use the C-function used to generate the results in Section 4.2. The user has to specify the order of the Taylor series, the number of points in an interior point cloud, the number points in a wall cloud, and the number of derivatives.
Results of applying higher-order terms to the solution of Laplace's equation are described below. We utilized the grid shown in Figure 3.11, which contained 3280 points, and solved Laplace's equation for flow over a 5:l ellipsoid. The boundary conditions were the usual for potential flow, no flow normal to the boundary. First, we ran the existing 2-D code which was second-order using 8 points in a cloud. Then we ran a number of cases for higher order solutions.
The results are disappointing. They indicate that by using a higher-order Taylor Figure 4.7 is the largest local error value in the domain. The second-order Taylor series calculation using eight nodes gives the best results, although the error reaches a minimum and then begins to grow before the convergence criterion is reached.
We observed that the largest local error in all cases occurred on the wall. So we tried one additional case. Without too much effort, we were able to modify the 2-D Laplacian code to allow one order of solution on the wall and another order of solution in the interior. This gave better results, as shown in Figure 4.8. We were able to continue to reduce the error as the residual was reduced when we selected third-order on the wall and second-order in the interior. Finally, it is important to note that a number of cases would not converge. For example, the case of 19 points on the interior, 16 points on the wall converged, however, for cases with the same order of Taylor series and less than 19 points on the interior, the solution did not converge.
These results suggest that we would probably have reduced error for higher-order terms if we would select a region "close" to the wall and incorporate a higher-order solution. This is not surprising. The flow field is changing rapidly near the wall and requires some form of adaption. In the far field, the flow is essentially freestream which has no higher order behavior. By trying to insist on higher-order functions for freestream behavior we are giving the solution too many degrees of freedom.

Summary and Recommendations for Future Work
We developed a generalized recursive relationship for an @-order, two-dimensional Taylor series which was incorporated into a C-program where the user is allowed to define the number of terms in the Taylor series. We used this code to evaluate behavior of higher-order terms on solution accuracy. We varied both the number of points in a particular cloud and also the cloud geometry. In general, we found the most accurate results were obtained if the user selected terms in the Taylor series which are two orders higher than the highest derivative to be evaluated. However, the accuracy is a strong function of cloud geometry and number of points in a cloud.
We then incorporated the higher-order Taylor series code into a Sandia developed code which solved the two-dimensional Laplace's equation. We evaluated the effects of the addition of terms in the Taylor series approximation on the convergence and accuracy of solutions. We found that simply fixing a higher-order Taylor series everywhere in the flow field did not give better accuracy because in the far field the higher-order terms provide too many degrees of freedom for a flow field that is essentially freestream. Since the largest error occurred near the body, we allowed the solution to be higher order on the wall and lower order in the interior.
These results provided close to the same accuracy as the simple second-order Taylor series results and better convergence properties. The error continued to diminish as the residual diminished.
Future efforts will be focused on incorporating both h and p adaption into the twodimensional Laplace's equation solver. The results found to date suggest that if we could select a region of points sufficiently close to the wall and apply a higher-order Taylor series, we would probably obtain better accuracy than the second order results. This will require a major rewrite of the two-dimensional code. We also plan to investigate the effects of h adaption on the solution accuracy. fw 0.

Introduction
Adaption is a procedure that complements and enhances the capability of numerical methods. The goal of adaption is twofold: provide a more accurate solution and minimize the computational effort. A more accurate solution will globally approach the true solution while capturing the important local features. Computational effort is minimized when calculations that provide negligible contributions to the accuracy or resolution of the solution are eliminated.
During the initial process of applying most numerical methods, a preprocessing procedure of "seeding" the probIem domain, whether by meshes or points, is required. Generally, this initial "seeding" is not an optimum one regarding the criteria of accuracy and computational effort. Hence the introduction and utility of adaption.
Adaption techniques can be categorized into three distinct types: movement, refinemenu coarsening and higher-order approximations. All three types have been investigated and utilized in traditional grid-based solvers. Adaption by movement is the most mature of the three, with the spring analogy perhaps the most widely used implementation. In this technique, the mesh structure is moved according to a local spring constant, which is itself a predefined function of the solution field. Adaption by mesh refinementkoarsening, sometimes referred to as hadaption, will locally refine or coarsen the mesh according to a programmed algorithm. Similarly, adaption by higher-order approximations, sometimes referred to as p-adaption, will locally modify the order of the approximating function to an appropriate level according to a programmed algorithm.
Extension of these standard adaption techniques for gridless methods is straightforward. In fact, adaption for gridless techniques should be easier since points are geometrically simpler structures than cells. The primary consideration for adaption, regardless of the type, is the composition and structure of the algorithm that drives the adaption procedure. This logic path must consider relevant issues such as: when to adapt; the proper forcing function for adaption; the order andor concurrency of adaption types; and perhaps most importantly, any imposed limitations on the adaption procedure. This is the primary focus of the research work in adaption for this LDRD.

Goals for Adaption
The goals for the adaption algorithm under development are stated as follows: 1. Achieve a level of semi-autonomy, allowing the user to intervene only to prescribe a general purpose accuracy parameter; 2. Incorporate generality, permitting the solution of a wide class of problem sets; 3. Achieve a state of robustness, in which the user-defined level of accuracy can be accomplished, regardless of the initial point seeding; 4. AIlow generdity of the adaption algorithm technology to accommodate multidimensional problems.
A general methodology was formulated and developed to achieve these goals. The basic rule set for the adaption algorithm is outlined in this chapter.

The Adaption Algorithm
Determination of the appropriate parameters on which to adapt is the primary consideration in constructing the adaption algorithm. A sampling of these parameters include: the local error from the "true" solution (generally not available), the approximate local error associated with truncation of the approximating function (which usually is dependent on the independent spatial variable), the local residual (generally a function of the rate of convergence which may not be directly associated with the error from the true solution), gradients of dependent variables and their derivatives, the condition number of the local neighboring point approximation matrix (which may affect accuracy through numerical round-off) and the interaction of the spatial resolution and grouping of points related to the approximating function. In some form, all of these parameters must be considered.
Because of their generality and their availability with this gridless technique, the current adaption algorithm utilizes the dependent variable, the gradient of the dependent variable and the derivative of the gradient. This implementation is a slight departure from adaption in commercial solvers, which primarily use the gradients of the dependent variables. By allowing differences in both the variable magnitudes and their first and second derivatives (changes in the slope and curvature of the solution) to influence the adaption of points, the methodology is both generalized and moderated.
The mathematical construction of the adaption parameter is shown in equations (5.1) and (5.2). The adaption parameter comprises two spring coefficients that are essentially linear sums of the uniformly weighted, normalized, absolute values of the function and its first and second derivatives. These linear spring coefficients are utilized in multiple roles in the hierarchical process of the adaption algorithm that is described below.
The first step in the adaption process is to identify the criteria under which adaption would be advantageous to the solution. Certainly, prior to engaging any form of adaption, the solution should have achieved some level of convergence such that the dependent variable's values and derivatives are reasonably accurate. This condition ensures that adaption will not contribute to solution divergence. Ideally, the local and global solution residuals should have converged to a specified level (e.g., two orders of magnitude on a normalized scale) and demonstrated a decrease in the rate of residual convergence.
Assuming the solution has converged to an appropriate level, the first adaption procedure implemented is point movement. Adaption by point movement is governed by spring forces, where a generalized spring force is given by equation (5.3) (the term Fi-~ refers to the force acting on point "i" by point '5-1"). The actual point movement is governed by equation (5.4), I c c which incorporates the force vectors and is normalized by the sum of the spring coefficients. Point movement continues until the discrete spring force imbalance declines below a user specified tolerance, as provided by equation (5.5).
The process of adaption by point movement is an iterative procedure. That is, multiple passes of each individual point may occur during a single adaption event, without the solution being updated during this process. This iterative procedure is necessary due to the discrete spring approach adopted for the system point movement algorithm. However, because of this iterative point movement approach, several constraints must be imposed on the process, as described below.
The first constraint imposed on this adaption technique is to limit the maximum distance that an individual point may move during a single iteration. This constraint is necessary for two reasons: it would be undesirable to have points leapfrog neighboring points and the algorithm that updates the dependent variables' value and derivatives (during the adaption procedure) is based on linear interpolation, which is relevant for short distance moves only. The second constraint is similar to the first, in that the cumulative movement of individual points during the entire iterative procedure of a single event of adaption by point movement is limited. Typically, an individual point is limited to a single iteration movement distance of 5% of the distance to the nearest neighboring point, while the cumulative constraint is restricted to twice that of the single iteration distance.
A third constraint imposed during the point movement algorithm is necessary to counteract the potential oscillatory nature of a system of individual, but coupled springs. That is, point movement is actually governed by equation (5.6), which is similar to equation (5.4), but incorporates a numerical viscosity dampening factor, a. The dampening factor ranges from 0 to 1, with a value of 0.9 to 0.95 typical. It is only invoked when an individual point movement reverses direction during a single adaption event, thereby dampening the oscillatory movement. A single constraint is imposed in conjunction with this refinementkoarsening adaption procedure. That is, points that are nearest neighbors to a point that has been refined or adapted cannot themselves be refined or adapted during this particular adaption cycle. This condition ensures that adaption by point movement has an opportunity to reposition the points appropriately, without the potential for extraneous local refinement or coarsening.
The aforementioned procedures and constraints reflect the 1-D point adaption scheme as currently implemented. This algorithm addresses the first three goals desired for the 1-D adaption implementation. It achieved the semi-autonomy desired in goal one, while providing a broad based from which to verify the generality of problem scope and robustness desired in goals two and three.
In summary, the current adaption scheme incorporates both point movement and point refinement and coarsening procedures. While the gridless code possesses the capability to globally utilize higher-order approximation functions, the extension to local p-type adaption has not been implemented, due to unresolved issues concerning the appropriate combination of pversus h-type adaption and the coexistence of neighboring points with an incongruent order of the approximation function.

Adaption Algorithm Performance
The performance of the adaption algorithm is now assessed through application of several test problems. These test problems are posed in the form of linear and nonlinear ODE'S. Each differential equation considered has a known closed-form solution, which allows an accurate assessment of the numerical gridless solution technique and the assistance provided by implementation of the current adaption algorithm.
The first problem investigated is the simple linear, 2nd order ODE, given by eq, (5.7): a'f =   Since this particular ODE is well behaved, it does not pose a severe challenge to the adaption algorithm in capturing large gradients in the solution. However, two other aspects of the adaption routine can be investigated with this simple ODE. The first considers the movement of points that have been intentionally clustered. If the adaption by point movement is configured correctly, the clustered points should be repositioned to a more uniform spacing. For this test case, the adaption by refinementkoarsening has been disabled, even though the latter type of adaption would be more computationally effective through judicious implementation of point deletion and insertion.
The result of the adaption by point movement test case is displayed in Figure 5.2. The top line shows the initial clustering of the points. The middle line demonstrates the final point locations subsequent to repositioning by movement adaption. As observed and desired, adaption by movement has repositioned the initially clustered points to a more uniform spacing. The bottom line provides an "ideal" point positioning for this method of adaption. This ideal distribution represents a point seeding in which the analytically determined kAx value is identical for each point. As observed, the adaption method has distributed the points adequately, when compared with the ideal distribution.
The second test case considers the robustness of adaption by refinementlcoarsening. Specifically, in solving the aforementioned linear ODE, two sets of initial point distributions were utilized. Within the tolerance and accuracy of the user-defined upper and lower bounds parameter, the adaption algorithm should be sufficiently robust to determine the approximate optimum number of points for this solution, regardless of the initial point seeding. In the first To further demonstrate the robustness and overall capability of the gridless technique, two additional differential equations were solved. For both of these solutions, adaption by point movement and refinementkoarsening were utilized.
The next example is a slight modification to the previous linear 2"d-order ODE, equation (5.7), given by (5.9) Utilizing the boundary conditions f(0) = 1 and f ( 1) = 0 results in the solution shown in Figure 5  where the wave speed is given by c = (f2 + f l ) / 2 . The variable $ asymptotically approaches the limits f l to the right and f2 to the left. By selecting the parameters f = 1 and f = -1 , the wave speed is set to zero, which simplifies the sohtion and eliminates the time dependency.
As a result of the choice of f, and f2, the boundary conditions are established as @(-1) = 1 and @( 1 ) = -1 . The Reynolds number for this solution is arbitrarily established at 50. The comparison between the numerical gridless solution and the analytic closed form solution is presented in Figure 5.4. As observed, the gridless technique has captured the true solution for this nonlinear ODE. However, it should be noted that the tolerance applied to the point movement algorithm was relaxed slightly, which stabilized the solution in the region near the discontinuity and reduced convergence times by a factor of three to four.

1-D Adaption Summary
Several example problems, including both linear and nonlinear ODE'S, were solved using the gridless technique developed in this LDRD project. Incorporated within this gridless technique is the capability to adaptively move points, which can enhance solution accuracy and minimize computational expense. Point adaption was demonstrated using both movement and coarsening/refinement techniques. The adaption algorithm was demonstrated to be autonomous and robust, which achieved the goals established for this task in the LDRD. Further, the unique formulation of the adaption parameter that utilizes local differences in the variable's value and

The Adaption Methodology in Multi-Dimensions
The extension of the l-D adaption methodology to 2-D (or higher dimensions) is not difficult. However, it encompasses additional degrees of freedom for implementation that generally manifest as more complicated limitations on the algorithm. The fundamental logic for point adaption is independent of the spatial dimension. That is, point movement precedes point refinementkoarsening and its complementary technology of higher order approximations. The user defined accuracy bounds defined for l-D are still applicable, as are the limitation imposed on the adaption procedures outlined for l-D. The departures from the l-D adaption methodology are outlined in the paragraphs below.
Regarding point movement, the consequence of the multi-dimensional spring system is not l-D translation, but a linear translation that incorporates the contributions of all the surrounding points of the cloud group. This requires the simultaneous solution of the cloud group spring forces to achieve the point adaption equilibrium during a single iteration of the movement adaption cycle. By necessity, the point distribution of the cloud group surrounding the focal point will be important to minimize directional dependence of the point adaption through movement.
Adaption by point deletion and insertion is less constrained than the analogous l-D case, however the basic rule set is still applicable. That is, adaption by point movement continues until a prescribed level of translation equilibrium is achieved, at which time the adaption by point deletion and insertion may be considered. As in the 1-D case, if the adaption parameter (specifically the linear sum of the difference in point values, first, and second derivatives detailed previously) exceeds a user-defined accuracy value for all the neighbors in a cloud group (or a majority of the neighboring points, depending on the degree of equilibrium achieved), then a point insertion process may be implemented. Conversely, if the adaption parameter declines below a user-defined lower bound for all the neighbors in a cloud group, then the point may be deleted. For the point insertion algorithm, the following procedure is proposed. Insert a point between the focal point and its neighbor with the highest adaption parameter magnitude. Continue this procedure with all of the other neighboring points. It may be preferable, however, to limit the number of points added by skipping every other neighbor point during a circumferential search. Other limitations, such as distance parameters for validity of the Taylor series approximation may be important and should be investigated either mathematically or empirically as a future research topic.

The Adaption Methodology for Time-Accurate Problems
The implementation of adaption for time-accurate problems is very complex. Whereas an optimum point location can be determined gradually during the iterative process associated with the solution of a steady-state problem, a time-accurate problem requires reasonable point placement at each successive time step. This necessitates an adaption algorithm that can rapidly react to the transient, even though it will inherently lag the solution. It would be more desirable to produce an algorithm that is not reactive, but predictive in nature. This might be accomplished utilizing a predictor-corrector scheme that is intrinsically coupled to the timeaccurate solver. By tracking trends in the point adaption, points could be more optimally positioned for future time steps. This task represents another future research topic, but certainly not one that represents a barrier to the generality of the adaption methodology under development.

Shell Elements
Frank Mello

Introduction
Recent advancements in computing technology have broadened the horizons of simulation capabilities by increasing the size and complexity of the models that can be considered. Frequently these simulations span different time scales as well as regions with different governing equations. Such is the case in coupled fluid structure modeling. A new technique for solving fluid dynamics equations on a random unstructured collection of points is being evaluated. In support of this gridless fluid modeling LDRD, plane strain and axisymmetric shells have been added to PRONT02D. The 3D shell element developed by Belytschko, et al."? l 2 has been specialized for 2D. The reduction in dimension dramatically simplifies the implementation. Rotation reduces to a vector of fixed orientation and the need for tracking nodal and element coordinates i s eliminated. Further, there are no hourglass modes to control in 2D shells. The same commands used to describe shell element attributes in PRONT03D are available in 2D, including shell thickness scaling and layered shell definitions. Much of the text in this chapter is borrowed from documentation of the 3D features.
The coupling between fluid calculations and structural calculations will occur by passing boundary condition data between the codes which will run simultaneously. This strategy has been effective in linlung PRONT03D and CTH. where T is the stress tensor, p is the material density, u is the acceleration vector and b is the body force vector. When a body is decomposed into a finite element mesh, the momentum balance equation is enforced in a weak sense over each element volume. The assembled weak form, after an integration by parts is then; where 6zii is an arbitrary virtual velocity field in the shell body. The first term in this expression is the rate of work of externally applied tractions and/or displacements. The following terms correspond to the rates of strain energy, work of the body forces, and work of the inertial forces, respectively. Since PRONTO uses a lumped mass approach, the above relation reduces to a set of uncoupled equations describing the motion of the nodal points. The acceleration of the Jth nodal mass may then be expressed as In shell elements the nodal points carry rotational translarions. The kinematical assumptions discussed in degrees of freedom in the following section  The simplicity of this single equation for rotation leads to the primary reduction in complexity for 2D shells. Since the only principal direction which admits rotation always has a fixed orientation the need to track nodal rotation coordinates is eliminated.

Shell Kinematics
Mindlin shell theoryI5 assumes that sections through the shell thickness remain straight, although they may rotate. Rotation of the section allows the element to model transverse shear strains. Since displacements are assumed to vary linearly through the thickness, the velocity at any point may be expressed in terms of midplane values as u = u + + e x o k (6.6) in which ri is the velocity of a point in the shell body, u is the velocity of the point on the midsurface which has the same normal, and ok is the rotational velocity of the section. The coordinate of the point through the shell thickness is 2 and the unit vector along the nonnal is e .
Since this is a 2D shell formulation, the unit vector k , indicating the axis of rotational velocity, is always normal to the plane of the problem.
The components of velocity strains in the corotational coordinate system are (6.7) Substitution of the above relation for zi into the velocity strain definitions yields the strain displacement relations for the shell. where the subscript 1 refers to the direction defined by progressing from the first node of the element to the second. The subscript 2 refers to the thickness direction and subscript 3 refers to the direction into the plane of the problem.

Constitutive Assumptions matrix as;
If the velocity strain components in the corotational system are arranged in a column then the conjugate Cauchy stress components can be written as Furthermore, the shell is assumed to be in a state of plane stress, so 62, = 0 (6.11) These representations are conjugate in the sense that the rate of internal strain energy density, W , can be expressed as w = aTs

Element Formulation
(6.12) The shell element in PRONT02D uses a two noded line segment to model the middle surface of a thin structure (or some other reference surface if an offset is given).

Coordinate Systems
In the three dimensional formulation three coordinate systems are needed. Translational motion is described in the global system. Element calculations are performed in an elementbased system and rotational motion is described in a nodal system aligned with the principle directions of the nodal inertia tensor. Approximating the nodes as spheres (i.e. homogenizing their inertias) is a reasonable approximation for well shaped elements and eliminates the need for a nodal coordinate system.
In the two dimensional framework we use a glob4 coordinate system and an element system. No nodal system is required since the inertia of each node is a scalar (motion is confined to the plane of the problem). Transformations between global and element systems is particularly simple since both have the same third axis.

Element Equations
The shell element is based on a two noded line segment with linear interpolation of reference surface coordinates and of both translational and rotational velocities. Coordinates in the reference surface of the shell are approximated as where $, are the shape functions. Repeated upper case indices indicate summation over the two nodes of the element. Similarly, the velocity of the reference surface and the angular velocity of the shell normal are interpolated as (6.14) using the same shape functions.
an orthogonal set of base vectors as For the case of two dimensional shells, these shape functions may be expanded in terms of (6.15) where 5 is the coordinate in the unit line into which the element is mapped. The base vectors represent deformation of the unit line. represents the stretching of the element and completely describes the element area. For this element, there are no neglected modes and therefore no hourglass control is required. Since we are only interested in quantities at the element center, it is easily seen that the velocities are the average of the nodal values and the velocity gradients are the differences between node 2 and node 1 values divided by the current element length.

Computation of Internal Force and Moment Resultants
As mentioned above the gradient operator for this element reduces to a difference between nodal quantities divided by the element length. So the velocity strains at the reference surface may be expressed as where K is the curvature rate. The constitutive relations operate on these velocity strains to produce stress increments. The total stress is then integrated to element forces and moments. These forces are transformed to global components and assembled. Moments are computed directly from the integration point forces and the offsets.

Hourglass Control
There is no hourglass control for two dimensional shells.

Calculation of Stable Time Increment
The central-difference operator is conditionally stable with the stability limit for a given system with no damping being 2 A t < -~r n l l x (6.17) where mmax is the maximum frequency of the system. Here the maximum frequencies of a two dimensional shell system are conservatively estimated based on the most heavily constrained collection of elements which still admits deformation. Frequencies associated with translations and rotation have been considered. The membrane frequency accounts for the hoop stiffness in axisymmetric problems. Further, the rotational frequency will always bound the transverse deformation case. So the critical time step is based on the largest of the membrane and rotational frequencies.

Constitutive Models
Currently two material models have been implemented, elastic and combined hardening elastic-plastic. More constitutive models will be included as needed.

Summary and Recommendations
A gridless method has been developed for the simulation of coupled fluidstructural interactions over arbitrary bodies. This method uses Eulerian-based points arbitrarily distributed over the computational domain with no formal connectivity as typically required for a traditional grid. Each point in the domain has associated with it a cloud of neighboring points. Derivatives in the governing PDE's are numerically approximated by assuming a functional dependence of each variable within the local cloud of points, where the number of points in the cloud is greater than or equal to the number of unknown derivatives. The governing partial differential equations are then solved locally at each point in the computational domain. This gives a solution to the governing equations without any a priori assumptions concerning grid or element topology, and allows for addition, deletion, or movement of points with no change in soIution algorithms. The gridless method is neither a finite difference, finite volume, nor a finite element extremely versatile in that the method permits one to solve the relevant equations based on the complexity of the problem at hand. Hence, the difficulty of simulating either a potential flow, Euler flow, or full Navier-Stokes flow, and the ability to model compressible or incompressible flows is dictated only by the set of PDE's solved and is easily amenable to a modular code structure. This modularity appears as a simple consequence of the fact that the equations of motion are independent of the point distribution and the solution algorithm is completely grid free. I I approach since differences, metrics, areas, or volumes, are not computed. The methodology is The gridless method had been implemented in a two-dimensional code written in ANSI standard C. Example calculations using Laplace's equation as the model test problem have shown that the current implementation of the gridless method gives reasonable accuracy when a second-order approximation is used to calculate local derivatives. Comparisons with optimized coupled linear methods for the solution of the equations have shown that the current basic implementation is very near optimum, given the point distribution and order of approximations. Accuracy can be improved by using a third-order approximation in regions with larger gradients, such as on the body surface. Point adaption has also shown promise for improved accuracy. However, point adaption has not yet been extended to two dimensions. The optimum balance between higher-order derivative approximations and point adaption has not been determined.
For the structural simulation side of the coupled fluidstructure problem, plane strain and axisymmetric shells have been added to the code PRONT02D. The three-dimensional shell element developed by Belytschko, et al, has been specialized for two dimensions. The reduction in dimension dramatically simplifies the implementation. Rotation reduces to a vector of fixed orientation and the need for tracking nodal and element coordinates is eliminated. Further, there are no hourglass modes to control in two dimensional shells. The same commands used to describe shell element attributes in PRONT03D are available in two dimensions, including shell thickness scaling and layered shell definitions. The coupling between fluid calculations and structural calculations occurs by passing boundary condition data between the codes which i, \ run simultaneously. This strategy has proven effective in linking the codes PRONT031.? and CTH.
More work needs to be done to bring the current gridless code to a production status. One of the more pressing needs is to extend the point adaption scheme to two dimensions. Once this is done, trade-off studies can be run to determine the balance between point adaption and higher order approximations for improved solution accuracy and speed. Work also is needed to improve the code efficiency so that run times do not increase as the square of the number of field points. Improvements need to be incorporated to extend the code beyond the current Laplacian model problem to more complicated fluid equations, such as the Euler and Navier-Stokes equations. Extensive coupled fluidhtructural calculations are needed to ensure that all of the problems of coupling the two codes are resolved. Finally, after the above improvements are made, detailed comparisons are needed between the gridless method code and existing fluid mechanics codes with similar capabilities. Even after the gridless code has been shown to work for all cases, it will not have been proven that this method is more efficient than, say, a fluid code that uses an unstructured grid and regrids as needed to handle a moving structural boundary.