Manipulate SIEs


This file calculates the applications in Y.-S. Chan, A. C. Fannjiang, and G. H. Paulino, Integral equations with hypersingular kernels -- theory and applications to fracture mechanics, Int. J. Eng. Sci., 41:683--720, 2003.


In [1]:
#Requires Gadfly
using ApproxFun, SingularIntegralEquations, Interact


Internal mode I crack in an infinite strip

For the crack opening displacement $\Delta v$: $$\frac{P.V.}{\pi}\int_{-1}^{+1}\frac{\Delta v(y)}{(y-x)^2}{\rm d}y + \frac{1}{\pi}\int_{-1}^{+1}K(x,y)\Delta v(y){\rm d}y = p(x), \qquad \Delta v(\pm1)=0,$$ where: $$K(x,y) = -\frac{1}{(x+y+2\epsilon)^2}+\frac{12(x+\epsilon)}{(x+y+2\epsilon)^3}-\frac{12(x+\epsilon)^2}{(x+y+2\epsilon)^4},\qquad p(x) = -1.$$


In [2]:
x = Fun(identity)
w = 1/sqrt(1-x^2)
d = domain(x)
d2 = d^2
B = dirichlet(d)
H2 = Hilbert(d,2)
S = DefiniteIntegral(d)
p = -Fun(one)

@manipulate for ϵ=1.05:0.01:1.5
    K = LowRankFun((x,y)->-1./(x+y+2ϵ).^2+12(x+ϵ)./(x+y+2ϵ).^3-12(x+ϵ).^2./(x+y+2ϵ).^4,d2)
    L = H2[w] + S[K*(w/π)]
    uSIE = [B,L]\[zeros(2),p]
    ApproxFun.plot(uSIE*w;axis=[0.0,1.5])
end


Out[2]:
x -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -1.50 -1.45 -1.40 -1.35 -1.30 -1.25 -1.20 -1.15 -1.10 -1.05 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30 2.35 2.40 2.45 2.50 2.55 2.60 2.65 2.70 2.75 2.80 2.85 2.90 2.95 3.00 -2 0 2 4 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 y

Mode III crack problem in nonhomogeneous medium

For the crack opening displacement $\Delta v$: $$G(x)\frac{P.V.}{\pi}\int_{-1}^{+1}K(x,y)\Delta v(y){\rm d}y = p(x), \qquad \Delta v(\pm1)=0,$$ where: $$K(x,y) = \frac{\beta e^{\beta(y-x)/2}}{2}\frac{K_1(\beta|y-x|)}{|y-x|},\qquad G(x) = e^{\beta x},\qquad p(x) = -1.$$


In [3]:
include("Fractureaux.jl")

d = Segment(-1.,1.)
sp = Space(d)
wsp = JacobiWeight(-.5,-.5,sp)
 = DefiniteLineIntegral(wsp)
x = Fun(identity,d)

B = dirichlet(d)
p = -Fun(one)

@manipulate for β = -2.0:0.1/π:2.0
    G = exp(β*x)
    K = ProductFun((x,y)->1/π*exp(β/2*(y-x))*β^2/2^2*(FK1(β*(y-x)/2)*(log(abs(β)/4) + γ) - GK1(β*(y-x)/2)/2),sp,wsp;method=:convolution)
    K0 = ProductFun((x,y)->exp(β/2*(y-x))*abs(β)/2*besseli(1,abs(β*(y-x)/2))/abs(y-x),CauchyWeight{0}(spwsp))
    K2 = ProductFun((x,y)->exp(β/2*(y-x)),CauchyWeight{2}(spwsp))

    L,p = [K2+K0+K],-Fun(one,sp)

    uSIE = [B,L]\[zeros(2),p/G]
    ApproxFun.plot(uSIE*w/G[1],axis=[0.0,12.0])
end


Out[3]:
x -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0 20.5 21.0 21.5 22.0 22.5 23.0 23.5 24.0 -20 0 20 40 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 y

