The xLab Series of Mathematica-based high precision 2-D elliptic PDE solvers

Michael Quinlan
University College Cork, Ireland
m.quinlan@ucc.ie

Abstract

A basis for solutions to two-dimensional elliptic partial differential equations is formed from closed-form solutions suggested by each segment of a region's boundary in isolation. Each component solution attains harmonic boundary conditions on the associated segment from which it decays into the region.

Harmonic matching techniques are used to compute the combinations of closed-form solutions that deliver harmonic values on a prescribed segment and zero boundary values on all other segments. These combinations permit the coefficients for an arbitrary set of boundary conditions to be quickly determined from their Fourier expansions on the various segments.

Navier's equation of equilibrium elasticity and Mathematica's big-number arithmetic are used as a particularly uncompromising development platform. The usual efficiencies arise when machine numbers are used. The regions may contain cracks and circular and elliptic holes and may either be enclosed by a polygonal boundary or attain a constant state of stress at infinity. The notebook interface caters to users with little experience with Mathematica.

Introduction

This paper gives the basic features of a family of solvers characterized by ElastostaticLab. The approach to solving Laplace's equation is similar with the potential and normal derivative in the role of displacement and traction. These solvers extend the range of  separation of variable techniques and counterpoint numerical techniques such as the Finite Element Method in much the same way as Mathematica's symbolic operators counterpoint their numerical counterparts to endow scientific projects with a hitherto unknown degree of quality assurance.

The Elastostatic Problem

The two dimensional elastostatic problem is to find a displacement field U = U[p] at each point p in a region of the plane which satisfies Navier's equation of elastic equilibrium

Δ w + 
        (1/(1 - 2 ν)) ∇ div[w] + ((Subscript[ρ, 0])/μ) b = 
      0         (1)

with w ≡ U in the region and takes on assigned boundary values for either the displacement U or traction τ at each point on the boundary where  Δ is the Laplacian,  ∇ is the gradient, div is the divergence, ν is Poisson's ratio,  (Subscript[ρ, 0]) is the density and b is the body force.

The traction  τ = τ[p] exerted across a curve with normal n is the vector

\(\(ℰ\/\(\(\ \)\((1 + ν)\)\)\) \((\(\(\ \)\(ν\)\)\/\(\(\ \)\((1 - 2\ ν)\)\)\ \ Tr[∇U]\ ℐ + 
            1\/\(\(2\)\(\ \)\)\ \((∇U\  + ∇U\^T)\))\) . 
        n\)

Here ℰ is Youngs modulus, Tr is the trace operator, ℐ is the two dimensional identity matrix and ∇U  is a 2×2 matrix of displacement components. The traction τ has units of force per unit area.

The boundary of the region is taken to consist of elliptic segments.  To fix ideas, the region exterior to segments 2 and 3 with prescribed tractions and to segment 1 with prescribed displacements is used.

Each segment is the locus of a point p with coordinates a Cos[ϕ] and b Sin[ϕ] relative to the principal axes \(x\&^\) and \(y\&^\) where ϕ is the eccentric angle.  Circles and cracks correspond to a = b and b = 0. Segments may be isolated or in the case of cracks be part of a polygonal boundary. It turns out to be convenient to work with traction boundary values τ[ϕ] via their local normal and tangential curvilinear components functions \(τ[ϕ] . 𝓃\&^[ϕ]\) and \(τ[ϕ] . 𝓉\&^[ϕ]\) respectively. Displacement boundary values U[ϕ] are conveniently dealt with via their local cartesian components functions \(U[ϕ] . x\&^\) and \(U[ϕ] . y\&^\) respectively.

The overall approach is to find a particular solution (Subscript[U, inh]) which satisfies the inhomogeneous equation (1) without necessarily satisfying the boundary conditions. The Fourier expansion of component functions corresponding to discontinuous boundary values on a segment σ converge slowly. In these circumstances,  particular solutions w ≡ (Subscript[U, h, σ]) are constructed to satisfy  (a) the homogeneous equation

\(\(\(		\ \)\(Δw\  + 
        1\/\(1 - 2\ ν\)\ ∇\ div[w] = 
      0\ \ \ \ \ \ \ \ \((2)\)\)\)\)

and  (b) the boundary condition on segment σ. Let

\(\(\(		\)\(u = 
      U - \((U\_inh + ∑\_σ U\_\(h, σ\))\)\)\)\)

and observe that the u ≡ w satisfies the homogeneous equation (2) and that the Fourier series of boundary component functions is enhanced. The function u is called a complementary solution.

Harmonic and Canonical Harmonic Solutions

Solutions (Subscript[u, 1]^σ), (Subscript[u, 2]^σ), … to the homogeneous equation (2) are called canonical harmonic solutions for the segment σ if (a) either the  (local cartesian ) components of the boundary displacements are

