Computing forces on interface elements exerted by dislocations in an elastically anisotropic crystalline material

3.9 MB·

by  Bing Liu

Computing forces on interface elements exerted by
dislocations in an elastically anisotropic crystalline
material
6
1
B. Liu1,*, A. Arsenlis1, S. Aubry1
0
2

b
1
e LawrenceLivermoreNationalLaboratory,Livermore,CA94550,USA
F *Correspondingauthor: bingliu a llnl.gov(B.Liu)
(cid:13)
0
2

February 23, 2016
]
i
c
s
-
l
r Abstract
t
m
.
t Drivenbythegrowinginterestinnumericalsimulationsofdislocation–interface
a
m interactionsingeneralcrystallinematerialswithelasticanisotropy,wedevelopalgo-
- rithmsfortheintegrationofinterfacetractionsneededto coupledislocationdynam-
d
ics with a finite element or boundary element solver. The dislocation stress fields
n
o in elastically anisotropic media are made analytically accessible through the spher-
c
ical harmonics expansion of the derivative of Green’s function, and analytical ex-
[

pressionsfortheforcesoninterfaceelementsarederivedbyanalyticallyintegrating
2
v the spherical harmonics series recursively. Compared with numerical integration by
8 Gaussianquadrature,thenewlydevelopedanalyticalalgorithmforinterfacetraction
3
integrationishighlybeneficialintermsofbothcomputationprecisionandspeed.
9
7
0
Keywords: Dislocation dynamics; anisotropic elasticity; finite domain; interface
.
1
tractionintegration;analyticalsolution
0
6
1
: 1 Introduction
v
i
X
r It is well known that a dislocation approaching a free surface experiences an attractive force,
a
so called image force [1]. In multiphase materials or single-phase elastically anisotropic poly-
crystalline materials, a dislocation near a phase or grain boundary is subject to a similar force
due to the change in elasticity across the interface, which can be either attractive or repulsive
[2–4]. As the magnitudes of these forces are inversely proportional to the dislocation–interface
1
distance, the effect of free surfaces on dislocation motion and multiplication becomes signifi-
cantinsubmicrometre-sizedcrystals[5–8],andtheelasticinteractionsbetweendislocationsand
phase/grainboundarieshaveastrongerimpactonthemechanicalpropertiesofthenanostructured
materials[9–12]thanthoseoftheircoarse-grainedcounterparts.
For free surfaces and phase/grain boundaries of elastically anisotropic half spaces, such vir-
tual forces (negative energy gradients) can be determined through the image force theorem of
Barnett and Lothe [13–17]. For finite domains of more complex shapes, this problem of elastic
interactionsbetweendislocationsandinterfacescanonlybesolvedthroughthecouplingofadis-
locationdynamics(DD)codeandafiniteelement(FE)orboundaryelement(BE)solver[18–28].
Insuchcases,thestressfieldofthedislocationinaninfiniteelasticallyhomogeneousmediumis
usedtocalculatethetractiononthesurfaceboundingtheelasticsolid,andthencorrectionfields
areaddedtoimposetheproperboundaryconditionsonthedomain. Forfreesurfacesthisresults
inanimpositionofanequalandoppositesurfacetractionsuchthatthenetresultisazerotraction
condition,andforphase/grainboundariesthereisatractionbalanceanddisplacementcontinuity
conditionthatmustbeimposed. ThisworkisfocusedonthefirstpartofDD-FE/BEsimulations
of dislocation–interface interactions, i.e. determination of forces on interface elements due to
tractionsimposedbydislocationsinaninfiniteelasticallyanisotropicmedium.
Interface tractions exerted by dislocation stress fields must be integrated into nodal forces in
aFEorBE solver. Thesenodalforces aresurfaceintegralsofthetractionfieldT (forceperunit
area)overtheindividualinterfaceFEorBEelements[20,24,28],
(cid:90) (cid:90)
F(n) = T (x)N(n)(x)dS = [σ(x) n]N(n)(x)dS, (1)
S · S
S S
where F(n) is the force on a FE or BE node n of an interface element S exerted by dislocation
stress field σ(x), n is the interface normal, and N(n)(x) are the FE or BE shape functions. In
s
this paper, we will refer to these nodal forces as traction forces. The traction force on a FE or
BEnodeisanalogoustotheinteractionforceduetodislocationstressfieldonadislocationnode
in nodal based (one-dimensional FE) dislocation dynamics models [20, 29–31]. The dislocation
interactionforcesarelineintegralsofthePeach-KoehlerforcefPK (forceperunitlength)along
theindividualdislocationsegments[20,31],
(cid:90) (cid:90)
f(n) = fPK (x)N(n)(x)dL = [σ(x) b t]N(n)(x)dL, (2)
L · × L
L L
wheref(n) istheforceonanodenofadislocationsegmentLexertedbydislocationstressfield
σ(x), with b, t, and N(n)(x) being the Burgers vector, line direction, and shape function of the
L
dislocation segment, respectively. While the dislocation interaction force calculation is the core
ofaDDmodel,theinterfacetractionforceevaluationisthemainlinkinDD-FE/BEfinitedomain
simulations. However, the compromise between accuracy and efficiency of the numerical inte-
grations of dislocation interaction forces and surface traction forces has been the bottleneck for
large-scale DD simulations [31] and DD–FE/BE simulations of the elastic interactions between
dislocations and free surfaces [20, 24, 32–34]. Due to the limitations of numerical integrations,
analternativeanalyticalintegrationofinterfacetractionforcesishighlydesirable.
2
For isotropic elastic media, the stress field of a dislocation loop can be calculated analyt-
ically through line integrations of the derivatives of Green’s function along piecewise straight
dislocation segments [35, 36]. Arsenlis et al. [31] developed analytical expressions for dislo-
cation interaction force calculations that involves double integrations along the individual pairs
of interacting dislocation segments, i.e. the first line integration along the source segment to
get its stress field and the second line integration along the receiving segment to obtain the in-
teraction force. These analytical expressions for dislocation interaction force calculations have
already been used in many DD simulations, e.g. discovery of ternary dislocation junctions in
body-centered cubic metals [37], interpretation of the size-dependent strength for micrometer-
sized crystals [38], determination of low-angle grain boundary penetration resistances [39, 40],
observationofstrainlocalizationviadefect-freechannelsinhighlyirradiatedmaterials[41],and
revealing the mechanisms of dynamic recovery during high temperature creep of single-crystal
superalloys [42]. Queyreau et al. [43] have recently formulated analytical expressions to calcu-
late surface traction forces induced by stress field of a dislocation in isotropic elastic media for
rectangularsurfaceelements,whichgivetheprecisesolutionstotripleintegrals,i.e. oneintegral
along the dislocation segments for the stress field and then a double integral over the surface
elementtoobtainthetractionforce.
Inanisotropicelasticitytheoryofdislocations,theanalyticalexpressionforthestressfieldof
anarbitrarydislocationloopdoesnotexist[44,45]. Thedislocationstressfieldcalculationbased
on Stroh’s formulism [35, 46] requires numerical integrations [47, 48] or solving an eigenvalue
problem for a six by six matrix [48] to obtain the associated matrices and their angular deriva-
tives. The dislocation stress field computations using Mura’s formula [49] has relied on direct
numerical integrationsof the derivatives ofGreen’s function [50]. Recently, Aubry andArsenlis
[51] used a truncated spherical harmonics expansion to approximate the derivatives of Green’s
function, and have analytically integrated the spherical harmonics series to calculate dislocation
stress field (single integrals) and interaction force (double integrals) for straight dislocation seg-
ments. Inthiswork,weusethesamesphericalharmonicsexpansionofthederivativesofGreen’s
functionformulatedinthepreviouswork[51],andanalyticallyintegratetheassociatedspherical
harmonics series to determine the interface traction force (triple integrals) in anisotropic elastic
mediaforquadrilateralsurfaceelements.
2 Method
In a homogeneous infinite linear elastic solid, the stress field of a dislocation loop can be ex-
pressedintermsofacontourintegralalongtheloop[49],
(cid:73)
∂G
σ (x) = (cid:15) C C b vp (x x)dx , (3)
js ngr jsvg pdwn w ∂x − (cid:48) (cid:48)r
L d
where C is the elastic stiffness tensor, (cid:15) is the permutation tensor, b is the Burgers vector of
ijkl
thedislocationloop,and∂G /∂x isthederivativeoftheGreen’sfunctionG (x x),which
vp d vp (cid:48)
−
is defined as the displacement in the x -direction at point x in response to a unit point force in
v
thex -directionappliedatpointx.
p (cid:48)
3
e
3
T
eˆ
x
ψ
ξ
eˆ
e y e
1 2
Figure1: Unitsphereinanisotropicelasticity
The Green’s function in an anisotropic elastic medium has been obtained as a single integral
[44],
1 (cid:90) π
G = M 1(ξ)dψ, (4)
vp 4π2R v−p
0
where
(cid:15) (cid:15) (ξξ) (ξξ)
M 1(ξ) = vsm prw sr mw ,
v−p 2(cid:15) (ξξ) (ξξ) (ξξ)
lgn 1l 2g 3n
with the notation (ξξ) = ξ C ξ . R is the norm of vector R = x x, i.e. R = R . T is
ij k kijl l − (cid:48) (cid:107) (cid:107)
the direction of vector R, i.e. T = R/R. ξ is a unit vector that varies in the plane ξ T = 0, eˆ
x
·
and eˆ are two orthogonal unit vectors in the same plane, and the angle between eˆ and ξ is ψ,
y x
Fig.1.
ThecorrespondingintegralexpressionforthederivativeofGreen’sfunctionis[52],
∂Gvp = 1 (cid:90) π(cid:0) T M 1 +ξ N (cid:1)dψ, (5)
∂x 4π2R2 − d v−p d vp
d 0
where N = C M 1M 1(ξ T +ξ T ). See the overview of Bacon et al. [44] for more
vp jrnw v−j n−p r w w r
details.
The derivative of the Green’s function is a product of a part depending only on 1/R2 and an
angularpartg dependingonlyonthedirectionT
(cid:90) π
(cid:0) (cid:1)
g (T) = g (θ,φ) = T M 1 +ξ N dψ, (6)
vpd vpd − d v−p d vp
0
where(θ,φ)arethesphericalcoordinatesofT.
4
There is no analytical expression for g , but the function g (T) is suitable for decompo-
vpd vpd
sitioninsphericalharmonics.
l
(cid:88)∞ (cid:88)
g(T) = glmYm(T) (7)
l
l=0 m= l
−
Theexpansioncoefficientsglm areindependentofT (θ,φ),andaredefinedas
(cid:90) 2π(cid:90) π
glm = g(θ,φ)Ym (θ,φ)sinθdθdφ (8)
l ∗
0 0
ThesphericalharmonicsYm aredefinedasthecomplexfunctions
l
(cid:115)
2l+1(l m)!
Ym(θ,φ) = − Pm(cosθ)eimφ, (9)
l 4π (l+m)! l
where Pm are the associated Legendre polynomials. To be consistent with the definition of
l
elastic stiffness tensor C , we rewrite Ym in the Cartesian coordinate system (e ,e ,e ), i.e.
ijkl l 1 2 3
x = T e ,y = T e ,andz = T e ,intheformof
1 2 3
· · ·
[(l m)/2]
−(cid:88)| |
Ylm(x,y,z) = fm(x,y) Q¯l|m|(k)zl−|m|−2k (10)
k=0
where
(cid:26)
(x+iy)m m 0
f (x,y) = ≥
m (x iy) m m < 0
−
−
and
(cid:115)
(cid:18) (cid:19)(cid:18) (cid:19)(cid:18) (cid:19)
( 1)m+k m! 2l+1(l m)! l 2l 2k l 2k
Q¯m(k) = − − − −
l 4π2 2l 4π (l+m)! k l m
Thefunctiong canthenbedefinedas
l [(l m)/2]
(cid:88)∞ (cid:88) (cid:0) (cid:1) −(cid:88)
g(x,y,z) = (x+iy)mglm Qm(k)zl m 2k (11)
(cid:60) l − −
l=0 m=0 k=0
where we note (x) is the real part of x, Q0(k) = Q¯0(k) when m = 0, and Qm(k) = 2Q¯m(k)
(cid:60) l l l l
whenm > 0.
Usingthesphericalharmonicsseriesexpansionofitsangularpartanddefininge = e +ie
12 1 2
withx+iy = R e andz = R e ,thederivativeoftheGreen’sfunctioncanbeevaluatedby
R · 12 R · 3
∂Gvp(R) = (cid:88)∞ (cid:88)l [(l−(cid:88)m)/2] (cid:18)Qm(k)glm (R·e12)m(R·e3)l−m−2k(cid:19). (12)
∂x (cid:60) l vpd Rl 2k+2
d −
l=0 m=0 k=0
5
This definition involves a quotient of terms that are a function of R, which depends only on two
variablesmandl 2k. Itcanbesimplifiedfurthertoobtain
−
2q+1 (cid:18) (cid:19)
∂Gvp(R) = (cid:88)∞ (cid:88) Sqm(R·e12)m(R·e3)2q+1−m (13)
∂x (cid:60) vpd R2q+3
d
q=0 m=0
where Sqm is a sum of products composed of Qm(k) and glm. The Green’s function and its
vpd l vpd
derivatives depend only on odd powers of 1/R. This property means that in the spherical har-
monicsexpansion,thenon-zerotermscorrespondtotheoddpowersof1/R.
Within linear elasticity, the total stress field of a dislocation network is a superposition of
the stress fields of individual dislocation segments in the network. For a straight dislocation
segment linking two nodes at Cartesian coordinates x and x , respectively, its stress field in an
1 2
anisotropicelasticmediumcanthenbeexpressedby
(cid:32) (cid:33)
σ (x) = (cid:15) C C b (cid:88)∞ 2(cid:88)q+1 Sqm (cid:90) x2 (R·e12)m(R·e3)2q+1−mdx (14)
js ngr jsvg pdwn w (cid:60) vpd R2q+3 (cid:48)r
q=0 m=0 x1
ConsideraquadrilateralelementdelimitedbyitsCartesiancoordinatesx ,x ,x ,andx as
3 4 5 6
showninFig.2,wherexisthecoordinatethatspanswithintheelement,andx isthecoordinate
(cid:48)
thatspansalongthedislocationsegmentdelimitedbyx andx . ThevectorR = x x canbe
1 2 (cid:48)
−
rewrittenasR = yt+rp+sq,witht = x2 x1 ,p = x4 x3 ,andq = x5 x3 .
x2−x1 x4−x3 x5−x3
Using n = p q and dS = p q(cid:107)dr−ds,(cid:107)with a(cid:107)lin−ear(cid:107)shape func(cid:107)tio−n in(cid:107) the form of a
×
p q (cid:107) × (cid:107)
(cid:107) × (cid:107)
four-termpolynomialN(n)(x) = a +a r+a s+a rs,thetractionforceF(n) canbeexpressed
0 1 2 3
as
(cid:90) s2(cid:90) r2
F(n) = σ(x) (p q)(a +a r+a s+a rs)drds (15)
0 1 2 3
· ×
s1 r1
Combining Eq. (14) and Eq. (15) with dx = tdy, the interface nodal force due to disloca-
(cid:48)
−
tiontractioninananisotropicelasticmediumcanbedeterminedthrough
Fj(n) = −tr(cid:15)sabpaqb(cid:15)ngrCjsvgCpdwnbw(cid:80)∞q=0(cid:80)2mq=+01
(16)
(cid:16) (cid:17)
(cid:60) Svqpmd(cid:82)ss12(cid:82)rr12(cid:82)yy12 (R·e12)mR(R2q·+e33)2q+1−m(a0 +a1r+a2s+a3rs)dydrds
TheunknownpartsofEq.(16)arethefollowingseriesofintegrals:
(cid:90) s2(cid:90) r2(cid:90) y2 (R e )i(R e )j
12 3
K = · · dydrds
ijk Rk
s1 r1 y1
(cid:90) s2(cid:90) r2(cid:90) y2 (R e )i(R e )j
Kr = · 12 · 3 rdydrds
ijk Rk
s1 r1 y1
(cid:90) s2(cid:90) r2(cid:90) y2 (R e )i(R e )j
Ks = · 12 · 3 sdydrds
ijk Rk
s1 r1 y1
(cid:90) s2(cid:90) r2(cid:90) y2 (R e )i(R e )j
Krs = · 12 · 3 rsdydrds,
ijk Rk
s1 r1 y1
6
p
x
4
x
s 3
q
R x r
x y
1
x
x 6
t 0
x
5
x
2
b
Figure2: Geometryandassociatedvariablesforinterfacetractionintegration
wherei+j = k 2,andthecommonpartoftheintegrandswillhereafterbereferredasI ,i.e.
ijk
−
(R e )i(R e )j
12 3
I = · · .
ijk Rk
The main effort of this work lies in solving these series of triple integrals. We find that with
the analytical solutions of a few seed integrals, all other integrals in the spherical harmonics
expansion can be calculated analytically through recurrence relations. In Fig. 3, the required
tripleintegralsforthetractionforcecalculationareindicatedwithabluecolor,andtheresttriple
integrals are necessary for the recursive integration to reach the required triple integrals of the
same expansion order q and the next expansion order q + 1. As illustrated by the blue dashed
lines in Fig. 3, the double integrals at the intermediate levels are needed for the calculations of
thetripleintegralsatthehighestlevels,andthesingleintegralsatthelowestlevelsmustbeused
tocalculatethedoubleintegralsattheintermediatelevels.
Therequiredrecurrencerelationsareconstructedusingthefollowingpartialderivatives:
∂ I = iαI +jφI k(R t)I (17)
y ijk (i 1)jk i(j 1)k ij(k 2)
− − − · −
∂ I = iβI +jθI k(R p)I (18)
r ijk (i 1)jk i(j 1)k ij(k 2)
− − − · −
∂ I = iγI +jψI k(R q)I (19)
s ijk (i 1)jk i(j 1)k ij(k 2)
− − − · −
∂ (yI ) = I +y∂ I (20)
y ijk ijk y ijk
∂ (rI ) = I +r∂ I (21)
r ijk ijk r ijk
7
Ky Kr Ks
Kijk K055 ijk ijk ijk
K K Kr
045 145 045
K K K Kr Kr
035 135 235 q = 1 035 135
K K K Kr Kr Kr
125 225 325 025 125 225
K K K Kr Kr Kr
215 315 415 115 215 315
K K K Kr Kr Kr
305 405 505 205 305 405
K
033
K K Kr
023 123 q = 0 023
K K K Kr Kr
013 113 213 013 113
K K K Kr Kr Kr
103 203 303 003 103 203
K
011 q = 1
K K − Kr
001 101 001
Hy Hr Fy Fs Er Es
Hijk Fijk Eijk H043 ijk ijk ijk ijk ijk ijk
H H Hy
033 133 q = 0 033
H H H Hy Hy
023 123 223 023 123
H H H Hy Hy Hy
113 213 313 013 113 213
H H H Hy Hy Hy
203 303 403 103 203 303
H
021
H H q = 1 Hy
011 111 − 011
H H H Hy Hy
001 101 201 001 101
H00 1 q = 2 Jy Jr Js
− − ijk ijk ijk
Jy
031
Jy Jy
021 121 q = 1
Jy Jy Jy −
011 111 211
Jy Jy Jy
101 201 301
Jy
Jy 01−J1y q = −2
00 1 10 1
− −
Figure3: Interfacetractionintegrationthroughhierarchicalrecurrencerelations
8
∂ (sI ) = I +s∂ I (22)
s ijk ijk s ijk
∂ (R tI ) = I +R t∂ I (23)
y ijk ijk y ijk
· ·
∂ (R pI ) = I +R p∂ I (24)
r ijk ijk r ijk
· ·
∂ (R qI ) = I +R q∂ I (25)
s ijk ijk s ijk
· ·
whereα = t e ,β = p e ,γ = q e ,φ = t e ,θ = p e ,andψ = q e .
12 12 12 3 3 3
· · · · · ·
ThefirstthreerecurrencerelationsEq.(26),Eq.(27),andEq.(28)forthetripleintegralsare
obtained by successive integrations over y, r, and s of the partial derivatives Eq. (17), Eq. (18),
andEq.(19),respectively. AnothertworecurrencerelationsEq.(29)andEq.(30)canbeverified
bythetripleintegraldefinitions.
Ky +cKr +dKs = 1 (cid:2)iαK +jφK E (cid:3) (26)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
cKy +Kr +fKs = 1 (cid:2)iβK +jθK F (cid:3) (27)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
dKy +fKr +Ks = 1 (cid:2)iγK +jψK H (cid:3) (28)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
K = αKy +βKr +γKs (29)
(i+1)jk ijk ijk ijk
K = φKy +θKr +ψKs , (30)
i(j+1)k ijk ijk ijk
where c = p t, d = q t, f = p q, Ky = (cid:82)s2(cid:82)r2(cid:82)y2I ydydrds, E = (cid:82)s2(cid:82)r2I drds,
F = (cid:82)s2(cid:82)y·2I dyds·, and H ·= (cid:82)r2ij(cid:82)ky2I s1dyrd1r.yT1heijdkouble integraijlks H s,1Fr1, aijnkd E
ijk s1 y1 ijk ijk r1 y1 ijk ijk ijk ijk
needtobepreviouslycalculatedusingthesecondsetofrecurrencerelationsgivenbelow.
The first two recurrence relations Eq. (31) and Eq. (32) for the double integrals H are
ijk
obtained by successive integrations over y and r of the partial derivatives Eq. (17) and Eq. (18),
respectively. The third recurrence relation Eq. (33) for the double integrals H is obtained by
ijk
successive integrations over y and r of the partial derivatives Eq. (20) and Eq. (21), and then
summation of these two integration equations. Another two recurrence relations Eq. (34) and
Eq. (35) can be verified by the double integral definitions. Analogously, the recurrence relations
for the double integrals F Eq. (36) to Eq. (40) and for the double integrals E Eq. (41) to
ijk ijk
Eq.(45)areobtainedusingthecorrespondingpartialderivativesanddoubleintegraldefinitions.
Hy +cHr +dsH = 1 (cid:2)iαH +jφH Jr (cid:3) (31)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
cHy +Hr +fsH = 1 (cid:2)iβH +jθH Jy (cid:3) (32)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
dsHy +fsHr +s2H = 1 (cid:2)iγsH +jψsH
ij(k+2) ij(k+2) +ij((kk+2) i kj 2)H(i−1)j+k rJy +iy(jJ−r1)k(cid:3) (33)
− − − ijk ijk ijk
H = αHy +βHr +γsH (34)
(i+1)jk ijk ijk ijk
H = φHy +θHr +ψsH (35)
i(j+1)k ijk ijk ijk
Fy +dFs +crF = 1 (cid:2)iαF +jφF Js (cid:3) (36)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
dFy +Fs +frF = 1 (cid:2)iγF +jψF Jy (cid:3) (37)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
9
crFy +frFs +r2F = 1 (cid:2)iβrF +jθrF
ij(k+2) ij(k+2) +ij((kk+2)i kj 2)F(i−1)+jkyJs +is(jJ−y1)k(cid:3) (38)
− − − ijk ijk ijk
F = αFy +βrF +γFs (39)
(i+1)jk ijk ijk ijk
F = φFy +θrF +ψFs (40)
i(j+1)k ijk ijk ijk
1 (cid:2) (cid:3)
Er +fEs +cyE = iβE +jθE Js (41)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
1 (cid:2) (cid:3)
fEr +Es +dyE = iγE +jψE Jr (42)
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k − ijk
(cid:2)
cyEr +dyEr +y2E = 1 iαyE +jφyE
ij(k+2) ij(k+2) ij(k+2) k (i−1)jk i(j−1)k(cid:3) (43)
+(k i j 2)E +sJr +rJs
− − − ijk ijk ijk
E = αyE +βEr +γEs (44)
(i+1)jk ijk ijk ijk
E = φyE +θEr +ψEs , (45)
i(j+1)k ijk ijk ijk
where Hy = (cid:82)r2(cid:82)y2I ydydr, Hr = (cid:82)r2(cid:82)y2I rdydr, Fy = (cid:82)s2(cid:82)y2I ydyds, Fs =
(cid:82)s2(cid:82)y2I ijksdydsr,1 Eyr1 ij=k (cid:82)s2(cid:82)r2Iijk rdrdrs1, Ey1s ijk= (cid:82)s2(cid:82)r2ijIk sdrs1ds,y1Jyijk = (cid:82)y2Iijkdy,
Jsr1 y=1 (cid:82)ijrk2I dr, ainjdk Js s=1 (cid:82)r1s2Iijk ds. Theijskingle ins1tegrr1alsijJky , Jr ,iajknd Js y1muijskt be
ijk r1 ijk ijk s1 ijk ijk ijk ijk
formerlycalculatedusingthethirdsetofrecurrencerelationsgivenbelow.
The first two recurrence relations Eq. (46) and Eq. (47) for the single integrals Jy are ob-
ijk
tained by integration over y of the partial derivative Eq. (17) and multiplying both sides of the
integrationequationbyαandφ,respectively. Thethirdrecurrencerelationforthesingleintegrals
Jy Eq. (48) is obtained by integration over y of the partial derivative Eq. (23). Analogously,
ijk
the recurrence relations for the single integrals Jr (from Eq. (49) to Eq. (51)) and Js (from
ijk ijk
Eq.(52)toEq.(54))areobtainedusingthecorrespondingpartialderivatives.
α (cid:104) (cid:105)
Jy = iαJy +jφJy I +Jy [(β αc)r+(γ αd)s] (46)
(i+1)j(k+2) k (i 1)jk i(j 1)k − ijk ij(k+2) − −
− −
φ (cid:104) (cid:105)
Jy = iαJy +jφJy I +Jy [(θ φc)r+(ψ φd)s] (47)
i(j+1)(k+2) k (i 1)jk i(j 1)k − ijk ij(k+2) − −
− −
(cid:104)
Jy = 1 R tI +i[(β αc)r+(γ αd)s]Jy
ij(k+2) k[R2−(R·t)2] · ijk − − (i−1)jk(cid:105) (48)
+j[(θ φc)r+(ψ φd)s]Jy +(k 1 i j)Jy
− − i(j 1)k − − − ijk
−
β (cid:2) (cid:3)
Jr = iβJr +jθJr I +Jr [(α βc)y +(γ βf)s] (49)
(i+1)j(k+2) k (i 1)jk i(j 1)k − ijk ij(k+2) − −
− −
θ (cid:2) (cid:3)
Jr = iβJr +jθJr I +Jr [(φ θc)y +(ψ θf)s] (50)
i(j+1)(k+2) k (i 1)jk i(j 1)k − ijk ij(k+2) − −
− −
(cid:104)
Jr = 1 R pI +i[(α βc)y +(γ βf)s]Jr
ij(k+2) k[R2−(R·p)2] · ijk − − (i−1)jk(cid:105) (51)
+j[(φ θc)y +(ψ θf)s]Jr +(k 1 i j)Jr
− − i(j 1)k − − − ijk
−
γ (cid:2) (cid:3)
Js = iγJs +jψJs I +Js [(α γd)y +(β γf)r] (52)
(i+1)j(k+2) k (i 1)jk i(j 1)k − ijk ij(k+2) − −
− −
ψ (cid:2) (cid:3)
Js = iγJs +jψJs I +Js [(φ ψd)y +(θ ψf)r] (53)
i(j+1)(k+2) k (i 1)jk i(j 1)k − ijk ij(k+2) − −
− −
10