Gradient elasticity applied to mode III cracks

For the crack opening displacement $\Delta v$: $$\frac{P.V.}{\pi}\int_{-1}^{+1}\left(\frac{-6\epsilon^2}{(y-x)^4} + \frac{1}{(y-x)^2}\right)\Delta v(y){\rm d}y = p(x), \qquad \Delta v(\pm1)=\Delta v'(\pm1)=0,\qquad {\rm where} \qquad p(x) = -1.$$


In [4]:
x = Fun(identity)
w = 1/sqrt(1-x^2)
d = domain(x)
d2 = d^2
B = [dirichlet(d),neumann(d)]
H2 = Hilbert(d,2)
H4 = Hilbert(d,4)
p = -Fun(one)
@manipulate for ϵ = 0.0005:0.0005:0.4
    L = -6ϵ^2*H4[w] + H2[w]
    uSIE = [B,L]\[zeros(4),p]
    ApproxFun.plot(uSIE*w;axis=[0.0,1.0])
end


Out[4]:
x -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 -1 0 1 2 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 y

The Faraday Cage

We seek that real function $\phi(z)$ that satisfies the Laplace equation: $$ -\Delta\phi = 0, $$ in the exterior domain, and the boundary condition: $$ \phi = \phi_0\qquad\mathrm{on~the~disks.} $$ such that: $$ \phi(z) = \log|z-z_s| + {\cal O}(1),\quad{\rm as}\quad z\to z_s, $$

$$ \phi(z) = \log|z| + o(1),\quad{\rm as}\quad z\to\infty. $$

with a source of $2\pi$ units at $z_s = 2$.

S. J. Chapman and D. P. Hewett and L. N. Trefethen, Mathematics of the Faraday cage, SIAM Rev., 2015.


In [5]:
require("Gadfly")
z_0 = 2.0
ui(x,y) = logabs(complex(x,y)-z_0)
g1(x,y) = 0.5

@manipulate for N in [4,6,8,10,12], r in [1e-1,5e-2,1e-2,5e-3,1e-3], domains in [:Circles,:Segments,:Alternating]
    cr = exp(im*2π*[0:N-1]/N)
    crl,crr = (1-2im*r)cr,(1+2im*r)cr

    if domains == :Circles
        dom = (Circle,cr,ones(length(cr))r)
    elseif domains == :Segments
        dom = (Segment,crl,crr)
    elseif domains == :Alternating
        dom = (Segment,crl[1:2:end],crr[1:2:end])  (Circle,cr[2:2:end],ones(length(cr[2:2:end]))r)
    end

    sp = Space(dom)
    cwsp = CauchyWeight(spsp,0)
    uiΓ, = Fun(t->ui(real(t),imag(t))+0im,sp),DefiniteLineIntegral(dom)

    @time G = GreensFun(g1,cwsp;method=:Cholesky)

    @time φ0,un=vec([0 ;1 [G]]\Any[0.,uiΓ])

    us(x,y) = -logkernel(un,complex(x,y))/2
    ut(x,y) = real(ui(x,y) + us(x,y))

    pl = Main.Gadfly.plot(ut,-2.,3.,-1.5,1.5,Main.Gadfly.Stat.contour(samples=75, levels=[-2:0.5:-1.0,-0.7:0.2:-0.3,-0.2:0.1:1.2]))
    pl.layers = [ApproxFun.layer(dom),pl.layers]
    Main.Gadfly.render(pl)
end


elapsed time: 0.226397618 seconds (49222992 bytes allocated, 26.80% gc time)
elapsed time: 0.165398966 seconds (41899136 bytes allocated, 40.52% gc time)
Out[5]:
x -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 -10 -5 0 5 10 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 1 0 2 -2 -1 f(x,y) -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 -4.5 -4.4 -4.3 -4.2 -4.1 -4.0 -3.9 -3.8 -3.7 -3.6 -3.5 -3.4 -3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 -5 0 5 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 y

In [ ]: