Resolving Non-Symmetry in Flows via Subdomain Shifts ∗

. In this study, non-symmetric ﬂow problems are modeled by selecting subdomains and shifting them in such a way that the symmetry is recovered. As a result, the domains are made of simple grid structures and re-generation of mesh is avoided. Three test problems with various decomposition characteristics, namely, translation, rotation and deformation are selected, and they are analyzed in diﬀerent ﬂow regimes. To study the internal ﬂow between eccentric cylinders, two cylindrical concentric subdomains are considered, one translated relative to the other. Hence, a simple polar-coordinates mesh can be utilized instead of generating a mesh for the solution domain between the eccentric cylinders of the original problem. External ﬂow around a curvature tube is studied shifting the subdomain around the object in rotation, relative to the outer domain thus avoiding a re-generation of the mesh as the angle-of-attack changes. A third example involves deformation of an object exposed to natural convection, and the shifting of the domain facilitates the iteration process as the object deﬂects. Systems of nonlinear equations are solved within Newton-Krylov framework using the matrix-free approach. Geometrical and physical parameters are used to improve the solution process. Several results are provided to show the applicability of proposed method.


Introduction
Symmetry is a crucial criterion in investigating physical phenomena. Scientists exploit the symmetry of a problem as much as they can whether they are using analytical or computational tools. From a numerical point of view, the physics can be examined using the same number of unknowns but with higher accuracy, when the symmetric model is compared with the full model. Unfortunately, such a simplification cannot be attained in general. In fact, most of the problems that we are dealing with are non-symmetric, especially in fluid mechanics.
In this study, flow within or around geometric shapes is studied. In order to solve such problems, several approaches can be utilized. If the domain is simple in shape, like a rectangle or a circle, then cartesian and polar coordinates can be used with ease, respectively. If the geometry of the object is more complicated, then the domain can either be mapped into a computational domain in which the problem is again structured or one can use an unstructured mesh around the obstacle. If the domain is mapped, then the computations are involved with calculations regarding the Jacobian of the transformations. If an unstructured grid is generated, then special care should be taken during the discretization of the governing equations. In both cases, the mesh generation process will take some additional time which increases the computational overhead. This is especially significant if more objects are studied or objects are in motion relative to each other. A good alternative to solve this problem is the Domain Decomposition Method (DDM).
Although DDM was first used by Schwarz [25] in the late 1800's to prove the existence of solution to elliptic equations in complex geometries, nothing significant was done on DDM until early 1980's besides the studies of Kron [16] and Przemieniecki [20]. In 1986, Bjorstad and Windlund [3] studied solutions of elliptic partial differential equation on decomposed domains. Smith et al. [26] presented a comprehensive study of DDM and provided several algorithms in terms of parallel multilevel methods. A comprehensive analysis of the application DDM to partial differential equations is given in Quarteroni and Valli [22]. Cai et al. [5,6] investigated use of Additive Schwarz techniques to accelerate the solution of incompressible flow problems.
Houzeaux and Codina [12] proposed a solution algorithm based on Dirichlet/Robin interface conditions for of chimera grids. Presented technique is tested for various advection problems on simple square grids. An attempt to analyze domain decomposition in terms of the shapes of the subdomains is realized by Yang and Du [30]. They investigated the Alternating Schwarz method to model geometries with concave angle.
In fluid-structure problems, use of Arbitrary Lagrangian-Eulerian formulation can be applied with the use of dynamic mesh that is allowed to change with deforming solid bodies as investigated by Felippa et al. [11]. Later, Wall et al. [28] used overlapping domain decomposition to model fluid-structure problems. ALE formulation is used to investigate the flow field around a vertical beam where a subdomain in the flow geometry deforms with the elastic strip. A more recent study on overlapping domain decomposition is performed by Landmann and Montagnac [17]. They stated an algorithm to automate the selection of overlapping grids where subdomains can be in arbitrary shapes depending on the geometry of the model.
Even though domain decomposition has other applications like the decom-position of a mesh for parallelization purposes (particularly in finite element analysis) or accommodation of different physics (Navier-Stokes Equations in the vicinity of the boundary matched with Euler Equations for the far field), it is primarily used to divide complicated domains into many simple domains. Sub-domains may represent some portion of the boundary of the global domain or they may not, however they all possess the so called artificial boundaries which are composed of nodes along fictitious curves (or surfaces in 3D) used to carry values from neighboring domains to the domain in concern. Values are transferred either directly if the domains are conforming or interpolated if non-conforming. Conforming domains share common nodes whereas nonconforming domains do not have matching vertices. Procedures in DDM can be simplified for conforming domains using the so-called additive and multiplicative procedures [26]. In the additive approach, all sub-domains are solved simultaneously (useful in parallel processing) and the values are interchanged at the end of the computations. In the multiplicative version, on the other hand, values are interchanged as soon as they are calculated. In this case, convergence can be achieved with less number of iterations, however, parallel implementation is not trivial and coloring of domains might be necessary. In terms of stationary iterative methods, additive and multiplicative versions are analogous to Jacobi and Gauss-Seidel methods, respectively.
Overlap is also important in the Domain Decomposition Method. Nooverlap is the case in which the global domain is divided into distinct pieces. If overlap is utilized, the convergence is improved with the degree of overlap. On the other hand, it is shown in [26] that overlap is more effective on Alternating Schwarz algorithm rather then Additive and Multiplicative versions. A discussion on overlap is also presented in [3]. DDM can also be used to solve simple domains; in that respect, it can act as a preconditioner to linear problems as well as to nonlinear problems [5].
The methodology followed in this study does not fit into a multiplicative framework as in [26] -since no commons nodes are present -the updates between domains can still be thought as Jacobi and GS analogs. Different domains can be either be solved at the same time or domains can wait their turn and then proceed with recently updated boundaries. More on domain decomposition can also be found in different books [18,27,29].
Parallelization, in general, is an important issue in scientific computing and that reflects to recent publications [8,9] in mathematical analysis.
The structure of the paper is as follows. In Section 2, essentials of subdomain shifting is presented. Later, test problems will be introduced along with the formulation used to simulate the flow field. The paper will continue with the details of implementation of domain decomposition and numerical techniques. In Section 6, results for all three model problems are given. The discussion on the proposed idea will be followed by conclusion and further remarks for future work.

Subdomain Shifts
In Domain Decomposition Method, selection of the sub-domains is straightforward for a model with simple geometry (Figure 1a). If the shape of the domain is much more complicated, then the division can still be made with simple cuts, yet sub-problems will vary in size since the distribution of the unknowns is virtually ignored (see Figure 1b. In addition, data exchange between the subdomains should also be reduced (this condition is especially important in parallel processing). In that respect, a more elegant solution to determine the decomposition might be to use spectral partitioning strategies like Metis [15] which is based on multilevel techniques to generate sub-domains for irregular geometries. A downside of this method is that, the procedure should be repeated if the mesh is modified, i.e. if the domain is changed. When the problem in concern is to be analyzed in different parameter spaces that depends on the geometry like the eccentric cylinders problem in Section 3.1, a direct correlation cannot be done with the subdomains. In this study, an idea is proposed to select the subdomains in such a way that separate problem resembles with each other and has a physical significance. The method is applied on non-symmetric problems which means that the geometry cannot be reduced into a smaller model using symmetry boundary conditions. In Figure 2 the geometry of a eccentric cylinder problem is given (Section 3.1). If the setup were concentric, then the centers of the outer and inner cylinder would coincide. However, the center is shifted in an arbitrary direction (which occurs in journal bearing when operated under load) so a computational grid should be generated to comply with the spacing between the cylinders. Instead, the problem can be solved with overlapping domain decomposition in such a way that the subdomains are formed like concentric cylinder using artificial boundary conditions. Domain I, Ω I , excludes the inner cylinder (denoted by dots) and includes a new boundary denoted by Γ I whose values are interpolated from Domain II, Ω II . Similarly, Ω II excludes the outer domain and works with artificial boundary Γ II . This time, values are interpolated from Ω I .
If the eccentricity is changing, i.e. inner cylinder is moving, then only the points that are used for interpolation will be changed. Hence, symmetry is preserved even if the domains are in motion. The geometric definitions of the subdomains will remain the same. Since the domains are now concentric, simple polar grids would be used with ease. If additional division are necessary, then decompositions can be performed on those two subdomains separately. Also, a recursive approach can be performed for example, if four domains to be used, than four concentric cylinders can be generated.
In this study, three different configurations are presented to demonstrate subdomain shifts. Although any restriction that applies on domain decomposition method also applies to this idea, main concern is geometric complications. Depending on the selected parameter, the geometry itself can be a constrained not to employ proposed method. In eccentric cylinder problems, the epsilon has a limit -which is also complication for domain decomposition technique, as well.

Mathematical Modeling
In this study, three model problems are selected to test the performance of the subdomain shifting method (SSM) in three different scenarios. In the first case, translation of subdomains are to be examined. For this purpose, an internal flow problem is selected where two cylinders are placed eccentrically (see Figure 3). The eccentricity, ε, is used as a geometric continuation parameter. In the second problem, rotation of subdomains is investigated. This time an external flow is studied at which flow around a curvature tube is to be found ( Figure 4). Different values of the angle of attack, α, are investigated by just rotating the inner domain. In the last situation, deformation of boundaries is examined. In this free convection problem, a strip is to be bent in the presence of a temperature difference ( Figure 5). The radius of curvature, ρ K , acts as the continuation parameter which defines the deflection of the beam. These problems are selected is such a way that they are mathematically simple yet related to engineering applications. Test problems are incompressible flow problems. The main difficulty in solving incompressible Navier Stokes equation is treatment of the pressure. The continuity equation lacks pressure terms so it behaves like a constraint to be fulfilled by the velocity field. To date, numerous methodologies are proposed to analyze incompressible flows, e.g. SIMPLE [19] or Artificial Compressibility [7].
In this study, however, stream function -vorticity formulation is selected which avoids the solution of pressure equation and reduces the numbers of unknowns to 2 (ψ, Ω) rather than 3 (u, v, p). Before giving the details of the models, it is useful to introduce the velocity and vorticity in both Cartesian and polar coordinate systems (see, Table 1). Stream functions are selected is such a way that the continuity is satisfied is respective coordinate systems. As seen from the table, velocities are defined as the derivative of stream function. Details of the models are given as follows.

Translation of subdomains -flow between eccentric cylinders
Analysis of flows between eccentric cylinders is significant in a wide range of problems in the field of engineering. The lubrication of bearings and shafts, viscosity measurements in torsional viscometer etc. are some of the basic examples. The flow is dependent on the velocities of the boundaries as well as the eccentricity of the geometry.
In the literature, most of the studies for eccentric cylinders are presented only for one moving cylinder [10]. Using the idea introduced in Section 2, rotation of both cylinders at the same time can be investigated. In this study computations are carried out for eccentricities up to 0.3 where the eccentricity, ε, denotes the amount of shift of the center of the inner cylinder off the outer cylinder. The analysis starts with concentric case where ε = 0. This is useful in different aspects. First, it has an analytic solution. Second, it will help us to provide a good initial guess while varying the eccentricity of problems. In this way, the problem is converted into a continuation problem over the parameter ε, in addition to the Re number which is the key physical non-dimensional parameter.
Geometry of eccentric cylinders can easily be modelled in polar coordinates. It can be seen from Figure 3 that the domain is divided into two overlapping parts. Domain 1, Ω 1 , is used to model the outer wall and excludes the inner cylinder. Its inner boundary is artificial and denoted by Γ 1 . On the other hand, domain 2, Ω 2 , includes the inner wall and exchanges data with Ω 1 over Γ 2 .
As shown in Figure 3, the subdomains are very similar in shape, in fact with the use of SSM, the whole domain is reduced into two concentric problems with artificial boundaries. Consequently, the number of nodes on both domains can be arranged to be the same, which simplifies load balancing of parallel processing where domains can be solved on different processors and exchange data at the end of each iteration (Jacobi Analogy). In this approach, load balancing is not a big issue since the dimensions of the nonlinear systems will be the same.
Throughout this study, incompressible flow problems are considered. Main issue in the analysis of incompressible flow problems is the treatment of the continuity equation which is related to the computation of the pressure field. Stream function (ψ) -vorticity (Ω) approach is selected to analyze the model problems although several other methodologies are available in the literature like SIMPLE [19] or Artificial Compressibility [7]. This approach reduces the number of unknowns and eliminates numerical complications arising because of the incompressibility constraint that is related to the lack of a pressure term in the continuity equation.
The non-dimensional form of the governing equations in polar coordinates is given by where the Reynolds number, Re, is defined as Re = ρu θ D/µ, where ρ is the density, u θ is the tangential velocity and D is the diameter of the rotating cylinder. On the wall of outer domain, (Ω 1 ), stream function is set to be zero whereas the vorticity is calculated as Different boundary conditions can also be found in [23].
The same equation is also used for inner domain, (Ω 1 ), however the stream function is evaluated in such a way, that the tangential velocity u θ , related to the rotation of the cylinder, is equal to 1: Artificial boundaries, Γ 1 and Γ 2 are treated as non-homogeneous Dirichlet data, and calculated with bilinear interpolation as explained in Section 4.

Rotation of subdomains -flow around a curvature tube
Our next model is an external flow problem at which the flow around a curved geometry is analyzed. This problem is a simple 2D model of a paraglider wing. Paragliders are composed of a wing called a canopy which deforms during flight since unlike a conventional aircraft wing, it does not have a rigid structure [1]. The canopy takes the form of a wing by the pressure of the air entering in through the openings at the leading edge. When the airfoil is moved at the allowed angle of attack range, lift and drag are created which counterbalance the gravitational force. Hence numerical analysis of paragliders is important in terms of performance related issues [1]. Additionally, the flexible nature of paragliders make them an excellent example of fluid-structure interaction problems.
Like the first problem, the whole domain is divided into two overlapping subdomains. Domain 1, Ω 1 , is used to model the far field conditions and keeps the curvature tube out while domain 2, Ω 2 , carries information to the first domain about the object. Ω 1 , has an artificial boundary, Γ 1 and studied in cartesian coordinates. On the other hand, Ω 2 , with interface Γ 2 is examined in polar coordinates. When the global domain is considered alone, then a mesh generation operation will be needed if SSM is not utilized. However, with the selection of those overlapping domains, the problem is split into two simple subproblems which have to be solved iteratively. The angle of attack (see, Figure 4) is set by just rotating Ω 2 .
In this model, the domain decomposition is used with two coordinate systems, polar and cartesian, so governing equations for both are needed. For the polar setup i.e. for Ω 2 , equations (3.1) and (3.2) are reconsidered. Cartesian equations for Ω 1 are given as The boundary conditions for polar domain (Ω 2 ) are the same as those used in model 1 (Section 3.1). In cartesian domain, inflow condition is specified properly with the stream function: Vorticity at the inlet and the far field is set to be zero, at the outflow equation

Deformation of subdomains -natural convection around a strip
In this last problem, deflection of a bimetallic strip is considered. It is composed of two metals, which differ in their coefficients of thermal expansion. If the strip gets heated, then two materials will elongate in some amounts, however their elongations will be different. Since they are also bonded together, the beam will deflect. The deflection of the solid will change the boundaries of the fluid domain. Normally, such a process needs a re-meshing if a single fluid domain is used, however, with the help of the domain decomposition regeneration of the grid is avoided. In fact, by using SSM the deformation is reflected only to the inner domain, Ω 2 , as observed in Figure 5.
In this model, the solid analysis is kept simple, i.e. the use of finite element method is avoided. The curvature, K, is calculated using equation (3.9), which is a formula introduced by Timoshenko. Hence, no real computations are performed to calculate the deflection of the strip and the curvature is found by analytical means. One can also consider this model as an introductory step to a more complex analysis in fluid-structure interaction where the deformation is calculated with thermoelasticity rather than analytic means. In this study, only an analytic methodology is followed in terms of the deformation mainly because we are interested in investigating subdomain shifting method on flow problems rather than the changes inside of the solid part. Since the deformation is not complicated, grid generation process is simplified. In fact, the shape of the grid is straight forward, i.e. it is also in polar coordinates since the bending can be determined by calculating the curvature of the final form in linear elastic theory. If elastic moduli of both metals are the same, E 1 = E 2 , and the parts have the same thickness then the curvature K can be found using equation Parameter ρ K given in equation (3.9) is the radius of curvature and this parameter is used to generate the inner domain which is proportional to the thickness of the beam, h, and inversely proportional to the difference of the temperature and the thermal coefficient of expansion. ∆α is the difference of thermal expansion coefficients (α 2 − α 1 ) and is used to adjust the deformed geometry for numerical tests. Deformed shape is used to perform a continuation strategy in such a way that it acts as an initial guess to another problem with a large temperature difference between the strip and the environment. So for different configurations, the analysis will not start from scratch, but from the last configuration.
The governing equations of natural convection are different from those of forced convection model problems introduced before. In natural convection, density changes are the main driving force of the flow in the presence of the body forces -which is the gravity in current study. Although the flow is incompressible, density changes can be analyzed using the Boussinesq approximation [13], where β is the coefficient of thermal expansion of fluid: In a natural convection problem, the energy equation should also be added to governing equations to solve for the temperature, T . The non-dimensional form of the governing equations are slightly different, and are given by equations (3.10)-(3.15): Here, g is the gravity, H is the height of the solution domain and α T is the thermal diffusivity. The stream function is zero on the boundary (ψ b = 0) all over the domain. Vorticity conditions for polar domain are as in equation (3.3). For the cartesian domain, vorticity at the walls are defined by In this equation, wall − 1 denotes the first interior node from the wall into the domain in the normal direction, n. ∆n is the spacing between the interior and boundary nodes. Now we have to introduce boundary conditions for the dimensionless temperature. In this problem, top and bottom walls are insulated whereas temperature on the left wall, T L is zero. At the initial configuration, when the strip is straight, temperature at the right wall and the strip, T R =1. As the T R increases, the strip bends and a new problem should be solved. It is assumed that an external heat source acts on the right of the domain and keeps the right wall and the strip at the same temperature during the analysis. When T R is varied, the desired curved geometry can be specified by changing with ∆α as a computational parameter.

Domain Decomposition
Domain decomposition is a division of a global solution domain to some number of sub-domains. In overlapping domain decomposition, sub-domains share common regions. On the contrary, in non-overlapping domain decomposition, each sub-domain represent a distinct and unique part of the global problem. Non-overlapping version is generally preferred in decomposing the domain for parallel processing especially in finite element analysis. In this study we have used non-conforming sub-domains.
Domain decomposition of overlapping type is also can be split into two types: the domains are either conforming or non-conforming. In conforming decompositions, some nodes belong to multiple subdomains. In non-conforming problems however, nodes are not necessarily on top of each other. Analysis on conforming domain decomposition can be simplified by additive and multiplicative versions as explained in [26]. In this study, one sub-domain is kept fixed while the other is allowed to be in motion either in terms of translation, rotation or deformation -so are the nodes that belong to those moving regions. Consequently a non-conforming domain decomposition is suitable to analyze the models hence nodes of different domains do not have to match.
In domain decomposition, data are exchanged through artificial boundaries. This exchange is mostly done with non-homogeneous Dirichlet boundary conditions. The Robin type boundary condition, as a combination of Dirichlet and Neumann type boundary conditions is also possible [26]. Data transfer is achieved through bilinear interpolations. Interpolation on cartesian coordinates is calculated using equation whereas that on the polar coordinates is handled with equation See also Figure 6). High order interpolation schemes are also possible, however, successful computations are performed with bilinear operators. a) b) Figure 6. Interpolation: a) on Cartesian, b) polar coordinate systems.

Numerical Methods
To simulate the models given in Section 3, the governing equations should be discretized and for that finite differences with second order accuracy is used. First order derivatives are also computed with central differences. Final form of the algebraic equations after the discretization are coupled and nonlinear, so proper solution approaches should be employed. The Newton-Krylov techniques are used to deal with the resulting nonlinear system of equations. Newton-Krylov methods allow the analyst to perform matrix-free computations. While the Newton method deals with the update of the solution vector, Krylov methods seek a solution of linearized equations by avoiding the formation of the Jacobian. Preconditioning is also utilized during the computations to enhance the convergence of the system of linear equations since the convergence of Krylov solver is limited without the implementation of a preconditioner. After discretization of governing equations and implementation of the boundary conditions the resulting system of equations can be written as Starting from a initial guess x 0 , the iterates can be found by the iterative process where λ is the damping parameter. In Newton's method the update ∆x is found using where J is the Jacobian Matrix such that J ij = ∂fi ∂xj . In solving equation (5.3) a Krylov method is preferred. In Krylov methods, like GMRES [24], only matrix-vector products (J v) are needed rather than the matrix J itself (except in preconditioning). This multiplication can be approximately evaluated using equation (5.1) through the following formulation Matrix-vector product is computed by perturbing x by an amount of ǫv. Accuracy of the differencing is dependent on the selection of the ǫ [4,21]. In this study, we take ǫ ≈ 1e − 7, which is accurate enough to proceed with the calculations. During the numerical tests, Jacobi preconditioning (namely diagonal scaling) is used primarily. This is the simplest preconditioner and the reason for its use is that the computations can still be treated in a matrix-free setting. For improved convergence properties, LU-SGS is used (in second test problem) which is also handled in a matrix-free fashion [14]. This is a variant of LU decomposition. Another variant is the incomplete LU (ILU). ILU(0) is utilized in the second test case, where the Jacobian is stored in compressed column storage to speed up the computations [2]. Here, the computations are not matrix-free anymore but can be treated as jacobian free because the preconditioner could be selected from a low order problem (first order accuracy rather than second order).

Results & Discussion
In this section, results for all model problems are presented separately in different aspects. In all three tests Newton's Method is utilized as in equation (5.2). GMRES is used to solve equation (5.3) for ∆x in a matrix-free context. Stopping criteria for nonlinear iterations was 1e − 6 in the absolute norm ( f 2 < 1e − 6) and for linear iterations it was 1e − 5 in the relative norm ( r n 2 / r 0 2 < 1e − 5). For the convergence of the domain decomposition the criteria given in [3] is used, i.e. subdomain visits are stopped if x k − x k−1 is reasonably small. In the third model problem a parametric study in terms of the radius of curvature, ρ, is performed and subdomain convergence limit is set to 1e − 8.

Translation of subdomains
In this internal flow problem, cases with varying Reynolds number, Re, and eccentricity, ε, are considered.  Figure 7 shows the contours of the stream function and u θ for Re = 100 and ε = 0.3. As shown in this figure, the streamlines are shifted with the geometry. For the case when outer cylinder is rotating with -0.5, the stream function plot is completely different as observed from Figure 8. Here, a cell is formed at the left of the inner cylinder. This solution is achieved using continuation on ε and the computations are started with concentric case.
In Figure 9, nonlinear residual history for both inner and outer cylinders are showed for Re = 100 and ε = 0.2 (inner cylinder rotating with u θ = 1 and outer cylinder rotating with u θ = −0.5). It can be seen that, the nonlinear norm at the beginning of each subdomain visit is reduced, generally. The graph is noisier for the inner cylinder especially for the first couple of visits.