({{{1, 0}, {0, 1}, }, {{Cos[ϕ], 0},  ... , {Cos[n ϕ], 
              0}, 
{0, Cos[ϕ]},  ... , {0, 
              Cos[n ϕ]}, }, {{Sin[ϕ], 0},  ... , {Sin[n ϕ], 
              0}, 
{0, Sin[ϕ]},  ... , {0, 
              Sin[n ϕ]}}})         (3)

or the (local curvilinear) components of the boundary tractions are

({{{0, 1}, {0, Cos[ϕ]}, {0, 
              Sin[ϕ]}, 
{1, 0}, {Cos[ϕ], 
              0}, {Sin[ϕ], 0}, }, {{0, Cos[2  ϕ]},  ... , {0, 
              Cos[n ϕ]}, 
{0, 
              Sin[2  ϕ]},  ... , {0, Sin[n ϕ]}, }, {{Cos[2  ϕ], 0},  ... , {Cos[n ϕ], 
              0}, 
{Sin[2  ϕ], 
              0},  ... , {Sin[n ϕ], 0}}})    (4)

depending to whether displacements or tractions are prescribed on segment σ and (b) the boundary values on all other segments are zero where n is the harmonic order appropriate to the required accuracy. The layouts facilitate functional programming [1] which deals with functions {0, Cos[𝒩 ϕ]} of a single list parameter 𝒩 = {1, …, n} rather {0, Cos[2  ϕ]},  ... , {0, Cos[n ϕ]}.

Any complimentary solution u can be written

u = (∑_σ)((Subscript[c, 1]^σ) 
              (Subscript[u, 1]^σ) + (Subscript[c, 2]^σ) 
              (Subscript[u, 2]^σ) + (Subscript[c, 3]^σ) 
              (Subscript[u, 3]^σ) + (Subscript[c, 4]^σ) 
              (Subscript[u, 4]^σ) + ⋯)

where (Subscript[c, 1]^σ), (Subscript[c, 2]^σ), … are permutations of either the Fourier coefficients of the boundary values \(u[ϕ] . x\&^\) and \(u[ϕ] . y\&^\) or of the Fourier coefficients of  \(τ[ϕ] . 𝓉\&^[ϕ]\) and \(τ[ϕ] . 𝓃\&^[ϕ]\) depending on whether displacements or tractions are prescribed on segment σ.

The problem of determining solutions  (Subscript[𝓊, 1]^σ), (Subscript[𝓊, 2]^σ), … to the homogeneous equation (2) with the same boundary values as the corresponding canonical solutions and vanishing stresses at infinity is outlined in standard classical treatises (e.g  [2]). These solutions are called harmonic solutions for the segment σ. For example, the harmonic solution (Subscript[𝓊, 3]^σ) has component fields \(u\_x = x\&^ . 𝓊\_3\%σ\) and \(u\_y = y\&^ . 𝓊\_3\%σ\)

with one dimensional plots confirming that (Subscript[u, x]) and (Subscript[u, y]) reduces to Cos[ϕ] and 0 at each eccentric angle ϕ on the perimeter of the elliptic hole. As the frequency increases, the fields make their transition from boundary values in the range -1 to 1 to zero at infinity over an increasingly smaller boundary layer.  The harmonic solution (Subscript[𝓊, 2 n + 7]^σ) has boundary displacements with component displacements

where (Subscript[u, x]) and (Subscript[u, y]) reduce to 0 and Sin[5 ϕ] at ϕ on the perimeter of the elliptic hole. Observe that a modest increase in the circular frequency results in a dramatic decrease in the transition region. Similar comments apply to the components of stress. Depending on the underlying precision, there is a circular frequency n beyond which harmonic and canonical harmonic solutions become indistinguishable.

Coefficients of the canonical harmonic solutions

Let 𝓃[1], 𝓃[2], …𝓃[S] be the number of harmonic solutions for segments 1, 2, ... S where S is the number of segments in the boundary. The key assumption is that each canonical harmonic solution (Subscript[u, i]^σ) can be represented as a linear combination of the combined set of harmonic solutions. That is

\(u\_i\%σ = \(C\_\(i, 1\)\%\(σ, 
              1\)\) 𝓊\_1\%1 + \(C\_\(i, 2\)\%\(σ, 
              1\)\) 𝓊\_2\%1 + ⋯ + \(C\_\(i, 𝓃[1]\)\%\(σ, 
              1\)\) 𝓊\_\(𝓃[1]\)\%1 + ⋯\)
+(Subscript[C, i, 1]^(σ, S)) (Subscript[𝓊, 1]^S) + (Subscript[C, i, 
            2]^(σ, 
            S)) (Subscript[𝓊, 2]^S) + ⋯ + (Subscript[C, i, 𝓃[
              S]]^(σ, 
            S)) (Subscript[𝓊, 𝓃[S]]^S)
