Surface Pressure Calculation Method of Multi-Field Coupling Mechanism under the Action of Flow Field

Abstract

At present, the majority of fluid mechanisms are multi-field coupling mechanisms, and their function is also achieved in the flow field. Therefore, calculating the aerodynamic characteristics of the multi-field coupling mechanism in a symmetric flow field is very important. However, at present, the strong coupling algorithm in the domain has the problems of low accuracy and computational efficiency, which make it more difficult to solve the coupling problem. This article obtains the vector potential of the law of conservation of momentum using the tensor analysis method in a Cartesian coordinate system. Meanwhile, the generalized operator of Navier–Stokes equations and the fundamental solution of the generalized operator are obtained on this basis. Then, this article proposes the boundary integral equation of the Navier–Stokes equations by combining the fundamental solution of the Laplace equation with generalized potential theory. Based on this boundary integral equation, this article has developed a new calculation method that can help achieve integral calculation without domains, greatly reducing the problem’s difficulty. Finally, by comparing the ellipsoid example solution with the experimental results, the algorithm’s reliability in solving the incompressible problem is verified.

1. Introduction

Since the multi-field coupling mechanism is widely used in engineering, it is the hot spot and focus of current research. In the multi-field coupling mechanism, fluid–structure interaction has substantial research and application value, such as aviation, aerospace, navigation and other fields. Numerical solutions for fluid–structure interaction have come a long way in the past few decades and have become one of the hottest topics in research. The core of the coupling solution process is calculating unsteady flow problems with moving boundaries and moving meshes because the flow domain’s size and shape change with the structure’s movement or deformation. Still, there are also a series of problems worth discussing.

Generally speaking, fluid–structure interaction is divided into strong coupling and weak coupling algorithms according to the control equation. Current studies of weak coupling (Figure 1a) have focused on arbitrary Lagrange–Euler (ALE) [1] and immersion boundaries (IB) [2]. Since these two methods rely on Euler’s formula, additional equations are required to study free-surface flow, resulting in low computational efficiency. In addition, ALE and IB methods require remeshing when dealing with significant deformation of the solid phase [3]; to accurately represent the interface, the IB method requires continuous mesh adaptation at the two-phase interface, reducing the computational efficiency. In other technical ideas, tracking reference space or evolving deformation gradients into field variables can solve the fluid–structure interface coupling problem in the Euler framework [4] but dramatically increase the calculation. Overall, the above weakly coupled method has a positive effect in improving the computational efficiency [5,6], but the computational accuracy is reduced due to grid instability [7]. The strongly coupled algorithm (Figure 1b) is discretized by the grid method [8], or particle method [9], and the fluid domain and the solid domain are solved directly, which can ensure high computational accuracy [10]. However, the strongly coupled algorithm of the grid class in the domain usually does not scale linearly with the size of the unknowns because the amount of computation usually does not change, and, even considering the sparsity of the matrix, a large set of unknowns can negatively affect the solution speed [11]. Moreover, the two phases will also bring algorithm stability problems at different scales [12]. Therefore, the strong coupling algorithm in the domain has the dual problem of computational accuracy and computational efficiency (Figure 2).

Since the boundary element method (BEM) is only discrete on the boundary, it can solve the problem of algorithm stability caused by different scales and display phenomena and effectively reduce the size of the unknown [13]. Therefore, the strong coupling method based on boundary elements effectively improves the computational efficiency and accuracy with great research value (Figure 2).

For the boundary element method, Yang D. et al. [13] used the fundamental solution of the Stokes equation to study incompressible viscous fluids. Qiang H. and Yao Z. et al. [14] established a correlation matrix between interface force and displacement according to the constitutive relationship of the inclusion material. Zhu R. et al. proposed a vertical integration method based on the multi-domain higher-order boundary element method [15] and a time-domain calculation method for wave propulsion [16]. Zhu J. et al. [17] found that if we want to calculate the regional integral accurately, the boundary element method needs to consider the discreteness of the internal elements of the computational domain and lose the advantage of dimensionality reduction. To this end, Nardini and Brebbdia [18] proposed a double reciprocal method for solid mechanics which uses the basic function to convert regional integrals into boundary integrals by expressing force effect quantities as basis functions. Power and Florez [19] used DRM technology to convert hydrodynamic area integrals into boundary integrals. Gao W. et al. [20,21] proposed the radial integral method (RIM) to solve the problem involving different types of regional integrals in the same integral equation and developed the singular integral boundary element method.

There are still many problems in using boundary element algorithms to deal with incompressible fluids. Greengard et al. [22] believe it is difficult to solve complex nonlinear phenomena with boundary elements because obtaining the solution of Green’s function is complex. In the DNS algorithm, the technology of directly simulating the Navier–Stokes equations using the initial value condition has matured [23]. This shows the characteristics of the Navier–Stokes equations that can solve the initial boundary value problem, the Green function, as a generalized function used to solve the non-homogeneous differential equation of the initial boundary condition [24]. The existence of its solution is unquestionable [25]. The difficulty is that nonlinear phenomena are often related to the stability and sensitivity of the resolution, even with standard methods such as Lie algebra, the perturbation method or geometry. Due to the existence of different applicable conditions, it is not easy to describe the nonlinear phenomenon mathematically; thus, it cannot be solved.

At the same time, due to the complexity of fluid flow, there are still many problems that cannot be solved at present in incompressible fluids. They involve boundary layer separation, flow transition, vortex generation, evolution and interaction between the fluids [26]. These problems make it very difficult to obtain the aerodynamic characteristics of object symmetric surfaces. Meanwhile, the problem of flow past an object is nonlinear, which makes it more difficult to obtain accurate solutions. Yue P. et al. [27] have proved the mathematical principle of flow field separation and transition, which makes it possible to further accurately solve fluid mechanics problems and obtain the aerodynamic characteristics of an object’s surface. This lays the foundation for calculating the pressure distribution in a complex flow field.

Since the Navier–Stokes equations can solve the initial boundary value problem, there are still many difficulties in solving the Navier–Stokes equations [28,29,30,31]. At present, for the solution of fluid flow, only the simplest parallel flow can obtain an accurate solution. In the case of non-parallel flow, and even more complicated flow, it is difficult to solve their accurate solutions [32,33]. Among them, Rostamzadeh A. et al. [34] proposed a new MACB method for the simulation of incompressible viscous flows with heat transfer based on the Navier–Stokes and energy equations in non-dimension form. However, this method cannot be applied to the solution of boundary value problems. Therefore, new methods are needed to achieve accurate boundary calculation of Navier–Stokes equations.

