Iterative circle-inserting algorithm CST3D-OC of truly orthogonal curvilinear grid for coastal or river modelling

A geometric method to generate orthogonal curvilinear grid is proposed here. Elliptic partial differential equations have frequently been solved to find orthogonal grid positions, but questions on orthogonality have remained so far. Algebraic methods have also been developed to improve orthogonality, but their applications have been limited to special situations. When two confronting boundary lines of the quadrilateral boundaries are straight, and their positions are known, and we assume that some degree of freedom exists on the other two confronting boundary curves under the condition that the each curve passes through a point, we can assign a set of latitudinal curves in the domain using polynomials. The curves are expected not to fold on their own. The grid positions along longitudinal curves are found by inserting circles between two neighbouring latitudinal curves one by one. If the two curves are straight, the new grid point above the grid point of interest can be found geometrically. This algorithm involves iterations because the curves are not straight lines. The present new algorithm is applied to a domain, and produced almost perfect orthogonality, and similar aspect ratio compared to an existing partial differential equation approach. The present algorithm also can express almost quadrant domain. The present algorithm seems useful for generation of orthogonal curvilinear grids along coasts or rivers. Some example grids are demonstrated.


Introduction
Some governing equations for fluid flow use other coordinates from Cartesian coordinates. When boundary shape is not straight, and still flow behaviour near boundary plays an unimportant role, boundary fitted grid becomes more useful. Flows at rivers or coasts are often horizontally twodimensional, and requires curvilinear grid. Furthermore, if terms for non-orthogonal grid are to be excluded, orthogonality of a curvilinear grid becomes very important, and affects the accuracy of the whole computation.
Methods to generate orthogonal curvilinear grid may be evaluated on three quantities, i.e. orthogonality, aspect ratio, and computational speed of generation. We don't pay attention to the computational speed here. The relative importance of the orthogonality and the aspect ratio depends on the situation or modellers preference, if one value erodes the other value. The uneven aspect ratio has been treated by numerical schemes, while the poor orthogonality cannot be recovered without adopting other governing equations which include terms for non-orthogonality. Here we assume a situation that the orthogonality is the ultimate target, the aspect ratio is the second priority.
Orthogonal grid generation has often been conducted by solving an elliptic partial differential equation (PDE) like Laplace equation or Poisson equation [1][2][3] On the other hand algebraic methods may have advantage to generate orthogonal grids due to its inherent mathematical attribute. An algebraic grid generation method by Chibisov et al [4] uses circlerolling method with invariant radiuses to produce almost orthogonal grid, but their method is limited to situations when orthogonal grid is needed around an obstacle or a boundary, and cannot choose the other outside or inside boundary limit. In that case, their method is not ideal for resolving coastal areas which are clearly defined by surrounding land and offshore boundaries. Conti et al [5] proposed a modified PDE method by including an interpolation, but their results still contain poor orthogonality.
We assume that the distribution of coordinate lines has some degree of freedom as Duraiswami and Prosperetti [6] argued. We divide the lines into longitudinal lines and latitudinal lines for convenience. If we can arrange non-intersecting latitudinal lines at first, then the following longitudinal lines will not intersect at all. We attempt to generate grid points for a quadrilateral domain, two confronting boundaries of which are straight, and each of the other two boundary lines passes through a known position.  Even a circular domain like coast surrounding an island can be described by four quadrants. We start from choosing two straight confronting boundary lines AB and CD, then, we choose another straight line EF between the two boundary straight lines, which provides mid points of the latitudinal lines. We call AB, CD and interim quasi-parallel lines "longitudinal lines", and all other orthogonal lines from AB to CD lines "latitudinal lines". We can choose arbitrary ratio for selection of points along the three straight line segments, AB, CD, and EF. Each latitudinal line has 5 constraints: point positions and gradients at two side boundaries, and a point position between them. Therefore, a fourth order polynomial function for each latitudinal line is adequate to contain these five constraints.  Because the lower and upper latitudinal lines are not strictly straight, the above correction is not complete at the first trial, and requires repetition until a given allowance is met. This geometric iteration could also be replaced by other algebraic iteration method like Newton-Rhapson method.

Applications
The present algorithm is applied to two simple domains, a square-like one, and a partial quadrant (donut). Domain A is surrounded by three straight lines, and a polynomial curve, see figure 5.  (0, 1)).
An existing PDE method [3] is compared with the present method for Domain A. Four properties are compared for evaluation of the two methods, i.e. the maximum deviation from orthogonality ( t ), the mean deviation from orthogonality ( tt ), the maximum grid aspect ratio ( t ), and the mean grid aspect ratio ( tt ).
The maximum deviation from orthogonality is computed by: The mean deviation from orthogonality is computed by: The maximum grid aspect ratio is computed by: And mean grid aspect ratio is computed by: th t arc cos 㜷 th Computed four properties of the two methods are shown in table 1.  Results demonstrate that the present circle-inserting method produces almost orthogonal grid. Conceptually the present method should produce perfectly orthogonal grid, but small deviation from perfect orthogonality appears due to the finite sizes of the grid cells. On the other hand both methods show similar order of aspect ratios, although the present method produces slightly larger aspect ratios for this specific case. Next the present method was applied to a circular domain, Domain B, see figure  6. The present method reproduces a quadrant, and the maximum error in radius is 0.354%, which is thought sufficiently small. This test result demonstrates a full circular domain can be expressed by four partial quadrants, each of which can be represented by the grid produced by the present method.
The present method is relatively free from choosing the internals latitudinal lines, see figure 7(a). If we want a zoomed zone in the vertical direction, an option of vertical zooming can be used, see figure  7(b). Similarly a horizontal zooming zone can be assigned, see figure 7(c). Zooming on both directions