\(\(\(=\)\(∑\+\(s = 1\)\%S\((\(C\_\(i, 1\)\%\(σ, 
                    s\)\) 𝓊\_1\%s + \(C\_\(i, 2\)\%\(σ, 
                    s\)\) 𝓊\_2\%s + ⋯ + \(C\_\(i, 𝓃[s]\)\%\(σ, 
                    s\)\) 𝓊\_\(𝓃[s]\)\%s)\) = 
∑\+\(s = 1\)\%S\((∑\+\(j = 1\)\%\(𝓃[s]\)\(C\_\(i, j\)\%\(σ, 
                      s\)\) 𝓊\_j\%s)\) \((5)\)\)\)\)

where (Subscript[C, i, j]^(σ, s)) is the coefficient of the j-th harmonic field  from segment s in the expansion of the i-th canonical harmonic solution from segment s.

Let v ↦ (Subscript[ℬ, i]^σ)[v] be the integral operator for fields v = v[ϕ] on segment σ with the orthogonality property \(ℬ\_i\%σ[𝓊\_j\%σ] = δ\_\(i, j\)\) (i.e 1 if i = j and 0 otherwise). For example, the ordering of fields in (3) and (4) give

\(ℬ\_\(n + 6\)\%σ[
        v] = \(-\(1\/π\)\) \(∫\_0\%\(\(-2\)\ π\)Cos[
              4  ϕ] 
            y\&^ . v[ϕ] dϕ\)\)

and

\(ℬ\_\(n + 6\)\%s[
        v] = \(-\(1\/π\)\) \(∫\_0\%\(\(-2\)\ π\)Sin[
              2  ϕ]\ 𝓃\&^[ϕ] . 
              v[ϕ] dϕ\)\)

in the case of segments σ and s with prescribed displacements and tractions respectively.  Note that (Subscript[ℬ, i]^σ)[
          (Subscript[u, j]^σ)] = (Subscript[δ, i, j]) because harmonic and canonical harmonic solutions have the same boundary values on their associated segments. Furthermore (Subscript[ℬ, i]^σ)[(Subscript[u, j]^s)] = 0 whenever σ ≠ s because the boundary values due to (Subscript[u, j]^s) vanish on segment σ. Then

\(ℬ\_h\%r[
        u\_i\%σ] = \(\(δ\_\(i, 
              j\)\) δ\_\(r, σ\) = ∑\+\(s = 1\)\%S\((∑\+\(j = 1\)\%\(𝓃[s]\)\(C\_\(i, j\)\%\(σ, 
                      s\)\) ℬ\_h\%r[𝓊\_j\%s])\)\ \ \ \ \ \ \((6)\)\)\)

which means that the array of canonical coefficients (Subscript[C, i, j]^(σ, s)) is the inverse of the array of Fourier coefficients (Subscript[ℬ, h]^r)[(Subscript[𝓊, j]^s)]. ElastostaticLab computes the Fourier coefficients and facilitates a graphical review of the Fourier coefficients {(Subscript[ℬ, h]^r)[(Subscript[𝓊, j]^s)] ∣ h = 1, 
        2, …} of the components of the boundary values due to the j-th harmonic solution from segment s over segment r.

An  Example

Let the number of harmonics associated with segments 1, 2, and 3 be 80, 60, and 80. The corresponding number of harmonic solutions are 322, 242, and 322. The associated number of Fourier coefficients is 886×886.  Let Mathematica's big number arithmetic correspond to 50 significant figures. On a 1.5GHz machine, ElastostaticLab takes 4 minutes (resp. 33 seconds) using high precision numbers (resp. machine numbers) to assemble the Fourier coefficients. The software includes graphical tools for reviewing the corresponding Fourier series. For example, the Fourier cosine coefficients of the boundary displacements (resp. tractions) on segment 1 (resp. 2) due to the harmonic solution with tractions \(Cos[3\ ϕ] 𝓉\&^\) on segment 2 are

The effects of various numbers of harmonic solutions on the Fourier coefficients can be observed and heuristically related to the accuracy of the canonical coefficients.

Mathematica's Inverse takes 80 minutes (resp. 46 seconds) using high precision numbers (resp.  machine numbers) to solve equation (6) for the canonical coefficients (Subscript[C, i, j]^(σ, s)) corresponding to the Fourier coefficients (Subscript[ℬ, h]^r)[(Subscript[𝓊, j]^s)].  The matrix associated with the Fourier coefficients has a determinant of 1.07.

The number of zeros in the maximum absolute error between the boundary values of the canonical harmonic solutions (Subscript[u, 1]^1), …, (Subscript[u, 322]^1), (Subscript[u, 323]^2), …, 
      (Subscript[u, 564]^2) and (Subscript[u, 565]^3), …, (Subscript[u, 886]^3) and the assigned boundary values for the displacements on segment 1 and tractions on 2 and 3 are

for big number arithmetic on the left and machine machine number arithmetic on the right. The presentation will include user interface details and  further numerical examples including solutions for the connected cracks that form an outer polygonal boundary.

References

[1] Wagner, D.B. Power Programming with Mathematica: The Kernel , McGraw Hill 1996
[2] Muskhelishvili, N.I., Some Basic Problems of the Mathematical Theory of Elasticity, translated by J.R.M. Radok, Groningen, P., Noordhoff, Ltd.,1963