In this regard, Lu Y. et al. [35] have derived basic differential operations and generalized Green’s integral theorems based on Cartesian tensor analysis. Their conclusion is applicable to multi-dimension space and solves some nonlinear problems more easily. This indicates that the new tensor analysis method can solve some complex nonlinear problems that cannot be solved at present. Meanwhile, according to the research of Vectors A. R. [36], it is of great significance to derive tensor analysis for solving aerodynamic characteristics. Based on the fundamental theorem, the generalized potential theory could be described in the form of a vector tensor. In other words, it is possible to propose the integral expression of Navier–Stokes equations (boundary integral equations) based on generalized potential theory and tensor basic solutions [37,38]; for example, a generalized Green formula. This would make it easier to solve a series of complex boundary value problems in fluid mechanics [39,40]. At the same time, the essential work of Erich R. et al. [41] reveals that functional theory plays a vital role in solving mathematical physics problems at the boundary integral equation. The introduction of generalized potential theory makes it possible to improve the Navier–Stokes equations and solve a series of complicated boundary value problems in fluid mechanics [42,43,44]. These works have laid a solid foundation for studying the mechanical properties of the multi-field coupling mechanism surface under the action of a complex flow field.

This article obtains a unique set of fundamental solutions through mathematical derivation that allows boundary integration of all terms in the Navier–Stokes equations. This domain-free integration method can eliminate the error of existing domain integration. It can directly calculate the surface pressure of the multi-field coupling mechanism’s symmetric surface when a fluid is acting upon it. Finally, this article verifies the accuracy and reliability of the algorithm in solving the incompressible problem by comparing the ellipsoid example solution with the experimental results. In the future, this method can transfer data from fluid to solid during strong coupling, avoiding some errors and influences caused by interface data transfer. At the same time, relying on the mathematical idea of symmetry, this can be transplanted to other disciplines to form a coupling calculation method for multi-field coupling mechanisms.

2. Mathematical Basics for Tensor Operation

According to the work of Lu Y. et al. [35], in order to simplify the mathematical operation process in the article, we define the calculation domain as follows.

Definition 1.

$\exists \mathbf{r}\left(x,y\right)\in \tau ,\tau \subseteq {\mathbb{R}}^2$ and our calculations are performed within τ, which is composed of closed outer boundary ∑ and closed inner loop L. In the physical sense, we can use τ as the flow field, ∑ as the outer boundary of the flow field, L as the ellipsoid in the flow field and $\mathbf{r}$ as the position vector of the ellipsoid profile L in ∑.

3. Generalized Operator and Fundamental Solution for Navier–Stokes Equations

In order to obtain the integral expression of Navier–Stokes equations, we need to transform Navier–Stokes equations into a tensor form and solve them [45,46], so the generalized tensor operator and fundamental solution for Navier–Stokes equations are needed.

3.1. Navier–Stokes Equations

One of the fundamental equations of fluid mechanics, the Navier–Stokes equations (N–S equations), refer to the equation of motion that describes the conservation of momentum in a viscous fluid. They have a crucial significance in fluid mechanics since resolving viscous fluid is essentially resolving Navier–Stokes equations. As partial differential equations, their solution is very difficult and complicated [47].

Before the solution idea or technology is further developed and breakthroughs are made, the accurate solution can only be obtained in some straightforward special cases of flow problems. However, in some cases, it is possible to simplify these equations to obtain an approximate solution. To obtain the integral expression of Navier–Stokes equations, we need to transform Navier–Stokes equations into a tensor form as follows.

Lemma 1.

Let $\Omega \subset {\mathbb{R}}^2$, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{V}\in {\mathcal{T}}_N\left(0,1\right)$. Navier–Stokes equations can be represented as the following equations:

$$\frac{\partial \rho }{\partial t}+\left(\nabla ·\rho \mathbf{V}\right)=0;$$

$$\rho \frac{\partial \mathbf{V}}{\partial t}+\left(\nabla ·\left(\nabla \rho \mathbf{V}\right)\right)=\left(\nabla ·\mathbf{T}\right),$$

(1)

where the tensor form of stress $\mathbf{T}$ is

$$\mathbf{T}=-\mathbf{I}p+\mu \nabla \mathbf{V},\mathbf{T}\in {\mathcal{T}}_M\left(1,1\right),$$

where $\mathbf{I}$ is the unit tensor, p is scalar pressure, ρ is the density of the medium and $\mathbf{V}$ is the velocity vector. In addition, the potential characteristics of given streamlines, which are pressure p and vorticity $\Omega \in {\mathbb{R}}^n$, are determined by the given boundary conditions:

$${\mathbf{V}}_{L_i}=0;\left|\left|{\mathbf{V}}_{\sum }\right|\right|=V_{\infty }.$$

The proof is omitted.

In addition, it is necessary to emphasize the basic fact that vector velocity $\mathbf{V}$ and vorticity $\Omega$ are solutions to the main problem of vector analysis [48]:

$$\left(\nabla ·\mathbf{V}\right)=q;\left(\nabla ·\Omega \right)=0,$$

(2)

where

$$q=-\frac{\partial \mathbf{V}}{\partial t}$$

is the heaving of the fluid.

3.2. Generalized Operator

According to the Navier–Stokes equations in tensor form given in Lemma 1, to solve these equations, the generalized operator should be derived and analyzed first.

In classical field theory, the continuum mechanics are analyzed using vector and field theory in Cartesian coordinates. Therefore, we can use the operator ∇ to perform tensor analysis, and the following generalized operator could be obtained.

Lemma 2.

Let $\Omega \subset {\mathbb{R}}^2$, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{a}\in {\mathcal{T}}_N\left(0,1\right)$. There is

$$\nabla \left(\nabla ·\mathbf{a}\right)=\left(\nabla ·{\nabla }^{*}\mathbf{a}\right)=\Delta \mathbf{a}+\left[\nabla \times \left[\nabla \times \mathbf{a}\right]\right],$$

(3)

where $\nabla \mathbf{a}$ and ${\nabla }^{*}\mathbf{a}$ are conjugate tensors, and they can be expressed as:

$$\nabla \mathbf{a}=\mathbf{i}\frac{\partial \mathbf{a}}{\partial x}+\mathbf{j}\frac{\partial \mathbf{a}}{\partial y};$$

$${\nabla }^{*}\mathbf{a}=\mathbf{i}\nabla a_x+\mathbf{j}\nabla a_y,$$

