Conformal Mappings with SymPy: Towards Python-driven Analytical Modeling in Physics

This contribution shows how the symbolic computing Python library SymPy can be used to improve flow force modeling due to a Couette-type flow, i.e. a flow of viscous fluid in the region between two bodies, where one body is in tangential motion relative to the other. This motion imposes shear stresses on the fluid and leads to a corresponding fluid flow. The flow forces exerted on the moving component are of interest in many applications, for example in system simulations of electrohydraulic valves. There, an eccentrically mounted cylindrical core (the armature) moves within an oil-filled tube (the polecap), experiencing fluid forces due to the viscous oil. SymPy can help to understand the range of validity as well as the limitations of analytical relations that are commonly used as standard approximations for these type of forces in many leading system simulation tools. In order to motivate these approaches, this contribution elucidates how the velocity of the flow is determined analytically by solving the Stokes equation in an eccentric annulus with a conformal mappingapproach. Afterwards analytical postprocessing leads to the corresponding flow force. The results obtained with SymPy are then checked against full 3D computational fluid dynamics (CFD) simulations. This work concludes with the combination of new Couette flow force approximations and similar results for the known Poiseuille flow (i.e. fluid flow induced by a pressure difference) to derive new relations for a combined Couette-Poiseuille flow force. This article is addressed to natural scientists and engineers that are interested in the application of conformal mappings and Taylor-expansions with the help of SymPy when solving partial differential equations analytically.


Introduction
In times of digitization and wide use of numerical methods in physics, the question emerges whether analytical tools, such as Taylor-expansions or conformal mappings, are still of interest and how they can be utilized in industrial and academic research.
Computational power has increased significantly in the last years and many physical problems, ranging from electromagnetism to fluid dynamics and structural mechanics, can be solved directly by numerically integrating the corresponding threedimensional PDEs, i.e. the Maxwell, Navier-Stokes or elasticity equations. However, when modeling physical systems such as hydraulic valves, transmission systems, engines, cars or planes, a direct 3D-approach for all relevant physical effects is still too difficult. In these situations, look-up tables containing a limited set of 3D-results are often included into 1D-system models for later interpolation. Alternatively, analytical approximations are used which are already included in the corresponding system simulation tools (e.g. Simcenter Amesim, [K19], [LGK21]). Figure 1 schematically shows such a valve system with input data from look-up tables and 1D-component symbols. Hence, in modern system modeling there are currently two main applications of analytical approximations: • Analytical approximations are included in the system simulation software components themselves, or • The user includes look-up tables for interpolation, entirely or partially generated with analytical approximations.
In this work we will focus on analytical approximations of flow forces that act upon the inner cylinder in an eccentric annular flow domain. Such forces are of interest in, for example, hydraulic valves that are electromagnetically actuated; see Figure 2. When the armature moves within the oil-filled interior of the polecap, that movement causes a Couette-type annular flow, i.e. a viscous flow due to motion of a solid body, between both components.
For an analytical treatment, this geometry has to be simplified considerably. Both armature and polecap are therefore modeled as solid and hollow cylinders respectively. Since in realistic scenarios, perfect concentricity between these two parts is rarely obtained, the armature can be supported eccentrically within the poletube. A cross-sectional cut perpendicular to the symmetry axes of both cylinders then leads to Figure 3. It shows the general case where an inner cylinder of radius R 1 is vertically displaced by a distance b from the center of an outer cylinder with radius R 2 . The eccentric annular flow domain is contained in the region between these two cylinders.
In leading system simulation tools, the flow force that acts upon the inner cylinder in Figure 3 is typically approximated by the relation Here µ denotes the viscosity of the fluid, l the common length of both cylinders, u R the velocity and δ = R 2 − R 1 the annular gap, i.e. the difference between outer and inner radius. Utilizing the capabilities of the open-source Computer Algebra System SymPy (as done e.g. in [MSP17]), we answer the following two questions: 1) How is Equation (1) related to the corresponding Stokes equation? 2) Does eccentricity ε = b/δ change this dependency and, if so, how exactly?
Furthermore, the velocities and forces obtained by solving the Stokes problem (i.e. the linear part of the Navier-Stokes system) with SymPy are compared to corresponding numerical solutions of the full, nonlinear Navier-Stokes equations, obtained from the commercially available Finite Volume tool ANSYS-CFX. Finally this article concludes with a note on the eccentric annular Poiseuille flow (that is a flow due to a pressure difference) and finishes with a comment on combined Couette-Poiseuille flow velocities and forces.

Material and methods
In order to solve the Stokes problem  [PHW33] and [BC09]) and Taylor-expansions, following [LGK21]. Equations (2) describe Couette flow when d p = 0 and u R = 0 and Poiseuille flow, when d p = 0 and u R = 0. Furthermore, Equations (2) describe Couette-Poiseuille flow when d p = 0 and u R = 0.

