1 Introduction

Neumann's method^[1] has been used in mathematical physics for a long time. The boundary problems for Laplace's equation have been shown to have a solution for the convex boundary of the domain^[1] with simple- and double-layer potentials. Once the resolving integral equations were built, Neumann^[1] suggested to use the successive approximation method for the numerical calculation (Neumann's method).

Fredholm^[2] has proven the existence of the solutions for an arbitrary smooth boundary. Developing the general theory of integral equations expressed by

(1)

on the [a, b] interval, Fredholm^[3] showed that the resolvent of the written equation was an analytical function within the circle |λ|＜λ₁, where λ₁ was the smallest eigenvalue of Eq. (1) . This result can be used to prove the convergence of the Neumann series for an arbitrary shape. Indeed, if the resolvent is regular within the circle |λ|< 1 while the point λ=1 on the boundary is not its pole, then

(2)

which means that the resolvent has no singularity for any |χ|≤1. Inserting Eq. (2) into Eq. (1) and re-expanding the series corresponding to the function φ(x) with respect to χ, we form a convergent series for the point λ=1. Let r=2. Then, we can obtain the Neumann series. This method has been used by Kantorovich and Krylov^[4] to solve integral equations.

With the potential method, the elasticity problems can be reduced to singular integral equations^[5]. Therefore, the transformation of Eq. (2) is applicable, and the Neumann series converges. These results have been published in Refs. ^[6]-^[8]. Additionally, it has been established that Neumann's numerical analysis can be interpreted only in terms of partial differential equations and without integral equations.

Differential equations allow for the natural variational formulation of elasticity problems. The corresponding function is symmetric. Once Korn's inequality has been proven, the variational issues of the elasticity theory can be resolved^[9]. The question has been raised concerning whether it might be possible to demonstrate the convergence of the Neumann series in elasticity problems without integral equations. The solution was offered in Ref. ^[10], and the embedding theorems in the Sobolev-Slobodetsky space have been extensively used ^[11-12].

This paper is organized as follows. Section 1 presents the mechanics of the shell theory in the form of differential equations and variational principles. Section 1 does not identify the areas where the solutions of the problem are being sought. Starting with Section 2, the needed precise mathematical definitions are introduced. The application of the Neumann series to the shell theory problems in a variational setting requires the proof of Korn's inequality. This inequality is proven in Refs. [13] and [14]. Building on this proof, Section 2 proves the convergence of the Neumann series for the manifold shells with holes within the distribution space. Section 3 proves the convergence of the Neumann series for mixed boundary value problems, while Section 4 proves the convergence acceleration of the successive approximations for the shells with holes.

2 Defining boundary problem for shells with holes

Let us consider a shell with the thickness h and a manifold middle surface D⁺ (see fig 1) . The outer contour S_o = S₁+S₂ consists of two parts, i.e., S₁ and S₂. There are no displacements for S₁, and there are no forces for S₂. The contours of the internal m holes are designated by the letters S_i^h(i =1, 2, …, m). Let the set of all contours be designated by the letter S. The surface is defined with the principal curvature orthogonal coordinates α₁ and α₂. Additionally, to simplify the following discussion, S is assumed to coincide with the line α₁ =const.

Let us define three vectors for an arbitrary point on the surface x∈ D⁺, where e₁ and e₂ designate the tangential plane, while e₃ designates the normal direction. The load is applied only to S. At the point of the contour where the vector s defines the position of the tangent line, the applied load reduces to the forces Q₁⁰, Q₂⁰, and Q₃⁰ directed along the vectors e₁, e₂, and e₃, respectively, and the moment Q₄⁰=M₃ directed along the vector s.

Under the static loads Q₁⁰, Q₂⁰ , Q₃⁰, and Q₄⁰, the points of the shell move within the local system of the coordinates by the vector u=(u₁, u₂, u₃). Let us demonstrate the variational statement of the linear static problem for the thin shell theory to the extent that will be required for the further proof. The introduced concepts are shown in Appendix A^[15-16]. There are several similar statements of the problem in question that do not differ much. Further results are given in Ref. [17].

The potential strain energy of the thin elastic shell is defined by

(3)

where

and

(4)

In the above equations, A₁ and A₂ are the coefficients of the first quadratic form of the surface, and E and υ are Young's modulus and Poisson's number for the shell material, respectively. The parameters R₁ and R₂ represent the main radii of the curvature. The turning angles of the normal are expressed via the displacements of the surface points as follows:

(5)

The work of the forces Q⁰=(Q₁⁰, Q₂⁰, Q₃⁰, Q₄⁰) on the contour S is designated with the integral

The static analysis of a shell with holes includes finding a vector u that can provide the minimum of the total energy

(6)

and converge to zero so as to coincide with the turning angle of the normal on the contour S₁, i.e.,

(7)

The variation of the functional ∏(u) is equal to zero, and the fixing of the boundary S₁ corresponds to three differential equations and four boundary conditions, i.e.,

(8)

(9)

(10)

The linear operators T_k⁺(u)(k=1, 2, 3, 4) include the first-order derivatives of u₁ and u₂ and the third-order derivative of u₃. Physically, these operators denote the limiting values of the internal forces when the point x∈ D⁺ approaches S from the domain D⁺. The details of the operators R(u) and T_k(u) are only provided in Appendix A because they are not needed for our purposes.

Let us designate E(u) as the set of the vector functions satisfying Eqs. (8) , (10) , and (11) and start looking for the vector in this set for which Eq. (9) is true.

The solution is sought with the following rule. The manifold D⁺ is supplemented with the simply connected domains D_i^- to form the simply connected shell D = D⁺ +∑D_i⁻. Within the domain D, the continuous vector u meeting the boundary conditions must be found by

(12)

The negative sign in the index in Eq. (12) denotes the limiting value when the direction of the vector aims from the interior of the domain D_i^- toward the boundary S.

The formulated problem can be solved by the successive approximation

(13)

(14)

Note that

(15)

The bilinear form W_D+^∗(w, u) in Eq. (15) is obtained from the quadratic function W_D-(u) by the replacement of each squared term dependent on one vector with one-half of the scalar product of the same functions dependent on two vectors, e.g., replacing (ε₁(u)+ε₂ (u))² with (ε₁(u)+ ε₂(u )) ((ε₁(v) +ε₂ (v))/2.

The other replacements are made similarly. Therefore, there is the convergence (14) in the norm induced by the scalar product W_D+^∗(w, u) if Eq. (16) holds for any vector u that is continuous on S and meets Eq. (8) when

(16)

For the remainder of the paper, c is used to represent various positive constants.

Indeed, the limiting equality follows from Eqs. (12) and (14) . For convergence, it is sufficient to satisfy

which is equivalent to Eq. (16) . The next section describes the conditions when Eq. (16) is satisfied.

3 Proving convergence of Neumann series (14)

Before proceeding, it must be determined that whether the area contained within the contour S₀ can be considered as a Lipschitz domain. For two-dimensional shells, a Lipschitz domain requires that the domains D⁺ and D^- can be broken down into a finite number of parts, each of which has a planar projection with the property

(17)

In Eq. (17) , the left side of the inequality denotes the Cartesian distance between the points z^* and z, which are situated on the surface, and the right side represents a distance between the projections of the same points.

The answer to the question is sought within the space E(u), which is induced by the scalar product W_D^*(w, u). The space E(u) is a result of the supplement with respect to

of the space K(u)=(K₁(u₁) , K₁ (u₂) , K₂(u₃) ), where u₁, u₂, and u₃ are components of the displacement vector u. The manifold K₁(u) is a space of continuous functions with the first-order derivatives, and K₂(u) is a space of functions with the second-order continuous derivatives.

Let us define the distribution class H within the domain D as follows:

where H¹ and H² are the Sobolev spaces with the respective norms

In the above equations, D^mu designates the norm of the mth-order distributional derivative within the space L₂(D). The norm in H(D) is calculated by

It is clear that

(18)

The 4-vector space u with the components u₁, u₂, u₃, and υ₁ is defined on the contours S_i (i=1, 2) such that

In the Sobolev-Slobodetsky class, the norm is set by

where $\hatu\left( t \right)$ is the Fourier transformation of the function u(x). The norm in H(S) is defined by

(19)

If the contours S_i (i=1, 2) meet the Lipschitz condition, then there is a trace of vectors from H(D) to H(S), satisfying^[18-19]

(20)

Because the vector u satisfies Eq. (8) within the domain D⁺, Korn's inequality is satisfied, i.e.,

(21)

Many authors, beginning with Friedrichs^[20], have been interested in proving Korn's inequality for various problems. The history is covered in Ref. ^[9]. Ciarlet and Mardare^[21] have proven it for elastic shells. From Eqs. (20) and (21) , we have

(22)

It is also the case that, within the domain D^-, the function W_D^-(u) reaches its minimum with any possible extension of the vector u from S inside the domain. The Lipschitz contours allow for the extension^{[12, 17]}. Then,

(23)

Combining Eq. (22) with Eq. (23) , we obtain the proof of the inequality (16) and the convergence of the Neumann series.

If the Green matrix G_ik is known, each member of the Neumann series in Eq. (14) can be explicitly calculated as a simple-layer potential from the domain D with the formula

(24)

The multiplication operator in Eq. (24) means that the matrix is multiplied by a column. Each first three columns of the matrix G_ik are vectors used to solve Eq. (8) when the zeros on the right-hand side of the equations are replaced by the three columns of the following matrix:

where δ designates the δ-function. The fourth column of G_ik corresponds to the following three functions:

(25)

The formulas in Eq. (25) assume that the coordinate line α₁=, const. coincides with S. If the condition is not met, the formulas in Eq. (25) become longer, but the components of the fourth column remain linear from the elements of the first three columns. Because the simple-layer potential satisfies Eq. (8) and is continuous on S while the operator T_k(w_n) having a jump on the contour S equals T_k( w_n-1}), Eq. (24) is an implementation of the Neumann series, as described by Eq. (14) .

The calculations for the problems of the elasticity theory for infinite domains (the Green tensor includes only the fundamental solution) and the calculations of the bending plates (when it is possible to explicitly write down the formulas for the Green tensor) are reduced to the quadratures of the singular integrals.

4 Application of Neumann series to mixed boundary problems of shell theory

The Neumann series (14) can also be built for other shell theory boundary problems, e.g., the mixed boundary problem, where within the domain D^-, the vector u satisfies Eq. (8) , a part of the contour S₁ is fixed, and S₂ is under the forces Q⁰=(Q₁⁰, Q₂⁰, Q₃⁰, Q₄⁰) (see fig 2) .

On the contour S₁, we identify a 4-vector line u with the components u₁, u₂, u₃ , and υ₁. Let us set a mixed problem of the shell theory, where

(26)

As in Section 2, we supplement the domain D^- with the surface D⁺ to the boundary S^*, and start searching for the vector u in the extended area D^-+D⁺=D such that

(27)

For the problem (27) , the Neumann series (14) is

(28)

Applying Eq. (15) to the contour S=S₁ +S₂, we observe that the convergence of the series (28) is assured if Eq. (16) is satisfied. This inequality is proven based on the assessment of W_D⁺(u) and W_D⁻(u) with the Hilbert-space norm $\left\| u \right\|$_S.

Problem (28) has boundary conditions on the contour S₁, similar to the terms of the Neumann series for the elliptic equation. Therefore, Eq. (28) cannot be formulated within the Hilbert space, but it does allow for an interpretation in the Banach space of the Sobolev functions with a negative index.

Let us define the norm of the vector u by

(29)

The second component can be found with Eq. (19) , while the first component is the norm of the 4-vector u in the space

which is calculated according to the following rule:

(30)

If Eq. (28) is interpreted as follows:

(31)

for any vector v∈ H(D) and

(32)

then, with Eqs. (31) and (32) , we can obtain the following estimate for the vector of any approximation u=w_n based on the embedding theorem for the domain D⁺:

(33)

Since the 1960s, many articles^[22-23] have been dedicated to the proof of the inequalities of the type in Eq. (33) for elliptic systems.

Korn's inequality^[14] for the domain D^- gives

(34)

Since Eqs. (33) and (34) are equivalent to Eq. (16) , the convergence of the Neumann series is proven in the energy norm.

Each term of the series (26) can be written explicitly if Green's adjoint matrix (equivalent of the double-layer potential) is introduced as follows:

(35)

The first integral is the same as that of Eq. (24) , while the matrix T^∧G_ik is

The elements of the above matrix are generated by

(36)

where g_k=G_1k, G_2k, G_3k (k=1, 2, 3).

Indeed, the first equation in Eq. (27) results from Green's tensor definition, and the second equation is due to a displacement jump w_n−1^∗ in the double-layer potential on S₂, while the operator T_k(w_n) is continuous (equivalent of Liapunov's theorem^[24] for the potential theory). The first integral in Eq. (35) is continuous, and has a jump equal to T_i(w_n−1) on the contour S₂ (equivalent to that of the single-layer potential). Thus, the solution to the mixed boundary problem of the theory can be found for an arbitrary boundary if Green's matrix is known for one standard domain.