where $\mathbf{i}$ and $\mathbf{j}$ are basic vectors of the two-dimension system. The proof is omitted.

Based on the analysis above, we use the steady incompressible generalized potential to analyze the generalized operator in the following theorem.

Theorem 1.

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of Ω, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{V}\in {\mathcal{T}}_N\left(0,1\right)$, $\forall \Omega \in {\mathcal{T}}_N\left(0,1\right)$. There is

$$\mathbf{VV}+\frac{p}{\rho }-\nu \nabla \mathbf{V}\equiv \mathbf{VV}+\frac{p}{\rho }+\nu [\mathbf{I}\times \Omega ]={\nabla }^{*}\Psi .$$

(4)

The generalized potential Ψ is the scalar basic solution of partial differential Equation (4):

$$\left(\nabla ·{\nabla }^{*}\Psi \right)\equiv \nabla \left(\nabla ·\Psi \right)=0.$$

Proof of Theorem 1.

From Equations (2) and (3), inserting Equation (3) into Equation (2), when $\Omega \subset {\mathbb{R}}^2$, K is the auxiliary space of $\Omega$, and manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{V}\in {\mathcal{T}}_N\left(0,1\right)$, $\forall \Omega \in {\mathcal{T}}_N\left(0,1\right)$, we can obtain

$$\nabla \left(\nabla ·\mathbf{V}\right)=\left(\nabla ·{\nabla }^{*}\mathbf{V}\right)=\nabla q,$$

$$\nabla \left(\nabla ·\Omega \right)=\left(\nabla ·{\nabla }^{*}\Omega \right)=0.$$

(5)

Regardless of the physicochemical changes of object surfaces, only the fluid properties are considered, which are q = 0, and the generalized operator (3) allows us to obtain a conservation of momentum in a conservative form:

$$\left(\nabla ·\left\{\mathbf{VV}+\frac{p}{\rho }-\nu \nabla \mathbf{V}\right\}\right)\equiv \left(\nabla ·\left\{\mathbf{VV}+\frac{p}{\rho }+\nu [\mathbf{I}\times \Omega ]\right\}\right)=0.$$

Then, introducing the vector potential to the above equation, we can obtain

$$\mathbf{VV}+\frac{p}{\rho }-\nu \nabla \mathbf{V}\equiv \mathbf{VV}+\frac{p}{\rho }+\nu [\mathbf{I}\times \Omega ]={\nabla }^{*}\Psi ,$$

and the proof is completed. □

3.3. Fundamental Solution of Generalized Operator

In this section, we will solve the generalized operator to obtain its fundamental solution using the previous definitions and analysis. Then, the integral expression of Navier–Stokes equations can be obtained in conjunction with tensor analysis after establishing the fundamental solution.

Lemma 3.

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of Ω, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{a}\in {\mathcal{T}}_N\left(0,1\right)$, $\forall \Gamma \in {\mathcal{T}}_N\left(1,1\right)$. There is

$$[\nabla \times [\mathbf{a}\times \Gamma ]]=\left(\nabla ·\left(\mathbf{a}·\Gamma \right)\right)-\left(\mathbf{a}·\nabla \right)\Gamma -\Gamma \left(\nabla ·\mathbf{a}\right)+\mathbf{a}\left(\nabla ·\Gamma \right).$$

(6)

Proof. 

In the two-dimensional case, for any vector $\mathbf{a}=a_x\mathbf{i}+a_y\mathbf{j}$ and $\Gamma ={\Gamma }_x\mathbf{i}+{\Gamma }_y\mathbf{j}$, there is

$$[\nabla \times [\mathbf{a}\times \Gamma ]]=\left[\nabla \times \mathbf{k}\left(a_x{\Gamma }_{\mathbf{y}}-a_y{\Gamma }_{\mathbf{x}}\right)\right]$$

$$=\mathbf{i}{\Gamma }_{\mathbf{y}}\frac{\partial a_x}{\partial y}+\mathbf{i}{a}_x\frac{\partial {\Gamma }_{\mathbf{y}}}{\partial y}-\mathbf{i}{\Gamma }_{\mathbf{x}}\frac{\partial a_y}{\partial y}-\mathbf{i}{a}_y\frac{\partial {\Gamma }_{\mathbf{x}}}{\partial y}-\mathbf{j}{a}_x\frac{\partial {\Gamma }_{\mathbf{y}}}{\partial x}-\mathbf{j}{\Gamma }_{\mathbf{y}}\frac{\partial a_x}{\partial x}+\mathbf{j}{a}_y\frac{\partial {\Gamma }_{\mathbf{x}}}{\partial x}+\mathbf{j}{\Gamma }_{\mathbf{x}}\frac{\partial a_y}{\partial x}$$

$$=\left(\mathbf{i}{\Gamma }_{\mathbf{y}}\frac{\partial a_x}{\partial y}+\mathbf{j}{\Gamma }_{\mathbf{x}}\frac{\partial a_y}{\partial x}+\mathbf{i}{\Gamma }_{\mathbf{x}}\frac{\partial a_x}{\partial x}+\mathbf{j}{\Gamma }_{\mathbf{y}}\frac{\partial a_y}{\partial y}\right)-\left(\mathbf{i}{\Gamma }_{\mathbf{x}}\frac{\partial a_x}{\partial x}+\mathbf{j}{\Gamma }_{\mathbf{y}}\frac{\partial a_y}{\partial y}+\mathbf{j}{\Gamma }_{\mathbf{y}}\frac{\partial a_x}{\partial x}+\mathbf{i}{\Gamma }_{\mathbf{x}}\frac{\partial a_y}{\partial y}\right)$$

$$-\left(\mathbf{j}{a}_x\frac{\partial {\Gamma }_{\mathbf{y}}}{\partial x}+\mathbf{i}{a}_y\frac{\partial {\Gamma }_{\mathbf{x}}}{\partial y}+\mathbf{i}{a}_x\frac{\partial {\Gamma }_{\mathbf{x}}}{\partial x}+\mathbf{j}{a}_y\frac{\partial {\Gamma }_{\mathbf{y}}}{\partial y}\right)+\left(\mathbf{i}{a}_x\frac{\partial {\Gamma }_{\mathbf{y}}}{\partial y}+\mathbf{j}{a}_y\frac{\partial {\Gamma }_{\mathbf{x}}}{\partial x}+\mathbf{i}{a}_x\frac{\partial {\Gamma }_{\mathbf{x}}}{\partial x}+\mathbf{j}{a}_y\frac{\partial {\Gamma }_{\mathbf{y}}}{\partial y}\right)$$

$$=\left(\nabla ·\left(\mathbf{a}·\Gamma \right)\right)-\left(\mathbf{a}·\nabla \right)\Gamma -\Gamma \left(\nabla ·\mathbf{a}\right)+\mathbf{a}\left(\nabla ·\Gamma \right),$$

where $\mathbf{i}$ and $\mathbf{j}$ are basic vectors of the two-dimension system. And the proof is completed. □

Based on Equation (6), the fundamental solution of the generalized operator could be established in the following theorem.

Theorem 2.

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of Ω, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \Gamma \in {\mathcal{T}}_N\left(1,1\right)$, $\forall \phi \in {\mathcal{T}}_N\left(0,0\right)$. There is

$$\nabla \left(\nabla ·\Gamma \right)=\Delta \Gamma +[\nabla \times [\nabla \times \Gamma ]]=\mathbf{I}\Delta \phi =0,$$

(7)

where Γ is the fundamental solution of the generalized operator.

Proof of Theorem 2.

According to the Cauchy–Riemann condition:

$$\frac{\partial \phi }{\partial x}=\frac{\partial \psi }{\partial y};\frac{\partial \phi }{\partial y}=-\frac{\partial \psi }{\partial x}.$$

(8)

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of $\Omega$, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\phi \in {\mathcal{T}}_N\left(0,0\right)$, $\psi \in {\mathcal{T}}_N\left(0,0\right)$ and $\mathbf{k}$ be the basic vector of auxiliary space K. Equation (8) can be represented in vector form:

$$\nabla \phi =[\nabla \times \mathbf{k}\psi ];\nabla \psi =-[\nabla \times \mathbf{k}\phi ].$$

(9)

The conjugate analytic functions $\phi$ and $\psi$ are basic solutions of the Laplace equation.

$$\phi =\frac{1}{2\pi }ln\left|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}\right|;$$

$$\psi =\frac{1}{2\pi }arctan\frac{y-y_0}{x-x_0},$$

(10)

where $\mathbf{r},{\mathbf{r}}_0\in \tau$, $\tau \subseteq {\mathbb{R}}^n$, and $\tau$ is the domain within outer boundary ∑. So, tensor

$$\Gamma \left(\left|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}\right|\right)=\mathbf{I}\phi -[\mathbf{I}\times \mathbf{k}\psi ]$$

(11)

is conservative.

Therefore,

$$\nabla \left(\nabla ·\Gamma \right)=\nabla \left(\nabla ·\left\{\mathbf{I}\phi -[\mathbf{I}\times \mathbf{k}\psi ]\right\}\right)=\nabla \left(\nabla ·\mathbf{I}\phi \right)-\nabla \left(\nabla ·\left[\mathbf{I}\times \mathbf{k}\psi \right]\right)$$

$$=\mathbf{I}\Delta \phi -\left[\mathbf{I}\times \Delta \mathbf{k}\psi \right]+\left[\nabla \times \left\{\left[\nabla \times \mathbf{I}\phi \right]-\left[\nabla \times \left[\mathbf{I}\times \mathbf{k}\psi \right]\right]\right\}\right]$$

$$=\mathbf{I}\Delta \phi -\left[\mathbf{I}\times \Delta \mathbf{k}\psi \right]+\left[\nabla \times \mathbf{I}\left(\nabla ·\mathbf{k}\psi \right)\right]=\mathbf{I}\Delta \phi +\left[\mathbf{I}\times \left\{\nabla \left(\nabla ·\mathbf{k}\psi \right)-\Delta \mathbf{k}\psi \right\}\right]$$

$$=\mathbf{I}\Delta \phi +\left[\mathbf{I}\times \left[\nabla \times \left[\nabla \times \mathbf{k}\psi \right]\right]\right]=\mathbf{I}\Delta \phi +\left[\mathbf{I}\times \left[\nabla \times \nabla \phi \right]\right]=\mathbf{I}\Delta \phi =0,$$

and the proof is completed. □

4. Direct Integral Expression of Generalized Operator for Navier–Stokes Equations

Here, to obtain the boundary integral equation, we re-express the generalized operator above in integral form. In addition, to facilitate calculation and analysis, we choose the sum of outer boundary $\Sigma$ and closed inner loop L, which is $\Sigma +L$, as the research domain $\partial \tau$.

According to the work of Lu Y. et al. [35], the following lemma could be obtained.

Lemma 4.

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of Ω, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{a}\in {\mathcal{T}}_N\left(0,1\right)$, $\forall \Gamma \in {\mathcal{T}}_N\left(1,1\right)$, and $\mathbf{n}$ be the normal vector of the boundary $\partial \tau$, so

$$\left(\mathbf{n}·\left[\nabla \times \left[\mathbf{a}\times \Gamma \right]\right]\right)=\left(\frac{{\partial }^{*}\mathbf{a}}{\partial n}·\Gamma \right)-\left(\nabla ·\mathbf{a}\right)\left(\mathbf{n}·\Gamma \right)-\left(\frac{{\partial }^{*}\Gamma }{\partial n}·\mathbf{a}\right)+\left(\nabla ·\Gamma \right)\left(\mathbf{n}·\mathbf{a}\right).$$

(12)

It can be proven by the basic vectors in a two-dimensional Cartesian space [35], and the proof is omitted.

Then, we apply the classical process of isolating singular points to integration based on a generalized Green formula. According to Gauss’s law, we can obtain the following lemma.

Lemma 5.

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of Ω, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{a}\in {\mathcal{T}}_N\left(0,1\right)$, $\forall \Gamma \in {\mathcal{T}}_N\left(1,1\right)$, and $\mathbf{n}$ be the normal vector of the boundary $\partial \tau$. We have

$$\underset{\tau }{\int }\left(\nabla ·\left[\nabla \times \left(\mathbf{a}·\Gamma \right)\right]\right)d\tau =\underset{\partial \tau }{\oint }\left(\mathbf{n}·\left[\nabla \times \left(\mathbf{a}·\Gamma \right)\right]\right)d\sigma .$$

(13)

The proof is omitted.

Then, we perform generalized integration in computational domain $\partial \tau$, and the results are as follows.

Theorem 3.

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of Ω, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{a}\in {\mathcal{T}}_N\left(0,1\right)$, $\forall \Gamma \in {\mathcal{T}}_N\left(1,1\right)$, and $\forall \phi \in {\mathcal{T}}_N\left(0,0\right)$. There is

$$\underset{L_{\epsilon }}{\oint }\left(\mathbf{a}·\frac{\partial \Gamma }{\partial n}\right)d{l}_{\epsilon }=\mathbf{a}.$$

(14)

Proof of Theorem 3.

Combining Equations (12) and (13), there is

$$\underset{L_{\epsilon }}{\oint }\left(\mathbf{n}·\left[\nabla \times \left(\mathbf{a}·\Gamma \right)\right]\right)d{l}_{\epsilon }=\underset{L_{\epsilon }}{\oint }\left\{\left(\frac{{\partial }^{*}\mathbf{a}}{\partial n}·\Gamma \right)-\left(\mathbf{a}·\frac{{\partial }^{*}\Gamma }{\partial n}\right)\right\}d{l}_{\epsilon }+\underset{L_{\epsilon }}{\oint }\left\{\left(\mathbf{n}·\mathbf{a}\right)\left(\nabla ·\Gamma \right)-\left(\nabla ·\mathbf{a}\right)\left(n·\Gamma \right)\right\}d{l}_{\epsilon }.$$

(15)

Performing vector analysis:

$$\underset{L_{\epsilon }}{\oint }\frac{{\partial }^{*}\mathbf{a}}{\partial n}d{l}_{\epsilon }=\underset{L_{\epsilon }}{\oint }\left\{\frac{\partial \mathbf{a}}{\partial n}+\left[\mathbf{n}\times \left[\nabla \times \mathbf{a}\right]\right]\right\}d{l}_{\epsilon }.$$

(16)

So,

$$\frac{{\partial }^{*}\mathbf{a}}{\partial n}=\frac{\partial \mathbf{a}}{\partial n}+\left[\mathbf{n}\times \left[\nabla \times \mathbf{a}\right]\right].$$

Then,

$$\underset{L_{\epsilon }}{\oint }\frac{{\partial }^{*}\Gamma }{\partial n}d{l}_{\epsilon }=\underset{L_{\epsilon }}{\oint }\left\{\frac{\partial \Gamma }{\partial n}+\left[\mathbf{n}\times \left[\nabla \times \Gamma \right]\right]\right\}d{l}_{\epsilon }.$$

(17)

According to the definition of tensor $\Gamma$, we can obtain the following formula from Equation (11):

$$\frac{\partial \Gamma }{\partial n}=\mathbf{I}\frac{\partial \phi }{\partial n}-\left[\mathbf{I}\times \mathbf{k}\right]\frac{\partial \psi }{\partial n}.$$

(18)

From the definition of function $\psi$ (see Equation (10)), there is

$$\frac{\partial \psi }{\partial n}=0.$$

Then,

$$\frac{\partial \Gamma }{\partial n}=\mathbf{I}\frac{\partial \phi }{\partial n}.$$

According to the generalized integral principle:

$$\left|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}\right|=\epsilon .$$

From the definition of function $\phi$ (see Equation (10)), we obtain

$$\frac{\partial \phi }{\partial n}=\frac{\partial \phi }{\partial \epsilon }=\frac{\partial }{\partial \epsilon }\left(\frac{ln\left|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}\right|}{2\pi }\right)=\frac{1}{2\pi \epsilon }.$$

Performing generalized integration [49,50]:

$$\underset{L_{\epsilon }}{\oint }\left(\mathbf{a}·\frac{\partial \Gamma }{\partial n}\right)d{l}_{\epsilon }=\underset{L_{\epsilon }}{\oint }\left(\mathbf{a}·\mathbf{I}\frac{\partial \phi }{\partial n}\right)d{l}_{\epsilon }$$

$$=\underset{L_{\epsilon }}{\oint }\left(\mathbf{a}·\frac{\partial \phi }{\partial n}\right)d{l}_{\epsilon }={\int }_0^{2\pi }\mathbf{a}·\frac{1}{2\pi \epsilon }·\epsilon d\theta =\mathbf{a}·\frac{1}{2\pi }{\int }_0^{2\pi }d\theta =\mathbf{a},$$

and the proof is completed. □

By deriving and integrating the above formulas, and calculating in $\partial \tau$, we can establish the integral expression of Navier–Stokes equations as the following equation.

Theorem 4.

Let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of Ω, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\forall \mathbf{a}\in {\mathcal{T}}_N\left(0,1\right)$, $\forall \Gamma \in {\mathcal{T}}_N\left(1,1\right)$, ${\mathbf{r}}_0\in \tau$ and $\mathbf{n}$ be the normal vector of the boundary $\partial \tau$. There is

$$\mathbf{a}\left({\mathbf{r}}_{\mathbf{0}}\right)=-\underset{\tau }{\int }\left(\nabla q·\Gamma \right)d\tau +\underset{\partial \tau }{\oint }\left[\left(\frac{\partial \mathbf{a}}{\partial n}·\Gamma \right)-\left(\frac{\partial \Gamma }{\partial n}·\mathbf{a}\right)\right]d\sigma -\underset{\partial \tau }{\oint }\left(\left[\mathbf{n}\times \left[\nabla \times \mathbf{a}\right]\right]·\Gamma \right)d\sigma .$$

(19)

In the simplest case where the velocity divergence is zero (q = 0) in the region, the non-thermally conductive incompressible fluid has an integral representation, such as a velocity vector:

$$\mathbf{V}\left({\mathbf{r}}_{\mathbf{0}}\right)=\underset{\partial \tau }{\oint }\left[\left(\frac{\partial \mathbf{V}}{\partial n}·\Gamma \right)-\left(\frac{\partial \Gamma }{\partial n}·\mathbf{V}\right)\right]d\sigma -\underset{\partial \tau }{\oint }\left(\left[\mathbf{n}\times \left[\nabla \times \mathbf{V}\right]\right]·\Gamma \right)d\sigma ,$$

(20)

where Σ is the outer boundary of the flow field τ, L is the closed inner loop in the flow field τ and $\partial \tau$ is the the research domain, which is the sum of Σ and L.

Proof of Theorem 4.

Integrate Equations (14), (16) and (17), and put the result into Equation (15):

$$\underset{\tau }{\int }\left(\nabla ·\left[\nabla \times \left(\mathbf{a}·\Gamma \right)\right]\right)d\tau =-\mathbf{a}+\underset{\partial \tau }{\oint }\left\{\left(\frac{\partial \mathbf{a}}{\partial n}·\Gamma \right)-\left(\mathbf{a}·\frac{\partial \Gamma }{\partial n}\right)\right\}d\sigma$$

$$+\underset{\partial \tau }{\oint }\left\{\left[\mathbf{n}\times \left[\nabla \times \mathbf{a}\right]\right]-\left[\mathbf{n}\times \left[\nabla \times \Gamma \right]\right]\right\}d\sigma +\underset{\partial \tau }{\oint }\left\{\left(\mathbf{n}·\mathbf{a}\right)\left(\nabla ·\Gamma \right)-\left(\nabla ·\mathbf{a}\right)\left(n·\Gamma \right)\right\}d\sigma .$$

(21)

According to Green’s theorem [21], we can obtain

$$\underset{\partial \tau }{\oint }\left\{\left(\nabla ·\mathbf{a}\right)\left(n·\Gamma \right)-\left(\mathbf{n}·\mathbf{a}\right)\left(\nabla ·\Gamma \right)\right\}d\sigma =\underset{\tau }{\int }\left\{\left(\nabla \left(\nabla ·\mathbf{a}\right)·\Gamma \right)-\left(\nabla \left(\nabla ·\Gamma \right)·\mathbf{a}\right)\right\}d\tau .$$

(22)

With Equation (7), there is

$$\nabla \left(\nabla ·\Gamma \right)=0.$$

Similarly, an integer of zero will be zero, so

$$\underset{\tau }{\int }\left(\nabla \left(\nabla ·\Gamma \right)·\mathbf{a}\right)d\tau =0.$$

(23)

Integrating Equations (22) and (23), we can obtain

$$\underset{\partial \tau }{\oint }\left\{\left(\mathbf{n}·\mathbf{a}\right)\left(\nabla ·\Gamma \right)-\left(\nabla ·\mathbf{a}\right)\left(n·\Gamma \right)\right\}d\sigma =-\underset{\tau }{\int }\left(\nabla \left(\nabla ·\mathbf{a}\right)·\Gamma \right)d\tau .$$

(24)

This is because $\left(\nabla ·\left[\nabla \times \mathbf{X}\right]\right)=0,\underset{\tau }{\int }\left(\nabla ·\left[\nabla \times \left(\mathbf{a}·\Gamma \right)\right]\right)d\tau =0$, where $\mathbf{X}$ is any vector or tensor (see Equation (7)). Then, through Equation (21), the following formula is easy to obtain.

$$\underset{\tau }{\int }\left(\nabla ·\left[\nabla \times \left(\mathbf{a}·\Gamma \right)\right]\right)d\tau =\underset{\partial \tau }{\oint }\left(\frac{\partial \mathbf{a}}{\partial n}·\Gamma \right)d\sigma -\underset{\partial \tau }{\oint }\left(\mathbf{a}·\frac{\partial \Gamma }{\partial n}\right)d\sigma -\mathbf{a}$$

$$+\underset{\partial \tau }{\oint }\left\{\left[\mathbf{n}\times \left[\nabla \times \mathbf{a}\right]\right]-\left[\mathbf{n}\times \left[\nabla \times \Gamma \right]\right]\right\}d\sigma -\underset{\tau }{\int }\left(\nabla \left(\nabla ·\mathbf{a}\right)·\Gamma \right)d\tau =0.$$

(25)

So,

$$\mathbf{a}\left({\mathbf{r}}_{\mathbf{0}}\right)=-\underset{\tau }{\int }\left(\nabla \left(\nabla ·\mathbf{a}\right)·\Gamma \right)d\tau +\underset{\partial \tau }{\oint }\left\{\left(\frac{\partial \mathbf{a}}{\partial n}·\Gamma \right)+\left[\mathbf{n}\times \left[\nabla \times \mathbf{a}\right]\right]\right\}d\sigma -\underset{\partial \tau }{\oint }\left\{\left(\mathbf{a}·\frac{\partial \Gamma }{\partial n}\right)+\left[\mathbf{n}\times \left[\nabla \times \Gamma \right]\right]\right\}d\sigma .$$

(26)

We perform integral expression on the solution of operator $\nabla \left(\nabla ·\mathbf{a}\right)=\nabla q$ of any vector $\mathbf{a}$ in the region of the control boundary [51]:

$$\mathbf{a}\left({\mathbf{r}}_{\mathbf{0}}\right)=-\underset{\tau }{\int }\left(\nabla q·\Gamma \right)d\tau +\underset{\partial \tau }{\oint }\left\{\left(\frac{\partial \mathbf{a}}{\partial n}·\Gamma \right)+\left[\mathbf{n}\times \left[\nabla \times \mathbf{a}\right]\right]\right\}d\sigma -\underset{\partial \tau }{\oint }\left\{\left(\mathbf{a}·\frac{\partial \Gamma }{\partial n}\right)+\left[\mathbf{n}\times \left[\nabla \times \Gamma \right]\right]\right\}d\sigma .$$

(27)

So, in the simplest case where the velocity divergence is zero ($q=0$), let $\Omega \subset {\mathbb{R}}^2$, K be the auxiliary space of $\Omega$, manifold $N=\left(\Omega ,\mathcal{S}\right)$, $\Gamma \in {\mathcal{T}}_N\left(1,1\right)$, $\mathbf{V}\in {\mathcal{T}}_N\left(0,1\right)$, ${\mathbf{r}}_0\in \tau$, and $\mathbf{n}$ be the normal vector of the boundary $\partial \tau$. The non-conductive incompressible fluid has an integral expression, such as a velocity vector:

$$\mathbf{V}\left({\mathbf{r}}_{\mathbf{0}}\right)=\underset{\partial \tau }{\oint }\left\{\left(\frac{\partial \mathbf{V}}{\partial n}·\Gamma \right)+\left[\mathbf{n}\times \left[\nabla \times \mathbf{V}\right]\right]\right\}d\sigma -\underset{\partial \tau }{\oint }\left\{\left(\mathbf{V}·\frac{\partial \Gamma }{\partial n}\right)+\left[\mathbf{n}\times \left[\nabla \times \Gamma \right]\right]\right\}d\sigma ,$$

(28)

where $\Sigma$ is the outer boundary of the flow field $\tau$, L is the closed inner loop in the flow field $\tau$ and $\partial \tau$ is the the research domain, which is the sum of $\Sigma$ and L.

According to analysis above, Equations (27) and (28) could numerically implement the closed-loop integral by using their singular integral and weak singular integral. Then, we can use the above conclusions to accurately calculate the aerodynamic characteristics on object surfaces. The proof is completed. □

5. Discrete Form of the Boundary Integral Equation

Equation (28) must be discretized before being applied to the calculation since the aerodynamic surface and its calculation process are complex. Effective discretization can reduce the time and space cost of the solution process and improve the accuracy of the calculation results. Additionally, BEM also conducts discretization processing to simplify problems by reducing their dimension [52]. According to the work of Xu K. et al. [53], we can discretize the equations obtained above to apply them to BEM.

According to the analysis above, BEM can be used to determine the aerodynamic properties of an ellipsoid’s symmetric boundary. Meanwhile, using numerical implementation could make the boundary problem simple to solve. During the solving process, boundary discretization is primary, and this helps to reduce the difficulty of the aerodynamic problem solving process. According to the work of Xu K. et al. [53], the linear interpolation polynomial we would choose here is

$$x=\frac{1-\xi }{2}{x}_1+\frac{1+\xi }{2}{x}_2,$$

and it is a typical discretization method. We will conduct subsequent analysis on the basis of this discrete result.

Based on the Equation (28) and the results in Theorem 2, the following parameters required for the solution are obtained with introducing the concept of vector potential [54] into the discrete process.

$$u_{L_i}=\frac{1}{2}\sum _{k=1}^K\left(-\frac{1}{\upsilon \rho }{p}_k{S}_k^i+{\omega }_k{T}_k^i\right)+\frac{3}{2}\sum _{m=1}^M\left(-\frac{1}{\upsilon \rho }{p}_m{S}_m^i+{\omega }_m{T}_m^i\right)+\frac{3}{2}{C}_i;$$

$$v_{L_i}=\frac{1}{2}\sum _{k=1}^K\left(-\frac{1}{\upsilon \rho }{p}_k{T}_k^i-{\omega }_k{W}_k^i\right)+\frac{3}{2}\sum _{m=1}^M\left(-\frac{1}{\upsilon \rho }{p}_m{T}_m^i-{\omega }_m{W}_m^i\right)+\frac{3}{2}{D}_i;$$

$$u_{{\Sigma }_i}=\frac{3}{2}\sum _{k=1}^K\left(-\frac{1}{\upsilon \rho }{p}_k{S}_k^i+{\omega }_k{T}_k^i\right)+\frac{1}{2}\sum _{m=1}^M\left(-\frac{1}{\upsilon \rho }{p}_m{S}_m^i+{\omega }_m{T}_m^i\right)+\frac{1}{2}{C}_i;$$

$$v_{{\Sigma }_i}=\frac{3}{2}\sum _{k=1}^K\left(-\frac{1}{\upsilon \rho }{p}_k{T}_k^i-{\omega }_k{W}_k^i\right)+\frac{1}{2}\sum _{m=1}^M\left(-\frac{1}{\upsilon \rho }{p}_m{T}_m^i-{\omega }_m{W}_m^i\right)+\frac{1}{2}{D}_i,$$

(29)

where K and M represent the number of elements on L and $\Sigma$, respectively; p and $\omega$ represent pressure and vortex, respectively [55]. Furthermore,

$$C_i=\sum _{m=M_2+1}^{M_3}{R}_m^i+\sum _{m=M_1+1}^{M_2}\left(R_m^i-\frac{1}{\upsilon }{u}_m^2{F}_m^i\right)+\sum _{m=M_3+1}^M\left(R_m^i+\frac{1}{\upsilon }{u}_m^2{F}_m^i\right);$$

$$D_i=\sum _{m=M_2+1}^{M_3}{Q}_m^i+\sum _{m=M_1+1}^{M_2}\left(Q_m^i+\frac{1}{\upsilon }{u}_m^2{F}_m^i\right)+\sum _{m=M_3+1}^M\left(Q_m^i-\frac{1}{\upsilon }{u}_m^2{F}_m^i\right);$$

$$S_j^i=n_{xj}{F}_j^i-n_{yj}{G}_j^i;T_j^i=n_{yj}{F}_j^i-n_{xj}{G}_j^i;W_j^i=n_{xj}{F}_j^i+n_{yj}{G}_j^i;$$

$$R_j^i=-\frac{u_j^2}{\upsilon }{n}_{xj}{F}_j^i+\frac{u_j^2}{\upsilon }{n}_{yj}{G}_j^i-u_j{H}_j^i;Q_j^i=-\frac{u_j^2}{\upsilon }{n}_{yj}{F}_j^i+\frac{u_j^2}{\upsilon }{n}_{xj}{G}_j^i-u_j{H}_j^i;$$

$$F_j^i={\int }_{a_j}^{b_j}\phi d\sigma ;G_j^i={\int }_{a_j}^{b_j}\psi d\sigma ;H_j^i={\int }_{a_j}^{b_j}\frac{\partial \phi }{\partial n}d\sigma .$$

The above discrete result is complex in form, but the core of the algorithm is very simple. That is, the matrix A containing the integration parameters and the vector x containing the boundary conditions are combined to obtain the final result b, as shown in Equation (30).

$$\left\{A\right\}·\left\{x\right\}=\left\{b\right\}.$$

(30)

6. Results

6.1. Analysis of Results

According to above mathematical derivation, we selected the simplest case, that being that there is only one multi-field coupling mechanism (ellipsoid) in the symmetric flow field, and the pressure coefficient is calculated under the angle of attack of $-5^{\circ }$, $0^{\circ }$, $4^{\circ }$, $7^{\circ }$ and ${10}^{\circ }$. The results are as follows.

According to above calculation results, we have obtained the ellipsoid’s symmetric surface pressure distribution curve at different angles of attack. The ellipsoid’s symmetric surface pressure distribution results are shown in Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7. In fluid mechanics, since both the lift and drag coefficients are integral to the pressure coefficient curve along the chord, they can be calculated from the pressure coefficients. At the same time, borrowing the calculation idea in the study of Sumbatyan A. M. [56] and the following integral formula, we can calculate the aerodynamic coefficient $C_a$.

$$C_a=\frac{1}{x_{TE}-x_{LE}}{\int }_{x_{LE}}^{x_{TE}}\left({\overline{p}}_l\left(x\right)-{\overline{p}}_u\left(x\right)\right)dx,$$

(31)

where ${\overline{p}}_l$ is the pressure coefficient of the lower surface, ${\overline{p}}_u$ is the pressure coefficient of the upper surface and $x_{LE}$ is the leading edge position of the airfoil. $x_{TE}$ is the trailing edge position of the airfoil. Furthermore, considering the direction of pressure induction, the lift coefficient and drag coefficient can be obtained by decomposing the aerodynamic coefficient. To compare with the results of Ejeh C. J. [57], we selected the angles of attack as −5${}^{\circ }$ and 7${}^{\circ }$. Under the premise that other conditions are consistent, after calculation, the errors between the method in this article and the results of Ejeh C. J. (Figure 8) are

$${\Delta }_{\alpha =-5^{\circ },C_l}=\frac{|-0.0227-(-0.0230\left)\right|}{0.0230}=1.3\%;{\Delta }_{\alpha =-5^{\circ },C_d}=\frac{|0.035-0.036|}{0.036}=2.8\%,$$

$${\Delta }_{\alpha =7^{\circ },C_l}=\frac{|0.0402-0.0390|}{0.0390}=3.1\%;{\Delta }_{\alpha =7^{\circ },C_d}=\frac{|0.0371-0.0362|}{0.0362}=2.5\%.$$

Obviously, these errors are all around 3%. Meanwhile, given the influence of data estimation, environmental conditions and computer power, we consider these errors acceptable. This shows that the results have been experimentally verified and are in good agreement with the experimental results. These results can help us quickly determine the pressure distribution at a specific angle of attack, which is of great significance for improving the aerodynamic characteristics of the multi-field coupling mechanism in a symmetric flow field.

At the same time, we have also obtained the influence curve of the number of ellipsoid computing grids, which is of great significance for us in determining the optimal number of computing grids and the optimal solution. As shown in Figure 9, when the number of computing grids is 2500, the calculation result deviates significantly from the actual situation. The result is much better when the number of computing grids is 3000 to 3500, but there is still a certain deviation from the actual result. When the number of computing grids is 5000, the result still has some errors with the actual experimental data. When the number of computing grids is 8000, the calculation results are consistent with the actual experimental results, and the error is within the acceptable range. The optimal number of computing grids should be between 6000 and 8000.

At the same time, in order to test the correctness of this method, we have made a comparison with the method of Belotserkovsky O. M. [58]. The comparison is carried out on a symmetric sphere which is in a flowing fluid with Re = 40 or 100. The results are as follows:

From Figure 10, under the condition of Re = 100, the calculation result of this method is highly consistent with the calculation result of Belotserkovsky O. M., and the error is small. Under the condition of Re = 40, the calculation result of this method is consistent with the calculation of Belotserkovsky O. M. From Figure 10, the biggest gap appears at $\theta ={120}^{\circ }$. If the values from Figure 10 are taken, the prediction accuracy could be calculated as

$${\Delta }_{Re=40,\theta ={120}^{\circ }}=\frac{|0.9-0.92|}{0.9}=2.22\%<2.5\%.$$

It is within the acceptable range. In general, this method’s calculation results are in line with the actual experimental results. In other words, this method is correct in calculating the vortex distribution on the sphere’s surface.

Meanwhile, we also calculated for a symmetric ellipsoid with an aspect ratio of 5 ($\lambda =5$) and compared the results with the results of Lyona’s experiment [59]. The comparison is performed with $\alpha =0^{\circ }$, Re = $2.04\times {10}^6$, $\lambda =5$, and the results are shown in Figure 11.

It can be seen from the above results that this method has very little error with the classical method in calculating the pressure distribution on an ellipsoid’s surface when $\overline{x}\in \left(0,0.5\right)$. It is basically consistent with the result of the classical method. When $\overline{x}\in \left(0.5,1\right)$, the results of this method have a big error when used with the classical method, which might be due to the separation of fluids. Taking the values from Figure 11, the prediction accuracy could be calculated as

$${\Delta }_{\overline{x}=1}=\frac{|0.22-0.21|}{0.22}=4.55\%<5\%.$$

For the cause of the error of the above results, on the one hand, there will be a certain error in the experiment itself, whether it is in the setting of the experimental conditions or the collection of experimental data. On the other hand, since there will be numerical dissipation in the calculation process, and the BEM has certain errors when solving nonlinear problems and complex vortex structures [60], the above errors will occur between the two results. They are within the acceptable range, so we believe that the two results are consistent.

At the same time, we also calculated the lift and drag coefficients of a symmetric ellipsoid with varying angles of attack and aspect ratios $\lambda$. The results are shown in Figure 12 and Figure 13.

From the results in Figure 12 and Figure 13, under the condition of a positive angle of attack, when $\lambda =3.5$, the ellipsoid is subjected to greater lift coefficient and drag coefficient. In addition, the influence of $\lambda$ on lift is much greater than the influence on drag. The above calculation results have been verified through experiments.

6.2. Discussion of Results

According to the above results, whether it is the calculation of surface pressure at different angles of attack or the calculation of surface pressure under different calculation grids, the method proposed in this paper is consistent with the previous experimental results. It shows the correctness of the new calculation methods developed in this paper. This paper makes it possible to use BEM calculation methods for incompressible fluids by proposing fundamental tensor solutions, which are entirely new methods. In addition, it is easier to couple solids, achieving weak coupling calculations similar to strong coupling.

7. Conclusions

This article first derives and proposes a boundary integral equation mathematically according to a generalized operator and a generalized Green integral theorem and discredits it. Then, we apply it to the algorithm basis and numerical realization of the boundary element method (BEM) to determine the aerodynamic characteristics of the multi-field coupling mechanism (symmetric ellipsoid) when it is in viscosity incompressible fluid flow. The discrete geometric boundary algorithm created in this article calculates the integrals (such as single-layer or two-layer potential form), and singular solutions can be obtained from the integral expression of the Navier–Stokes equations.

Compared with experimental results and results from other methods, the method in this article dramatically simplifies the calculation process, improves calculation efficiency, and can determine the unknown value on the boundary of a two-dimensional region under certain situations. At the same time, this method transforms the nonlinear boundary value problem into a linear boundary integral equation, which can effectively reduce the dimension of high-dimension problems, reduce the difficulty of calculation, and effectively solve complex dimensional mechanical problems. The calculated results of this method are in good agreement with the experimental results, and the error is small.

At the same time, the calculation method in this article for aerodynamic characteristics on a multi-field coupling mechanism’s symmetric surface can also be applied to other fluid analysis, such as boundary analysis of a complex flow field, aerodynamic characteristic analysis of an aircraft, fluid–solid interaction analysis and other symmetric flow problems. In the future, this method can transfer data from fluid to solid during strong coupling, avoiding some errors and influences caused by interface data transfer. At the same time, relying on the mathematical idea of symmetry, this can be transplanted to other disciplines to form a coupling calculation method for multi-field coupling mechanisms. Meanwhile, the new method proposed in this article can effectively reduce the difficulty of analyzing these problems, which is of great significance for the theoretical development and engineering application of fluid mechanics.