python - 使用 FG 方法将圆盘映射为正方形

Question

上下文:

我正在研究计算流体动力学中的网格。我想围绕一个圆生成一个结构化网格。我的计划是生成一个极坐标网格(下图左网格)，然后使用 FG 公式得到最终网格(下图右网格)。我正在使用 this article 第 4 页中的 FG 方法将圆盘映射到正方形。不幸的是，这篇文章没有提到如何处理公式中的奇点。表达式如下:

x = sgn(uv)/(v*sqrt 2)*sqrt(u**2+v**2-sqrt((u**2+v**2)(u**2+v**2-4u**2v**2)))

y = sgn(uv)/(u*sqrt 2)*sqrt(u**2+v**2-sqrt((u**2+v**2)(u**2+v**2-4u**2v**2)))

在编程之前，我正在努力解决这些公式的一些问题

问题

• 为什么这些公式映射以下点:(1,0)、(0,1)、(-1,0)、(-1,-1) 和 (0,0) 到点 (0,0) ？

• 我应该如何获得圆形和正方形之间的中间形状，如下图所示。

• 是否可以提供一种算法来从左侧网格中获取右侧网格？

这是我的尝试:

"""Map a circular computational domain with structured mesh around a circle (circular cylinder in 3D) to
Rectangular domain"""

import numpy as np
from numpy import sqrt, sign, pi, cos, sin
import matplotlib.pyplot as plt

def FGsquircle(u, v):
    SMALL = 1e-15
    t0 = u**2+v**2
    t1 = (u**2+v**2)*(u**2+v**2-4*u**2*v**2)
    t2 = u**2+v**2
    t3 = (u**2+v**2)*(u**2+v**2-4*u**2*v**2)

    x = sign(u*v)/(v*sqrt(2.0)+SMALL)*sqrt(t0-sqrt(t1))
    y = sign(u*v)/(u*sqrt(2.0)+SMALL)*sqrt(t2-sqrt(t3))
    return x, y

R0 = 1.0 # radius of the disc
RMAX = 5.0 # the radius of the outer circle in the domain
NT = 360 # num of division in the theta direction
NR = 10 # num of radial divisions
r = [R0+(RMAX-R0)/NR*k for k in range(NR)] # the radii of circles
theta = np.array([2*pi/NT*k for k in range(NT+1)])

u = [r[k]*cos(theta) for k in range(NR)]
v = [r[k]*sin(theta) for k in range(NR)]

u = np.array(u)
v = np.array(v)

x, y = FGsquircle(u, v)

我收到以下错误:

utils.py:21: RuntimeWarning: invalid value encountered in sqrt
  x = sign(u*v)/(v*sqrt(2.0)+SMALL)*sqrt(t0-sqrt(t1))
utils.py:22: RuntimeWarning: invalid value encountered in sqrt
  y = sign(u*v)/(u*sqrt(2.0)+SMALL)*sqrt(t2-sqrt(t3))

感谢任何帮助。

Answer 1

Why do these formulae map the following points: (1,0), (0,1), (-1,0), (-1,-1), and (0,0) to the point (0,0) ?

在这两个表达式中，您都除以 u 或 v，因此当两者之一为 0 时，表达式将变为未定义(不是零)。

How I am supposed to get the intermediate shape between a circle and a square as shown in the figure below.

只需绘制一个半径小于 1 的圆即可。

Is it possible to provide an algorithm to get the right mesh from the left one?

您只需要变换网格点。单元格可以保持不变。

<小时/>

示例代码:

from numpy import sqrt, sign
import numpy
import matplotlib.pyplot as plt


def f(x):
    u, v = x
    alpha = sqrt(
        u ** 2
        + v ** 2
        - sqrt((u ** 2 + v ** 2) * (u ** 2 + v ** 2 - 4 * u ** 2 * v ** 2))
    )
    return numpy.array(
        [sign(u * v) / (v * sqrt(2)) * alpha, sign(u * v) / (u * sqrt(2)) * alpha]
    )


for r in numpy.linspace(0.1, 1.0, 10):
    theta = numpy.linspace(0.0, 2 * numpy.pi, 1000, endpoint=True)
    uv = r * numpy.array([numpy.cos(theta), numpy.sin(theta)])
    xy = f(uv)

    plt.plot(xy[0], xy[1], "-")

plt.gca().set_aspect("equal")
plt.show()