As it has been demonstrated for the plates in Ref. [8], the explicit solution to the problem (27) per Eq. (28) with Eqs. (35) and (36) may also seem efficient in practical calculations because many authors have devised a decent method for the construction of fundamental solutions and the Green functions for elliptic systems. Hormander^[25] covered the fundamental solutions to elliptic systems (especially useful for the equations with constant coefficients) based on the Fourier transformation. Antosik et al.^[26] used the fundamental solutions as a series of differentiable functions with preliminary chosen singularities. Based on the distributions, Chen et al.^[27] and Lukasiewicz^[28] presented the fundamental solutions and the Green functions g_k=G_1k, G_2k, G_3k(k=1, 2, 3) for a lot of interesting cases (spherical, conical, cylindrical, and shallow shells). Ultimately, Eq. (35) provides a solution to the singular pseudo differential equation

where w^*=w if s∈ S₁, and w^*=T_i⁺ (w) if s∈ S₂.

The kernel of the above equation is calculated per Eq. (36) . This equation solves the mixed problem of the shell theory within the boundary element method^[29] when only the values of the functions on the contours S₁ and S₂ are used for the calculations. Many problems benefit from the successive approximation method because, unlike the boundary element method, it requires much less memory for the inversion of the high order matrices.

5 Improvement in convergence of Neumann series for shells with holes

Equations (24) and (35) demonstrate that the elastic stress-strain state of complex-shaped shells can be described based on the calculations of relatively simple 1-connected bodies. At the same time, the numerical analysis of multiply connected shells demonstrates that there is a strong stress concentration near the holes. Thus, the series of Eq. (14) should not converge too fast. The parameter c in Eq. (16) should be close to zero for small holes. The problem of improving the convergence of the series (14) for the shells with small holes arises.

Let us review a similar problem involving the need to find the eigenvalues λ_i and the eigen vector-functions q_iz for the system of Eq. (37) , i.e.,

(37)

Formula (15) leads to the question of the existence of the discrete spectrum and relevant eigenfunctions in the Hilbert space. Since this problem is solvable, any vector within the space E(u) can be expanded in a convergent series with the orthonormal system q_i, e.g.,

(38)

According to the symmetry of the form W^∗_D(w, u) and Eq. (16) , all the eigenvalues are real and are bigger than one, i.e.,

(39)

The transformation of each term of the series (38) via Eq. (14) , with (37) , leads to

(40)

This equality shows that after n iterations, the series (38)

(41)

In view of Eq. (39) , the series (41) asymptotically tends to the first eigenfunction of the system (37) , as in a geometrical progression.

(42)

The feature summarized in Eq. (42) is often used in the theory of symmetrical integral equations. It helps to accelerate the calculation process according to the summation of the infinite series of the geometrical progression, and is also a criterion for the proper functioning of the algorithm. Another advantage of Eq. (14) is that it allows for the possible use of parallel calculations while finding the integrals (24) and (35) .

6 Conclusions

This study proves that Neumann's method of successive approximations can be applied to the numerical analysis of complex-shaped shells when it is known how to calculate the relatively simple standard samples. For instance, a cylindrical shell can be used as a standard, while a cylindrical shell with holes of arbitrary shapes can be used as a complex shell.

The numerical implementation of Neumann's algorithm is associated with the calculations of the boundary integrals. This fact significantly reduces the computer memory space and calculation time required for the analysis compared with the method of finite two-dimensional elements. Another advantage of the proposed algorithm includes the possibility of speeding up the calculation process by selecting the rows with the geometrical progression properties. Thus, Neumanns' method is useful for the analysis of the elastic shells with multiple holes and unsmooth contour lines.

Appendix A

The thin-walled shell is a three-dimensional body with the following Cartesian coordinates of any X_i point:

(A1)

where x_i(α_k) is the Cartesian coordinate of a point on the middle surface D lying on the same normal with the point x_i. The parameter z represents the distance between the points x_i and x_i along the vector of the outer unit normal η_i(α_k)_i. The distance between the boundary points on the same normal defines the shell thickness h(α_k)=h. The shell is considered to be thin if the ratio h/L is small. The size L means the maximum distance between the points. The middle surface D is parameterized in the Lagrangian coordinate α_k. Let us define the first quadratic form of the middle surface as follows:

(A2)

where a_mn is the metric tensor of the middle surface. The form describes the distance between the infinitely close surface points. The differences between the Lagrangian coordinates of these points are dα₁ and dα₂. The summation over the repeated indices is performed in this formula and in all further ones.

It follows from Eq. (A1) that the differential squared of the Euclidean distance between the shell points provided with the coordinates α₁, α₂, and z are calculated by

(A3)

When m and n equal 1 or 2,

The remaining elements of the metric tensor g_mn of the shell are g₃₃=1, and g₃₁=g₃₂=g₁₃=g₂₃=0.

If the arbitrary vector u_i(α_k) is plotted at any point of the shell, the covariant derivative for it is calculated per the formula

(A4)

The Christoffel symbols are expressed through the metric tensor with the formulas

(A5)

In particular, for the unit vectors in Fig. 1,

The following designations are introduced in the orthogonal coordinates of the principal curvatures:

The Christoffel symbols that do not equal zero are written as follows:

where R₁ and R₂ mean the principal radii of curvature on the middle surface.

Let us assume that the points of the shell X_i move to the position X_i^* under the load q.

Based on Eq. (A1) , we have

(A6)

In Eq. (A6) , u_i(α_k) is the displacement of the points on the middle surface, and η^*_i is the new position of the normal to the middle surface.

Let us accept the Kirchhoff hypothesis that the fiber perpendicular to the middle surface before deformation remains to be perpendicular to the deformed surface and does not change its length. Then, we have

(A7)

where v_i allows for the calculation of the shell deformation based on the formulas for the three-dimensional medium

(A8)

The stress tensor σ _ij elastically (Hooke's law) depending on the deformation tensor is associated with each point of the shell, i.e.,

(A9)

The internal energy of the thin shell can be presented as follows:

(A10)

Substituting Eqs. (A8) and (A9) into Eq. (A10) and taking into account Eq. (6) of the main text, we can obtain the full energy of the shell by

The Lagrangian principle of the minimum full energy claims that the real displacement vector minimizes the functional Π. Note that the functional W rewritten in the coordinates of the main curvatures coincides with the functional W_D(u)=W_ε(u)+W_χ(u) (see Eqs. (3) and (4) ).

The necessary condition for the minimum δΠ =0 and the transition to the independent variations δ u₁, δ u₂, δ u₃, and δ θ₁ after partial integration of the obtained expression results in the following three equilibrium equations over the vector u (u_i, i=1, 2, 3) for the internal points of the domain:

(A11)

where the forces T^ij and M^ij are expressed through the vector u and its derivatives per

Thus, the left sides of the three equations in Eq. (A11) identically describe the operator R(u) in Eq. (8) .

During the integration by parts of the functional δ Π=0, the curvilinear integral appears

The operators included in P¹ look as follows:

(A12)

The functions σ_t and τ_t mean the normal curvature and the geodesic torsion, while the vectors υ, n, and t define the tangential normal, the normal to the middle surface, and the tangent line in the contour points. The remaining operators are calculated by contraction of the tensors, i.e.,

The given expressions demonstrate that the four lines of the boundary condition operator T_λ(u) from Eqs. (9) and (10) of the first section can be identified with the four expressions on the right sides of Eq. (A12) .

At the same time, it is established that the variational statement of the static problem for the thin shells, which aims to minimize the quadratic functional

equals solving Eq. (8) with the boundary conditions in Eqs. (9) -(11) .