Patchwise Mapping Method for Solving Elliptic Boundary Value Problems Containing Multiple Singularities ()
1. Introduction
It has been introduced to solve multiple crack problems by using various numerical methods. First, converting the multiple crack problems into Fredholm integral equation using two elementary solutions was introduced in [2] . A numerical method by using both Fredholm integral equation method and the weighted residual method was introduced in [3] . Numerical methods based on Galerkin approximation such as the finite element methods, boundary element methods, and meshless method were also applied to solve them [3] - [10] .
In this paper, we solve the elliptic boundary value problems with multiple singularities based on the mapping method [1] . But, the mapping technique proposed is not able to catch singularities emerging at multiple locations in a domain. In order to resolve the drawback, we introduced the enriched Isogeometric Analysis (IGA) in [11] . In the paper [11] , we approximate the solution on the small circular zone centered at the crack tip or point singularity by enriching the finite approximation space generated by the singular mapping introduced in the mapping method. However, it is hard to evaluate the inverse functions of the singular mapping, and the NURBS mapping so that tracking the domains of integrals whose integrand is involved both B-spline function from the singular mapping and NURBS function from the NURBS geometrical mapping, is an additional work. Also, the conditional number of the stiffness matrix could be an issue for the enriched IGA. In order to alleviate these problems, we design the geometrical map having multiple singularities by using the singular mappings only. To do that, we divide the physical domain into multiple patches which may have a singularity for each, and then design the geometrical maps by the mapping methods for each patch having a singularity. Here, we consider the design of control points on the interface between/among patches. Because this is related to the smoothness of the global basis functions. Also, we modify the B-spline functions whose supports include the interface between/among them due to the compatibility condition. In this paper, the potential of the mapping method proposed with multiple patches regarding to handling the multiple fatigue-cracks propagation in various types of plate will be shown by solving the elliptic boundary value problems with multiple singularities or cracks.
In Section 1 and 2, we briefly review definitions and terminologies that are used throughout this paper. We follow those in the book [12] , and we thus refer to these texts for details. And then we introduce the mapping method that generates singular functions by using B-spline or NURBS in Section 3. In Section 4, we introduced the patchwise mapping method by solving elliptic boundary value problems containing multiple singularities. Finally, the conclusions is in Section 5.
2. Nomenclature
In this section, we introduce B-spline, NURBS, and design geometrical mappings referring to [12] .
2.1. B-Splines
A knot vector
is a nondecreasing sequence of real numbers in the parameter space
, and the components
are called knots. An open knot vector of order
is a knot vector that satisfies
in which the first and the last
knots are repeated.
The B-spline functions
of order
corresponding to the knot vector
are piecewise polynomials of degree p which are constructed recursively by the Cox-de Boor recursion formula:
for
For example, the piecewise quadratic polynomial B-spline functions
corresponding to the open knot vector
are depicted in Figure 1.
The B-spline functions are useful in design as well as finite element analysis because they have the following properties: variation diminishing, convex hull, non-negativity, piecewise polynomial, compact support, and partition of unity.
A B-spline curve is defined as follows:
where
and
are control points that make B-spline functions draw a desired curve as shown in Figure 2(a).
Let
be an open knot vector and let
and
, respectively, be the polynomial degree and order of B-spline functions
. Then a B-spline surface is defined by
Figure 1. B-spline functions
,
of order
corresponding to the knot vector
.
Figure 2. (a) B-spline curve and control points on the open knot vector
. (b) B-spline basis functions corresponding to the B-spline curve shown in (a). (a) B-spline curve and control points; (b) B-spline basis functions.
where
and
are control points that make a bidirectional control net as shown in Figure 2(b).
2.2. Nonuniform Rational B-Spline (NURBS)
In this section, we review the non-uniform rational B-splines (NURBS), NURBS curve and surface, and primary properties of NURBS.
2.2.1. NURBS Curve
A pth-degree NURBS curve is defined by
(1)
where the
are the control points, the
are the weights,
, and the
are the pth-degree B-spline basis functions defined on the nonperiodic (and non-uniform) knot vector
We assume that
, and
for all i. Setting
(2)
allows us to rewrite Equation (1) in the form
(3)
The
are the rational basis functions; they are piecewise rational functions on
.
2.2.2. NURBS Surface
A NURBS surface of degree
in the
direction and degree
in the
direction is a bivariate vector-valued piecewise rational function of the form
(4)
The
form a bidirectional control net, the
are the weights, and the
and
are the nonrational B-spline basis functions defined on the knot vectors
Introducing the piecewise rational basis functions
the surface Equation (4) can be written as
An example of the NURBS surface is shown in Figure 3.
2.3. Weak Solution in Sobolev Space
Let
be a connected open subset of
. We define the vector space
Figure 3. An example of B-spline surface with control net in three dimensional space.
to consist of all those functions
which, together with all their partial derivatives
of orders
, are continuous on
. A function
is said to be a
-function. If
is a function defined on
, we define the support of
as
For an integer
, we also use the usual Sobolev space denoted by
. For
, the norm and the semi-norm, respectively, are
Suppose we are concerned with an elliptic boundary value problem on a domain
with Dirichlet boundary condition
along the boundary
. Let
The variational formulation of the Dirichlet boundary value problem can be written as: Find
such that
(5)
where
is a continuous bilinear form that is
-elliptic ( [13] ) and
is a linear functional. The solution to (5) is called a weak solution which is equivalent to the strong (classical) solution corresponding elliptic PDE whenever u is smooth enough. The energy norm of the trial function u is defined by
Let
,
be finite dimensional subspaces. Since the NURBS basis functions do not satisfy the Kronecker delta property, in this paper we approximate the nonhomogenuous Dirichlet boundary condition by the least squares method as follows:
such that
We can write the Galerkin form (a discrete variational equation) of (5) as follows: Given
, find
, where
, such that
which can be rewritten as: Find the trial function
such that
(6)
2.4. Variational Formulation of Equilibrium Equations of Elasticity
In elasticity, the displacement field is denoted by
, and the stress field is denoted by
. Let
be the strain field. Then the strain-displacement and the stress-strain relations are given by
(7)
respectively, where
is the differential operator matrix,
and
is the
symmetric positive definite matrix of material constants. Material constants are classified by the property of the material. For an isotropic elastic body,
Here,
where E is the Young’s modulus of elasticity and
is Poisson’s ratio.
The equilibrium equations of elasticity are
(8)
where
is the vector of internal sources representing the body force per unit volume.
The equilibrium Equation (8) can be expressed in terms of the displacement field
through the relations (7). Then we consider the following system of elliptic differential equations in terms of the displacement field,
(9)
subject to the boundary conditions,
(10)
where
,
is a unit vector normal to the boundary
of the domain
.
For the Galerkin approximation to the equilibrium equations in terms of the displacement field (9), the variational form of (9) through (10) is:
find the vector
such that
on
, and
(11)
where
(12)
The finite element approximation of the solution of (11) is to construct approximations of each component of the vector
.
3. Mapping Method
We introduce a geometrical mapping from the parameter space
to
that generates singular basis functions [1] .
3.1. B-Spline Curves That Generates Singular Basis Functions
In particular, we first show how a B-spline curve
handles effectively one-dimensional singularities. Let
be an open knot vector of order
. Then the B-spline functions
corresponding to
are
(13)
Here,
,
, are also called the Bernstein polynomials of degree
. Let
be control points, for a constant
. Then the B-spline geometrical mapping
(14)
(15)
(16)
maps the parameter space
onto the physical space
and its inverse is
Thus, the approximation space
, where
is an integer greater than or equal to
and
are the Bernstein polynomials (B-spline functions) of degree
and contain the following singular as well as smooth functions:
In other words, the geometrical mapping
is able to generate the singularity of type
, where
.
For example, if
, then the Bernstein polynomials of degree 2 are
and
are control points. Then the geometrical mapping obtained by these control points and its inverse, respectively, are
Suppose
where
are the Bernstein polynomials corresponding the the open knot vector
of order 5, then
contains
. Hence the approximation space
for isogeometric analysis contains
3.2. NURBS Surface That Generates Singular Basis Functions
The argument which is the construction of geometrical mapping that generates singular basis functions, can be extended to NURBS surface from the parameter space
to
. To end this, we construct a NURBS surface to deal with monotone singularity of type
, where q is a rational number with
,
is a piecewise smooth function,
is the polar coordinates. the construction of the NURBS surface from
to the unit disk in [1] is generalized in this section. We refer to this reference for the details.
We now consider a NURBS surface from the parameter space
to the physical domain
. Consider the knot vectors:
(17)
where
,
,
,
, and
.
Let m and
be the number of knots in the knot vectors
and
, respectively. Also, let k and
be
and
, respectively. Here, if the function to be approximated has a singularity of type
with
, where
, then the polynomial degree of B-spline functions corresponding to
is
.
Let
,
be the B-splines corresponding to the knot vector
and let
,
be the B-splines corresponding to the knot vector
. Here, the B-spline functions
,
, corresponding to the open knot vector
are the Bernstein polynomials of degree
.
Consider the control points
and the weights
for
,
, that are listed in Table 1. We now construct a NURBS surface from the parameter space
onto
as follows:
(18)
Here
,
,
, are NURBS basis functions defined by
(19)
where
Since
is the Bernstein polynomial and
unless
, substituting Equations (13) into (19) the NURBS surface mapping (18) becomes
Table 1. Control points
and weights
.
Now, we derive the derivative of the mapping
by using formulas in Chapter 4.5 in [12] .
Let
where
is the numerator of
.
Denoting
, the derivative of
can be expressed as
Then
(20)
Solving the Equation (20) for
, we obtain
(21)
We employ the lemma below from Chapter 3 in [12] in order to determine
and
.
Lemma 1 Let
, and
. Then
with
and the knot vector corresponding to
is
Applying the lemma 1 into (21), we have
where
The derivative of the total weight function
, also, can be described in detail by substituting the Bernstein polynomial into
.
3.3. Numerical Example of Mapping Method for Solving an Isotropic Elasticity Containing Single Singularity
The mapping method proposed was implemented in the paper [1] , and the paper showed that the mapping technique using NURBS geometrical mappings constructed by an unconventional choice of control points are effective for numerical solutions of elliptic boundary value problems containing a single singularity. In this subsection, we solve an elastic problem containing a singularity to show that the proposed method is also applicable to implement elastic problems.
Throughout this paper, we measure the error
of the computed solutions obtained by isogeomtric analysis using the proposed mapping method in the following norms:
· The relative error in the maximum norm in %:
· The relative error in
norm in %:
· The relative error in energy norm in %:
Assuming that the Young’s modulus
and the Poisson’s ratio
in a sector of the unit circle whose the central angle is 270˚, plane strain isotropic elastic body, we consider that the following analytic stress field,
(22)
where
, and
. Then, the stress field (22) satisfies the equilibrium equations of elasticity on the sector shaped domain
. And the displacement field has the singularity of the form
where
is a smooth function.
For the design of the physical domain
, we set
and
in the knot vector (17) so that the open knot vector corresponding to
-direction is as follows:
(23)
We construct the open knot vector corresponding to
-direction using the form of the knot vector
in (17):
(24)
which make Bernstein polynomials in
on the parameter space. We choose
for control points
, and set the other control points as depicted in Figure 4(a). Then the NURBS surface mapping
and the inverse of the design mapping generates the singularity of the form
along the radial direction on the
in Figure 4(a).
In order to enrich the NURBS or B-spline basis functions without failing the structure of the mapping technique, we employ refinements [14] , [15] in the NURBS functions which are used to design the physical domain as depicted in Figure 4(a). In particular, we use p-refinement to enrich the basis functions corresponding to both the open knot vectors (24) and (23). Note that inserting new knots to increase the number of basis functions along the
-direction may cause malfunction regarding the production of singular functions [11] .
Figure 4(b) and Table 2 depict the relative errors of the displacement u and v in the maximum norm(blue and red line, respectively), and in the
norm (green and purple, respectively) versus the number of degrees of freedom. Figure 4(c) and Table 3 depict that the relative errors of the stress field
in the maximum norm versus the number of degrees of freedom. Both Figure 4(b) and Figure 4(c) show that the proposed mapping method captures the singularity effectively, and follows the theoretical results in [1] . Figure 5 exhibits the relative errors of the strain energy of the stress field.
4. Patchwise NURBS Mapping Method for Solving Elliptic Boundary Value Problems Containing Multiple Singularities
In the case of that a physical domain contains multiple cracks, we re-design the
Figure 4. (a) The NURBS surface
maps from the parameter space
to the sector shaped domain
. The coordinates of the primary control points
are described. (b) Relative errors of the displacement field
in the maximum norm and
norm for (22) (c) Relative errors of the stress field
in the maximum norm for (22). (a) The scheme of the NURBS surface
that generates singular functions; (b) Relative errors of the displacement field
; (c) Relative errors of the stress field
.
Table 2. The relative errors (%) of the elasticity containing singularity on the sector shaped elastic material (22): relative errors (%) of displacement field.
Table 3. The relative errors (%) of the elasticity containing singularity on the sector shaped elastic material (22): The computed strain energy and the relative errors (%) of the stress field
. The row “
” indicates the exact values.
Figure 5. Relative errors of the strain energy of the stress field (22) on the sector shaped elastic body.
geometrical mapping by using both standard NURBS mappings and the proposed mappings. Then it simplifies things to describe these sub-domains by different patches. We describe how to construct the set of global basis functions crossing interfaces between patches. Throughout the following examples, we show that the patchwise mapping method is effective in dealing with a problem containing multiple singularities.
First, We apply the mapping method for the elliptic boundary value problems with multiple singularities of type
, where
, and
’s are smooth functions.
Example 1. Let
where
(25)
Then
is the analytic solution of the Poisson equation:
(26)
and has two singularities at
and
.
4.1. Patchwise NURBS Mappings and Interfaces
In Example 1, we divide the physical domain into three patches:
Let
’s are physical patches. We construct NURBS geometrical mappings
and
that generate singularities of the type
and
, respectively. They are also the design maps from the parameter space
to the physical patch
, for each
. To build up
we use the following knot vectors: For
,
(27)
where
,
, and
. For
,
where
.
In the design mappings, we observe the following:
1) We employ control points and weights from Example 5.3 in [1] to build
, and primary control points are shown in Figure 6(a).
2) We design the NURBS geometrical mapping
that generates a singularity
using the control points as depicted in Figure 6(b).
3) Using the affine transformation we define
In Figure 6, the quasi-physical patch is a physical patch translated.
4) Since a singularity does not appear in the patch
, we employ the standard NURBS design technique to build the mapping
from the parameter space
to
.
4.2. Construction of Global Basis Functions over Interfaces and Approximation Space
Now, we construct an approximation space by using B-spline functions which were used in the design mapping
. First, we consider connectivity among B-spline functions defined on different patches and are nonzero along the interface as depicted in Figure 7(a). To obtain
global basis functions crossing
Figure 6. Primary control points of the NURBS geometrical mapping
and
from
to (a) the physical patch
(b) the quasi-physical patch
which has the singularity at the origin.
is constructed by composition of
with an affine transformation, and
maps from the parameter space,
to the physical space
. (a) Primary control points and design of the physical patch
; (b) Primary control points and design of the quasi-physical patch
.
an interface between two different patches, we merge two B-spline local basis functions defined on different patches that have the same nonzero value on the interface between the two patches. In Figure 8(a),
’s represent intervals corresponding to the interface on the physical domain in Figure 7(a). In
, the index i means the index of the patch belonging to the interval, the other index j indicate the index of patch such that
Figure 7. The control points for NURBS surface with (a) three patches for the Example 1, and (b) three patches for the Example 2, respectively. (a) The control points with two singularities; (b) The control points with two singularities.
Figure 8. (a)
represents the interval along ξ- or η-directions corresponding to the interface in the physical domain. We follow directions of arrows when we set the global index. (b) We merge two B-spline basis functions which are interpolants for each parameter space. (a) Parameter spaces of
’s; (b) Construction of new basis function from two interpolant B-spline basis functions defined on two distinct parameter space.
To construct global basis functions which are nonzero functions on
, we merge the nonzero basis function in
and the nonzero basis function in
, where
such that they are reflection about the interface
in the physical domain. In Figure 8(b), for example, let
and
be B-spline basis functions of
in
and
, in Figure 8(b), respectively. Note that i and j in the bracket
represent indices of patches. So
and
. Let
and
be B-spline basis functions of
such that
Because
·
and
·
,
·
,
·
on the inter face
.
We merge two B-spline basis functions
and
so that we count the new function as one global basis function. The new function has a nonzero value on two distinct patches. Here, we should carefully set
, and we apply the same degree of p-refinement into each parameter space.
For a space with the non-homogeneous boundary condition in Example 1,
(28)
We decompose the space (28) into
and
The finite dimensional subspace i.e. approximation space of the Poisson equation (26) is
(29)
where
1)
and
are the number of B-spline functions in ξ- and η-direction of the patch
, respectively.
2)
are new global basis functions by merging two B-spline functions in
and
, respectively of
and
, respectively.
3)
and
are the set of B-spline basis functions composition with the inverse of NURBS surface mapping
on the physical domain
satisfying homogeneous boundary condition.
4)
and
are the set of B-spline basis functions composition with the inverse of NURBS surface mapping
on the physical domain
satisfying non-homogeneous boundary condition.
Figure 9 shows the relative errors (%) versus DOFs. In Figure 9(a) and Table 4, we enrich the set of basis functions by p-refinement and increase the degree of polynomial
and
up to 14. The DOF is 3061 when
. We can see that the proposed mapping method is effective to capture multiple singularities as well as a single crack or singularity.
Example 2. Let
be the unit disk, where
’s are minor sectors whose central angles are 120˚ for each
, and
Figure 9. The relative error (%) in the maximum norm, the L2-norm, and the energy norm of the computed solutions of the Poisson equation (a) (26) (b) (31) versus number of degrees of freedom, respectively. (a) Rel error vs number of degrees of freedom; (b) Rel error vs number of degrees of freedom.
Table 4. The relative errors (%) of the Poisson equation on the rectangle (26): The computed strain energy and the relative errors (%) of the approximate solution
. The row “
” indicates the exact values.
where
(30)
Then
solves the following elliptic boundary value problem:
(31)
and the solution
has three crack singularities at
.
In Example 2, we divide the unit disk into three sectors
including each crack as depicted in Figure 7(b).
Then, we build a design mapping
from the parameter space
to a quasi-physical sector
using the proposed mapping method, for
. Here, we define three physical patches
by using quasi-physical sectors
as follows:
which means that
’s are sectors having the same radii and the central angles as these of
’s through the transformation (30) but the position of the crack tip in
is the origin other than
in
for
. A structural drawing detailed for
and
is shown in Figure 10, and
is designed by rotating
to −90˚. Finally, we define the NURBS geometrical mapping that generates singularities as follows:
Similar to that of Example 1, considering the continuity of the basis functions and the construction of the basis functions on interfaces, we merge two basis
Figure 10. The NURBS geometrical mapping
from
to the quasi-physical patch
whose crack tip is the origin, generates the singularity of the type
, for
. Once we design the mappings, we update them by composition with transformation so that the final NURBS geometrical mapping preserving the mapping technique in
maps from the parameter space,
to the physical sector
, for
. (a) Primary control points and design of the quasi-physical patch
; (b) Primary control points and design of the quasi-physical patch
.
Figure 11. The directions of arrows represent the order of the index t of
along the boundaries
corresponding to interfaces on the physical domain. We need to consider these directions when we set the global index. Both right and left sides colored by red of each parameter spaces are lines corresponding to each cracks in the physical domain.
functions defined on different patches that are nonzero along the interface as depicted in Figure 7(b). Figure 11 shows the order of the index t of
when we set the global index. Similar to the approximation space (29), the finite dimensional subspace of the weak solution of (31) has the form
where
1)
and
are the set of B-spline basis functions composited with the inverse of the NURBS surface mapping
on the physical domain
satisfying the homogeneous boundary condition.
2)
and
are the set of B-spline basis functions composited with the inverse of the NURBS surface mapping
on the physical domain
satisfying the non-homogeneous boundary condition.
Figure 9(b) and Table 5 show that the relative errors (%) in the maximum norm, the
-norm, and the energy norm of the computed solution for equation (31). We can see that the relative errors in the maximum norm and
-norm decrease exponentially, and the error in the energy norm decrease almost linearly. In Figure 9(b), we enrich the basis by p-refinement along
and
-directions, and the number of degrees of freedom and the degree of polynomial were increased up to 14 and 4915, respectively. Figure 12(a) and Figure 12(b) show the graph of the computed solution of Equation (26) and (31), respectively.
5 Conclusions
In this paper, the physical domains of the elliptic boundary value problems containing multiple singularities, were re-designed by the patchwise mapping methods. In the patchwise mapping method, the construction of the approximation space is different from that in the conventional mapping method [16] due to the use of multiple singular functions.
One of the advantages of the patchwise NURBS mapping method including the NURBS mapping technique is to not only generate singular functions but also preserve the properties of IGA. The properties are the followings [15] :
Table 5. The relative errors (%) of the Poisson equation on the unit disk (31): The computed strain energy and the relative errors (%) of the approximate solution
. The row “
” indicates the exact values.
Figure 12. (a) The graph of approximate solution of the equation (26). The DOF was used 2675 with
. Under these parameters, we obtained the rel. maximum norm error 2.232E−08%, rel.
norm error 4.300E−09%, and rel. energy norm error 4.330E−04%. (b) The graph of approximate solution of the equation (31). The DOF was used 4915 with
. Under these parameters, we obtained the rel. maximum norm error 3.022E−06%, rel.
norm error 4.563E−07%, and rel. energy norm error 9.146E−04%. (a) Approximate solution of the equation (26); (b)Approximate solution of the equation (31).
1) The capability of more precise geometric representation of complex objects than conventional Finite Element Methods.
2) Mesh refinement without altering the geometry throughout the refinement process.
Thus, we expect that the patchwise mapping method will be effective for dealing with multiple curved [17] , angled, branched, or radiating cracks. Also, the proposed method can be applied to compute the stress intensity factor and energy release rate in the plate theory [18] [19] . These will be introduced in the subsequent paper.
On the other hand, the drawback of the mapping method is that it is not eligible to use control points and weights imported from Computer Aided Design (CAD) whereas the conventional IGA is available. To overcome this drawback, the approximation space of the standard IGA can be enriched by the mapping method to deal with singularities [11] . The enrichment of IGA by the mapping method is more practical because there are several advantages in the view of engineering designers and IGA users. First, the original design mapping is not needed to be changed. The enriched NURBS approximation space in IGA can generate singular functions through the proposed mapping method on a singular zone. In the mapping method, k- and h-refinement do not produce optimal results. But k-refinement is applicable in the space of NURBS basis in the enriched IGA for improved computational solution. So the enriched IGA for multiple singularities or cracks would be the future work.