Solution of the Stokes problem within a concentric annulus for Couette-type flow
The solution of the Stokes problem within a concentric annulus for a Couette-type flow is well known, e.g. [LL87], and given by where r = x 2 + y 2 . This can easily be checked by using the diff function of SymPy. Keep in mind, that the natural logarithm in Equation (3) is denoted by log there.
import sympy as sym u_R, R1, R2, x, y = sym.symbols('u_r, R1, R2, x, y', real=True) u = u_R * sym.log(sym.sqrt(x ** 2 + y ** 2)/R2) / sym.log(R1/R2) laplacian = sym.diff(u, x, 2) + sym.diff(u, y, 2) It then follows that >>> sym.simplify(laplacian) 0 as expected. Further analytical solutions to the Laplace problem for other simple domains such as circles or rectangles can be found in e.g. [G13], [CB81] or [PP12]. In the following two Sections we will show with SymPy how the Couette flow problem within an eccentric annular domain can be transformed into a problem within a concentric annular region or within a rectangle. In these simple geometries analytical solutions to this problem are well-known. In order to transform the domains we make use of complex analysis, inspired by the French mathematician Jacques Hadamard (1865-1963): The shortest path between two truths in the real domain passes through the complex domain.
The ideas and strategies of conformal mappings using SymPy are mostly described in the following Section, where a Möbius transform is used.

Transformation to a concentric annulus with Möbius transforms
Using a Möbius transform (also called a bilinear transformation) in the form of an eccentric annulus in the complex z-plane can be mapped onto a concentric annulus in the corresponding w-plane. The Möbius transform used here is a slightly adapted version of the one presented in [BC09]; a is a constant (given in [BC09]) and will be defined further down in this Section. First of all, we will need some additional symbols for working with complex numbers and for the constant a. z, a = sym.symbols('z, a', real=True) Scaling the geometry in such a way that the outer circle ends up having a radius of 1 w = (z + sym.I * a)/(a * z + sym.I) w = w.subs(z, x/R2 + sym.I * y/R2) and separating real and imaginary part with SymPy functions xi_ = sym.simplify(re(w)) eta_ = sym.simplify(im(w)) one arrives at The latex rendering in the Jupyter Notebook shows directly the result of code in proper mathematical symbols, for instance >>> sym.simplify(im(w)) x (−R 2 − ay + a (R 2 a + y)) a 2 x 2 + (R 2 + ay) 2 After the scaling, the Möbius transform constant a reads as with c M given by Applying the Möbius transformation (4) to the boundaries leads to a concentric annular flow domain in the w-plane with inner radius 1 and outer radius R, given by (9) This new flow domain is depicted in Figure 4. Conformal mappings preserve harmonic functions, so the Stokes equation in the w-plane is of the same form as in the zplane. However, Equation (4) interchanges inner and outer boundaries. This will affect the corresponding boundary conditions one needs to specify there so that the Stokes-problem in the w-plane is given by Using the structure of Equation (3), the velocity in the w-plane is given by where ρ = ξ 2 + η 2 .
With the parameters specified in Table 1, the velocity in the w-plane (i.e. Equation (11)) can be used as an example for visualization and further evaluation.
The very convenient SymPy function lambdify is used to compute numerical values that are postprocessed by Matplotlib and depicted in Figure 5. The term R_ in the following code block denotes the numerical expression of the outer radius in the w-plane (see Equation (9)).

Mapping rectangles onto eccentric annuli by bipolar coordinate transformations
Another way of solving this problem utilizes conformal mappings related to bipolar coordinates. These coordinates are described in [PHW33] and are commonly used in elasticity theory (e.g. [L13] and [TG10]). For this contribution, we slighty adapted this  (2) transformation in such a way that it can be applied to the eccentric annulus of Figure 3. The mapping is given by where γ, c are constants from [PHW33] which are explicitly given in [W06] and [SL78]; the term i γ is added by the authors. Using this transformation, a properly chosen rectangular domain gets mapped onto an eccentric annulus; see Figure 7 for the domain in the w-plane. The boundaries are color-coded in order to visualize how the mapped borders are traversed in the z-plane. In addition the vertices are labelled and some coordinate lines are highlighted as well. Furthermore the left and right vertical boundaries in the wplane are identified in the z-plane, so periodic boundary conditions need to be applied to any PDE one wants to solve on the simple rectangle. Please note that for demonstrational purposes the radius of the inner circle in Figure 8 is reduced in order to indicate how the coordinate lines are distorted. For conformal mappings however, although distances between corresponding points and lengths of curves are changing, the intersecting angle between any two curves is preserved.
Further details on the relation between conformal mappings and bipolar coordinates can be found in e.g. [CTL09]. Inverting Equation (12) and separating real and imaginary parts as in the previous Section one gets Here, arctan 2 (y, x) is the 2-argument arctangent which returns the polar angle of a point with Cartesian coordinates (x, y). The constants from [W06] and [SL78] read as In the w-plane the corresponding Stokes-problem within the rectangular domain of Figure 7 is then prescribed by The last two equations specify the periodic boundary conditions one has to supply additionally. The solution to the system of equations (20) is easily obtained and given by the simple relation (21) Figure 9 shows a Matplotlib-visualization of the velocity distribution in the w-plane which is constant along ξ and increases linearly with η.

m/s is applied to the upper boundary
By again expressing η in terms of x and y, one obtains the very same velocity distribution in the eccentric annulus (in the z-plane) as already depicted in Figure 6.
It is interesting to remark, that Equations (11) and (21) look somehow related to each other due to the logarithm in both relations. However it is not immediately evident that they are actually identical. Nevertheless, due to existence and uniqueness theorems for the Stokes equation from [L14], one knows that relations (11) and (21) are in fact the same. Figure 10 compares these two analytically obtained velocities with results from a 3D computational fluid dynamics simulation (using ANSYS CFX) solving the full Navier-Stokes system. For these computations a velocity of u R = −0.4 m/s is prescribed onto the inner cylinder as boundary condition. All obtained velocities are evaluated along the symmetry axis of the annulus across the larger gap. The inner boundary is then reached on the left side, the outer boundary is hit on the right side of this Figure. As one can see, the two analytical approaches lead to the same velocity distribution across the larger gap and both boundary conditions are met exactly. On the other hand, due to the finite mesh size particularly at the outer radius R 2 , the boundary condition there is only approximately satisfied.
In the next Section, the corresponding flow force is obtained with SymPy-driven postprocessing and then compared again to the forces obtained by 3D-CFD and numerical evaluation.

Force calculation and comparison with 3D-CFD
The relation for the annular flow force that acts upon the armature in Figure This equation can be implemented in SymPy using the velocity distribution from Equation (11).
>>> u_w = u_R * sym.log(rho)/sym.log(R) >>> u_w u R ln(ρ) ln(R) Using the diff, subs and integrate functions from SymPy then leads to >>> Fe = mu * sym.diff(u_w, rho) >>> Fe = (rho * Fe).subs(rho, R) >>> Fe = sym.integrate(Fe, (z, 0, l)) >>> Fe = -sym.integrate(Fe, (phi, 0, 2 * pi)) >>> Fe −2π lµu R ln(R) Substituting the relation for R into F e , the flow force of the eccentric annular Couette flow is obtained. It can be manually adapated to the esthetic preferences of the authors, e.g. (23) Equation (23) therefore answers the second question posed in the Introduction: The flow force is decisively influenced by the eccentricity. Alternatively, the Couette flow force can be derived from Equation (21), which is obtained from solving the equivalent Stokes-problem in bipolar coordinates and for this case it is given by  With the data in Table 1 and Table 2, Figure 11 shows a comparison between the analytically obtained relations (23) and (24)   Again, both analytical relations agree perfectly but since the numerical CFD-results for the velocity slightly diverge from the analytical solution especially towards the outer boundary (as seen in Figure 10), the flow force computed from this data also shows smaller deviations.
Taylor-expansions and small gaps Equation (23) is even defined for the concentric case. Substituting b = 0 into this relation and simplifying the resulting expression leads to In order to finally answer the first question of the Introduction, i.e. how Equation (1) The answer to the aforementioned question then is: (1) is the leading term of a Taylor-expansion of the concentric annular Couette flow force around δ = 0. The contribution of this article closes with some additional remarks on eccentric annular Poiseuille flow and new possibilities of combining the results of the last Sections with results from [PHW33] and [LGK21].

Eccentric annular Poiseuille flow velocity
In various circumstances Couette flow may also induce a secondary flow driven by a pressure difference; a so-called Poiseuille flow. This particular type is of interest in many areas and we'll briefly show how the corresponding solution presented in [PHW33] is derived conceptually as well as how it can be implemented with the help of SymPy.
As far as we know, most of the current literature either refers to the aforementioned paper only by using its derived results (e.g. the volume flow relation found in [W06]) or by solving the Poiseuille problem numerically (as done in [TKM19]). The fact, that in the current context blood coagulation and hemodynamics are omnipresent in the media, eccentric annular blood flow in arteries is extensively studied ([TKM19]) and flow forces that act upon the arteries are of great medical interest (e.g. [S11]), makes it even more interesting to retrace the existing formulae of [PHW33], which are tedious to use when implemented by hand.
In the case of Poiseuille flow, the righthand-side of the corresponding Stokes equation is non-homogeneous (d p = 0; u R = 0); see also Equation (2). Hence, we need to deal with a different mathematical problem here compared to the previous Sections.
However, it possible to reduce the Poiseuille problem to an equivalent Couette problem with prescribed velocities on the boundaries (e.g. [M96]). That is the idea followed by [PHW33], who seek a solution of the form Here, Ψ is a harmonic function in the w-plane found by solving Laplace's equation in ξ and η. By using the conformal mapping of Equation (12) an appropriately chosen rectangle in the w-plane gets mapped onto an eccentric annulus in the z-plane, thereby preserving the harmonicity of Ψ. It then follows that ∆u = d p/(µl) in the z-plane and the boundary conditions for Ψ result from the task of eliminating the auxiliary term − d p 4µl (x 2 + y 2 ) on the boundaries associated with inner and outer radius.
The method described here is not only restricted to fluid dynamics. In elasticity theory, which inspired the work of [PHW33], Ψ is the harmonic conjugate of the so-called warping-or St. Venant torsion-function φ (see [L13] or [M77]), specified by The warping function helps to describe the elongation of an elastic cylinder that is also twisted. A practical implementation of φ can be found in e.g. [B14] and [BPO16] where it is called n inner 1,4 and where analytical approximations are compared to results from 3Dsimulations obtained with COMSOL. The velocity for eccentric Couette-Poiseuille flow can easily be found by superposing Equation (29) with one of the two Couette flow velocities derived in this contribution by utilizing SymPy.
The following relation shows such a superposed Couette-Poiseuille flow velocity, where both velocities where obtained by using the bipolar coordinate transformation (12) that maps rectangles onto eccentric annuli. Combining Equation (24) with the flow force from [PHW33], the overall exact analytical eccentric annular Couette-Poiseuille flow force that acts upon the inner cylinder is given by where Since the conformal mapping (12) is not defined for the concentric case b = 0, this drawback also translates to the corresponding forces in Equations (32) and (33). The relation above therefore is only defined for eccentric cases. However, the Couette flow force obtained with the Möbius transform, i.e. Equation (23), is defined for the concentric case as well. But since, to our knowledge, no one has ever constructed the Poiseuille flow velocity using a Möbius transform, the equivalent flow force (most likely defined for b = 0 too) is not available.
Therefore, the best analytical approximation for the eccentric Couette-Poiseuille flow force, defined both for the eccentric and concentric case, that we can present here, is a combination of Equation (23) and a Taylor-expansion of Equation (33) in the relative eccentricity ε = b/(R 2 − R 1 ) around ε = 0.
Here, F c is the well known Poiseuille flow force that acts upon the inner cylinder in the concentric case (e.g. [BSL07]) and a(κ) is a function of the ratio κ = R 1 /R 2 given by The particular approximation for the eccentric flow force due to a pressure gradient, i.e. F Piercy ≈ F c 1 + a(κ) ε 2 , was obtained for the first time in [LGK21].
To conclude this Section it is remarked, that again the useful SymPy function series can help in figuring out how a(κ) is approximated in the relevant practical case where R 1 ≈ R 2 .

Conclusion
This article showed that classical tools from mathematical physics, such as conformal mappings and Taylor-expansions, are still relevant and indispensable in times of digitization and wide use of numerics.
As an example, SymPy was used as a tool for symbolic mathematics in order to demonstrate that a popular approximation of the eccentric annular Couette flow force in modern system simulation tools is actually the leading-order term in a Taylorexpansion of the corresponding concentric annular force.
This force is calculated as special case of the more general eccentric annular Couette flow by postprocessing the resulting velocity distribution. Here, the velocity profile is analytically obtained by solving the equivalent Stokes problem with the help of conformal mappings, i.e. holomorphic functions in the complex plane.
The utilization of analytical methods is not solely restricted to fluid dynamics. Another application of SymPy in the context of PDEs in general could be homogenization. There, asymptotic expansions are substituted into the PDE and limiting problems are obtained in an algorithmical way, so SymPy might prove to be a valuable supporting tool. A starting point could be the introductory example from [BP89], which is worked out and compared to a FEM-solution obtained by COMSOL in [B14]. Furthermore, due to similar equations in axisymmetric electromagnetic problems (e.g. [LL84]), corresponding usage of conformal mappings and Taylor-expansions with SymPy is certainly possible there.
The authors think, that these methods may not only be applicable to mathematical physics but could be helpful in other areas as well, e.g. for understanding neural networks. Already available work described in [H10] and [H12] points in that direction and SymPy might be of great help in such areas, too.