Rotation of subdomains
This second problem is the first model where external flow is examined and two sets of coordinate systems are used. In this part, flow around of a curvature tube at Re = 160 is studied for two different angle of attack, α, values. Outcome of the analysis with α = 0 o is given in figure 10. On the left, stream function is given whereas vorticity is on the right. We should pay attention to the fact that although this is a zero angle of attack problem, the contours are not symmetric since the geometry is curved. A symmetric case might be generated with α = 90 o .
To fit into the context of non-symmetric flows, another case with α = 42 o is investigated instead of this symmetric case. Results are given in Figure 11 at which wakes are observed easily downstream of the curvature tube. Rotated sub-domain is also apparent in this second instance.
In this model, effect of different preconditioners on linear solves is also investigated. As seen in Figure 12, LU-SGS is competitive to ILU(0) in cylindrical domain however in all domains ILU is superior to all other preconditioners. Not much improvement is attained with Jacobi preconditioner on GMRES.
A note on LU-SGS should be said before proceeding, although LU-SGS is as successful as ILU(0), the calculations take a lot of time since matrix-free computations are done for 2N times for N number of unknowns. This is because both forward and backward substitutions are needed in incomplete decomposition preconditioning, respectively. On the other hand, if we compensate on matrix-free methodology, the evaluations can be accelerated with the use of compressed storage schemes. In this study, compressed row storage (crs) is applied and overhead of preconditioning is reduced. However, at this point it can be said that one is free to decide either on ILU or LU-SGS since both need storage for acceptable computational performance. To test the performance of different solver combination, Jacobi, LU-SGS and ILU(0) preconditioner are used with CGS, BiCGStab and GMRES(100). Table 2 gives the results for one typical visit of the interior domain (i.e. polar domain) for Re = 160, where NS denotes the number of Newton steps, TI denotes the total linear iterations and LI linear iterations per Newton step. Number of subspaces for GMRES is fixed at 100, that is why Jacobi-GMRES(100) fail to converge for a given amount of maximum linear steps (500) per Newton iteration. When the chart is examined, SGS and ILU(0) are superior to Jacobi preconditioning, where ILU(0) is slightly faster than LU-SGS. Among different combinations, ILU(0)-BiCGStab is the fastest solver. In this model problem, solvers fail to converge w/o a preconditioner.

Deformation of subdomains
In the last model problem, free convection for varying Rayleigh numbers, Ra, is analyzed while keeping the Prandtl number, P r, as unity. Different geometries are also investigated by changing the radius of curvature, ρ, of the bimetallic strip. In order to improve the convergence of the discretized system, continuation is used with both Ra and ρ. Computations are started with Ra = 10 4 which is a reasonable value to start without continuation. Stream function and temperature profiles for Ra = 10 5 and ρ = 6.5 are given in Figure 13. Similarly, another set of results for Ra = 10 6 is given in Figure 14. Continuation is used for both cases to improve the convergence.
A typical linear convergence history (see, Figure 15) as well as nonlinear residual history (see, Figure 16) are also shown. In the computations, total number of nonlinear correction steps increase as Ra increases which is expected since the nonlinearity of the governing equations depends on Ra number. Another observation is that total number of subdomain visits decreases as the Ra increases, however, Newton steps per visit shows an increasing trend. When radius of curvature, ρ, is in concern it can be said that more visits are needed as we reduce ρ. This is also what would be predicted before the computation since a large value of ρ mimic the fully cartesian setup and this is simpler to analyze compared to polar geometry. In addition, number of subdomain visits can be cut down if continuation is applied in terms of the radius of curvature as seen in Table 3, where numbers above denote the number of total subdomain visits and numbers below show the number of total number of nonlinear iterations. When Table 3 is examined, it can be said that continuation in Ra and ρ reduces the number of subdomain visits. In order to get the result for Ra = 10 6 and ρ = 6.5, two different paths can be followed. Either one can first fix ρ as ρ = 999.5 and start with Ra = 10 4 to reach Ra = 10 6 and then change the radius of curvature. The other path is also possible, first compute for ρ = 6.5 at Ra = 10 4 and then continue to Ra = 10 6 . It is seen from the table that the former approach requires 194 subdomain visits, whereas the latter path needs 292 visits. However, when the total Newton steps are considered, the picture is just the opposite: 9129 to 5934. The second option is also more efficient in terms of the average Newton steps: 47.1 to 20.1. As a result it can be concluded that we need two times less Newton steps per visit if continuation in ρ is favored.

Conclusion
In this study, a decomposition method on non-symmetric flow problems is investigated. The domains are kept simple as a result discretizations are done in standard way. Also, the burden of generating different meshes for each geometry is replaced by the iterations between subdomains. It is showed that, by relevant selection of subdomains and continuation parameters, subdomain shifting can be used as a solver.
Problems examined in this study are simple models of different engineering problems. First model was an eccentric cylinder problem where subdomain shifting is tested in such a way that translation of subdomains is investigated. The translation is achieved over the eccentricity, ε. Next, flow around a curvature tube is analyzed which represents a simple model of a paraglider wing. In this problem, two domains are solved in two different coordinate systems namely, cartesian and polar coordinates. This approach is tested for rotation of domains, so different values of angle of attack, α, are given by rotating the inner domain. The last problem is designed to analyze the deformation of subdomains. A bimetallic strip is bent with the presence of a temperature difference from the initial configuration. The deformation of the domain is reflected only onto the inner domain which is in polar coordinates, so grid generation all over the domain is avoided and the outer domain is kept in cartesian coordinates. Radius of curvature, ρ K , is the geometric parameter used to model the deformation.
In terms of fluid mechanics, three different flow problems are examined in three test cases. An internal flow problem is investigated in the first case, whereas in the second model, an external flow problem is examined. Reynolds, Re, number was the key physical parameter for both models. These problems can be considered as the first step of solving a forced convection problem, however, the solution of the energy equation is avoided for these two cases because the velocity field and the temperature field are decoupled. On the contrary, the third case was a free convection problem in which the solution of energy equation to calculate the temperature distribution was essential since the equations are coupled through the buoyancy term. For this model, various Rayleigh number, Ra, flows are analyzed.
Stream Function Vorticity approach with relevant boundary conditions is utilized to model the dynamics of the flow. During the discretization process of the governing equations, we took the advantage of simple domains and finite difference algorithms. Systems of nonlinear equations are solved taking Newton-Krylov approach and using the matrix-free methodology. In matrixfree approach, directional differencing is used to replace matrix-vector multiplications in Krylov solvers (e.g. GMRES), while avoiding the formation of the Jacobian in Newton's method. In order to improve the convergence of the linear solvers, different preconditioners are tested. Diagonal scaling is applied in all test problems within the matrix-free context which provided mild improvement in linear solver. In addition, ILU(0) and LU-SGS preconditioners are tested and it is observed that LU-SGS is comparable to ILU(0) in number of iterations and it can still be applied in matrix-free approach however, the computations are taking too long and compressed storage schemes should be used to reduce the time overhead.
Parallelism is an important feature sought in today's numerical computations. Since the subdomains are similar in shape, load balancing issues are simplified. Also, with the use of subdomain shifts, mesh generation for moving objects is avoided so parallelization can be still be applied as stationary objects. Further analysis can be carried out in terms of the parallel processing to take the advantage of the proposed idea.