Displacement Field and Elastic Energy of a Circular Twist Disclination for Large Deformations - an Example how to Treat Nonlinear Boundary Value Problems with Computer Algebra Systems

11 Pages·0.3 MB·English

by  Alexander Unzicker

Displacement Field and Elastic Energy of a Circular Twist Disclination
for Large Deformations - an Example how to Treat Nonlinear Boundary
Value Problems with Computer Algebra Systems
3
0
0 Alexander Unzicker
2
Pestalozzi-Gymnasium
n
Melssheimerstr. 11, D-81247 Mu¨nchen, Germany
a
J [email protected]

8
2 Karl Fabian

] Department of Geosciences, University of Bremen
i
c Box 330440, D-28334 Bremen, Germany
s
- [email protected]
l
r
t February 2, 2008
m
.
t
a
m Abstract
-
d
A circular twist disclination is a nontrivial example of a defect in an elastic continuum that causes large de-
n
formations. The minimal potential energy and the corresponding displacement field is calculated by solving the
o
c Euler-Lagrange-equations.The nonlinear incompressibility constraint is rigorously taken into account.
[ By using an appropriate curvilinear coordinate system a finer resolution in the regions of large deformations is

obtained and the dimension of the arising nonlinear PDE’s is reduced to two. The extensive algebraic calculations
1
v that arise are done by a computer algebra system (CAS). The PDE’s are then solvedby a difference scheme using
1 the Newton-Raphson algorithm of successive approximations for multidimensional equations. Additional features
3 for global convergence are implemented. To obtain basic states that are sufficiently close to the solution, a one
5
dimensionallinearizedversionoftheequationissolvedwithanumericalcomputationthatreproducestheanalytical
1
results of Huang and Mura (1970).
0
3 With this method, rigorous solutions of the nonlinear equations without any additional simplifications can be
0 obtained. The numerical results show a contraction of the singularity line which corresponds to the well-known
/ Poynting effect in nonlinear elasticity. This combination of analytical and numerical computations proves to be a
t
a versatile method to solve nonlinear boundary value problems in complicated geometries.
m
-
d
n 1 Introduction Bullough, and Smith (1955) that dislocation density is
o
related to the differential geometric concept of torsion.
c Topologicaldefectsinelasticmediausuallygeneratelarge After the result of Frank (1949)1 the elastic energy of
:
v deformations which require finite elasticity for their de- topological defects in elastic continua has been consid-
Xi scription.Thetheoryofnonlinearelasticitygoesbackto ered in various publications on dislocations and discli-
Cauchy (1827), and significant contributions came from nations (Eshelby 1949; Kroupa 1960; Huang and Mura
r
a Love(1927),Signorini(1930),andfortheincompressible 1970; Nabarro 1979). Puntigam (1996) discussed topo-
case,fromRivlin(1948a).Foranintroductiontononlin- logical defects in a field theoretic context.
ear elasticity, the reader may consider Beatty (1987) or However,calculationsoftheelasticenergyoftopolog-
thestandardtextbooksTruesdellandToupin(1960)and icaldefects have alwaysbeen restrictedto the linearized
Truesdell and Noll (1965). equations.Inthis case the problemcanbe treated using
Topological defects, namely dislocations and discli- tensor analysis and stress function tensors of first and
nations, have been investigated by means of continuum second degree (Kr¨oner 1959). Two assumptions usually
theories (Kr¨oner 1959; Kr¨oner 1960) that were devel-
1He showed that the energy increases relativistically with the
oped after the discovery of Kondo (1952) and Bilby,
dislocationvelocity.
1
enter these approximations: lartwistdisclination.Theappropriatecoordinatesystem
Firstly, one restricts to linear elasticity in the sense and the respective transformation of the the displace-
that shearing stresses are assumed to cause a state of ment vector is also explained there. All the numerical
simple shear in the material (as in Huang and Mura issuesareaddressedin section3,andthe results aredis-
1970). cussed in section 4.
Secondly, the energy is assumed to be a function of
the distorsion tensor that can be approximated by the
2 Analytical description
gradient of the displacement field. Then, by applying
Stoke’stheorem,anintegrationovertheboundariesonly
2.1 Basic concepts of nonlinear contin-
can be performed. Regarding the first point, however,it
uum mechanics
can be shown (e.g. Rivlin 1948c, p. 467), that in the
general case shearing stresses alone cannot maintain a
Nonlinear elasticity was founded in 1827 by Cauchy.
stateofsimpleshearinthematerial.Alsothesecondas-
Oneassumesanundeformed,euclidiccontinuumwith
sumptionisjustifiedforsmall deformationsonly(Rivlin
Cartesian coordinates X = (X,Y,Z) (the so-called ‘ref-
1948c, p. 476). In the region close to the defect core,
erence configuration’) and attaches in every point a dis-
where the deformations are large the linear approxima-
placement vector u that points to the coordinates x =
tion fails to predict finite energies. Thus, these results
(x,y,z) of the deformedstate (‘configuration’): u=x−
are reasonable only outside the core regionin where the
X.
energy density diverges.
From u one deduces a quantity that transforms the
Largedeformationsincrystalshavebeeninvestigated
coordinatesX of the undeformed state to those x of the
rarely (Frank 1951), and if, not by means of analytical
deformed state: the deformation gradient
but by generaltopologicalmethods (Rogula 1976)with-
out considering the elastic energy. ∂x
F:= ,
There are several proposals how to treat nonlinear ∂X
effectsinelasticmediawithdefects.Teodosiu(1982)ob-
or, in components,
tained second-order effects for a straight screw disloca-
tion by applying Willis’ (1967;1970)scheme which goes 1+u u u
X Y Z
 
back toSignorini(1930).This requires,however,certain F:= v 1+v v , (1)
X Y Z
physical assumptions for every additional order of ap-  wX vY 1+wZ 
proximation.Guo(2001)modellednumericallythelarge
deformations of a hyperelastic membrane. The cylindri- where(u,v,w)denotethecomponentsofuandthe sub-
cal symmetry allowed a the reduction to a onedimen- scripts differentiation. The symmetrical tensor
sional PDE.
T
B:=FF , (2)
Lazar(2002a;2002b)obtainedsolutionsforedgeand
screw dislocations in an elastoplastic theory.
iscalled(right)Cauchy-Greentensor.Itisconvenientto
Povstenko (1995; 2000) has treated the twist discli-
introducethe so-calledprincipal invariants (I1,I2,I3)of
nation with a nonlocal modulus. Even if his method in-
a tensor that are defined as follows:
volvesnumericalintegration,itisbasedonananalytical
derivation of the stresses, and does not allow a contrac- I1 =λ1+λ2+λ3 (trace) (3)
tion of the singularity line. I2 =λ1λ2+λ2λ3+λ3λ1 (4)
Here a moregeneralapproachis proposed:The total
I3 =λ1λ2λ3 (det) (5)
elastic energy, i.e. the integral over the energy density
mustbe aminimum undervariationofthe displacement (the λ are the eigenvalues of B). Then, in general, the
i
field u. Additional constraints - like in the present case energydensityW isafunctionoftheprincipalinvariants
incompressibility - are included by means of a Lagrange
multiplier. If the energy depends on first derivatives of W =W(I1,I2,I3) (6)
u only,this leads to Euler-Lagrangeequations ofsecond
oftheCauchytensorB(e.g.Beatty1987),andthestress
order,evenifthemethodgivenhereallowsthetreatment
tensor T (defined as traction per surface element) is
of higher order equations.
given by the constitutive equation
The general method outlined above is applied to a
circulartwist disclination(Huang and Mura 1970)in an −1
T=β1B+β0E+β−1B (7)
incompressible hyperelastic continuum. It is a numeri-
cally automatized process for obtaining rigorous solu- with the βi(I1,I2,I3) being functions of the principal
tions ofthe nonlinearequationswithoutspecialphysical invariants. For small deformations, i.e. for I1’s close to
assumptions. 1, and I2.I3 close to 3, the β’s have a fixed values - that
Section2 gives a briefintroduction to nonlinear con- can be related to the known elastic constants in linear
tinuummechanics,followedbyadescriptionofthecircu- elasticity.
2
2.2 Incompressibility TheEuler-Lagrangeequationsarethusobtainedfrom
The nonlinear condition of incompressibility is given by ∂L d ∂L d ∂L d ∂L
− − − =0 (16)
∂u dx∂u dy∂u dz∂u
I3 =λ1λ2λ3 =det B=(det F)2 =1. (8) ∂L d ∂Lx d ∂Ly d ∂Lz
− − − =0 (17)
∂v dx∂v dy∂v dz∂v
andnot,asmanytextsonlinearelasticitystate,divu=0. x y z
As a consequence, the elastic energy W depends on the ∂L d ∂L d ∂L d ∂L
− − − =0 (18)
principal invariants I1 and I2 of B only and eqn. 6 re- ∂w dx∂wx dy∂wy dz∂wz
duces to
and
W =W(I1,I2) (9) ∂(x+u,y+v,z+w)
−1=0, (19)
In the present paper, the special case ∂(x,y,z)
which is equivalent to (8). Using the special form of W
W(I1)=C(I1−3) (10) from (11) without a Lagrange multiplier still leads to a
is considered, which is called a Neo-Hookean material2 linear PDE. Adopting W from (11) is not an essential
simplification of the problem, since the condition (19)
where C is a constant.
is the main nonlinearity. The principles of the following
calculations apply in the same way to a nonlinear W.
2.3 Strain-Energy function expressed by
displacements
2.5 Circular twist disclination
Considering again a Cartesian coordinate system and
A circular twist disclination is a nontrivial example of
taking into account (1), (3) and (10), the energy per
a topological defect which enforces large deformations
unit volume for an incompressible, neo-Hookean mate-
of an elastic continuum. The twist disclination can be
rial reads
realized as follows (see fig. 1):
W =C(ǫ +ǫ +ǫ ), (11)
xx yy zz
3
where C corresponds to the shear modulus µ, and the
respective components of strain are given by (Cauchy
1827; Love 1927, p. 60 and Rivlin 1948b, p. 461)
R
1 z
2 2 2
ǫ =u + (u +v +w ) (12)
xx x 2 x x x r
1
2 2 2
ǫ =v + (u +y +w ) (13)
yy y 2 y x y
1
2 2 2
ǫ =w + (u +z +w ). (14)
zz z 2 z x z
Hereby u=(u,v,w) denotes the displacement vector in
x-, y-, and z-directions and u , u ... etc. the respective
x y
partial derivatives.
2.4 Euler-Lagrange equations
The energy of the topological defect can be found by
minimizing W under the variation of the displacement Figure 1: Schematic description how to produce a twist
field u which fulfill the correct spacial boundary condi- disclinationinanelasticcontinuum.Thesolidtorus(size
tions. In addition, the nonlinear constraint (8) is taken R,thicknessr)isremoved.Thenthematerialiscutalong
into account by means of a Lagrange multiplier λ. This the hatched surface. After twisting the cut faces by the
requires the minimization of the Lagrangian amountoftheFrankangleΩ,thematerialisrejoinedby
glueing.
L=W +λ(det F−1). (15)
2This corresponds to the condition βi = const. Note that in Imagine RI3 filled with elastic material and remove
finiteelasticitytheenergydensitydependsuponthefunctionsβi. a solid torus centered around the origin and with z as
Therefore,nogeneral‘canonical’energydensityexists. symmetryaxis(fig.1).Cutthematerialalongthesurface
3Inanincompressiblematerial,3C equalsYoung’smodulusE.
boundedbytheinnercirclewithradiusR−rofthetorus
in the x−y-plane (hatched disk in fig. 1). The cut faces
aretwistedwithrespecttoeachotherbytheFrankangle
Ω, and glued together again.
3
Ifoneshrinkstheremovedtorustoasingularityline,
in the limit r → 0 one obtains a circular twist disclina-
tion with radius R. In the subsequent computations, a
sequence of finite values of r is used.
Theabovedescriptionofthetwistdisclinationisequiv-
alent to the one given by Huang and Mura (1970), who
used however the term ’edge disclination’4 to indicate
thattheFrankvectorisperpendiculartothedisclination
line, in analogy to edge dislocations where the Burgers
vector is perpendicular to the defect line. They investi-
gatedanalyticallythetwistdisclinationinlinearapprox-
imation.
2.6 Transformation to curvilinear coor-
dinate systems
The boundary conditions for a specific problem can be
formulated most conveniently in an appropriately cho- Figure 2: Toroidal coordinate system obtained by ro-
sen curvilinear coordinate system. Furthermore, differ- tating a bipolar coordinate system around the z-axis
entscalefactorsinthiscoordinatesystemsallowtopave (η = 0). For symmetry reasons, only the region η > 0
the regions of interest - for example regions of large de- and 0 < ϑ < π (at φ = 0) is shown. The value η = ∞
formations - more densely with coordinate lines, which correspondstoafocus ofthebipolarsystemor,afterro-
is of greatnumerical advantage.Frequently,one can use tating aroundthe z-axis,to a circular singularity line in
symmentries of a problem and reduce it to a lower di- the toroidalsystem.η =ϑ=0 isthe infinitely farpoint.
mension. The variety of orthogonal coordinate systems
allows to design an appropriate system for nearly ev-
In the region in the vicinity of the circular singular-
ery problem, even for those with less symmetry. To give
ity the Jacobian determinant is very small. η = ϑ = 0
some examples, straight line defects may bemodelled in
corresponds to the infinitely far point. The components
a cylindrical system ( e.g. using a log-transformed ra-
of the displacement vector in η-, ϑ-, and φ-direction are
dial component), closedline defects with various elliptic
denoted as u, v and w.
coordinate systems.
The transformation of differential operators like div
Fortheabovedescribedtwistdisclination,thetoroidal
andcurl incurvilinearsystemsislongknown(Love1927,
system, as visualized in Fig. 2, is a suitable choice. To
p. 54) and is easily done by CAS5.
obtain the 3D-toroidal coordinate system from the 2-D
To transform the expressions for the energy W and
bipolar system which is actually shown in Fig. 2, one
for the condition (19) in curvilinear coordinates, how-
must rotate the latter around the vertical axis.
ever, this is not sufficient, since the displacement vector
The transformationis given by
components (u,v,w) are involved. One has to take into
cos(ϕ) sinh(η) accountthatthebasisvectorsinacurvilinearsystemare
x= (20)
−cos(ϑ)+cosh(η) subjecttochangeanddifferentiatethem.Thistechnique
is well-known in differential geometry since it is needed
sin(ϕ) sinh(η)
y = (21) foracurvedspacetimeinthegeneraltheoryofrelativity.
−cos(ϑ)+cosh(η)
The spatial derivatives of the displacement vector
z = −cos(ϑsi)n+(ϑc)osh(η) (22) tura=nsfuoirm=at(iuon,vo,fwt)hearbeas∂∂ixsukv6e.ctToarks,intgheinutsouaalccdoeurinvtattihvee
has to be replaced by the ‘covariant’ derivative:
Theparameteroftheazimuthalrotationisφ∈[0,2π[,
D d
whereas the polar angle ϑ ranges from 0 to π. The hy- ui = ui−Γ iuj (23)
Dxk dxk jk
perbolic coordinateη ranges from −∞ (left focus) to ∞
(right focus). The connection Γ i is calculated from the Jacobianma-
jk
4Toavoidconfusionwithscrewandedgedislocations,theterms trix Aji (A11 = ∂∂xη etc.) by
‘twist’ and ‘wedge’ are commonly used for disclinations (deWit
1973a; 1973b; Zubov 1997). Rogula(1976) callsthe circulartwist Γ i =B i ∂ A m, (24)
disclination a ‘third type of defect’ and Unzicker (1996; 2000) a jk k ∂xm j
‘screw dislocation loop’. In any case, the twist disclination re- 5Here the Mathematica(cid:13)c package ‘VectorAnalysis’ has been
ferred here causes locally a Volterra distortion of the 5th kind
used.
(seePuntigam1996). 6Incontrasttoeqns.(16-18),(u,v,w)and ∂ referstoη-,ϑ-,
∂xk
andφ-directionsinthetoroidalsystem.
4
B j istheinverseofA i(Schouten1954,p.121ff.).Eqn.(24)π,η = η )= u(ϑ= ϑ ,η =η ) implicitly imple-
i j max min max
has to be applied to all derivatives occurring in (16-19). ment this condition. Only for infinitesimal small values
This again is done by a CAS. ofΩtheucomponentwouldvanish.Inthelinearapprox-
Thetransformationconsiderablycomplicatestheequa- imation, u, v and w vanish at ϑ = 0. At the symmetry
tions. For example (11) becomes in toroidal coordinates axis η = 0 the condition of incompressibility causes the
vanishingof uandw without anexplicitsetting to zero.
W =u (coshη−cosϑ)+v (coshη−cosϑ)+ (25)
η ϑ
u(csch η−cosϑcoth η)−2vsinϑ−usinhη+
3 Numerical methods
2 2 2 2
w (coshη−cosϑ) +w (−cosϑ+coshη) +
(cid:0) η ϑ
w2(cosϑcoshη−1)2csch η2+w2sinϑ2+ 3.1 Discretization with a
(v (coshη−cosϑ)+usinϑ)2+ difference scheme
η
2
(uη(coshη−cosϑ)−vsinϑ) + The nonlinear Euler-Lagrange equations generated by
2
(u(csch η−(cosϑcothη))−vsinϑ) + the CAS are discretized on a twodimensional lattice by
2 evaluating the coefficients η,sinϑ etc. at every lattice
(v (coshη−cosϑ)−usinhη) +
ϑ
2 point. The displacements and their derivatives are still
(u (coshη−cosϑ)+vsinhη) /2
ϑ (cid:1) maintained in analytic form. One obtains 4 variables
(displacement vector u and Lagrange multiplier λ) at
Thetransformationofeqns.(16-19),i.e.thefullnon-
linearPDEtosolve,isnotwrittenouthere,7theyamount each of the n grid points. Therefore 4n nonlinear equa-
tions have to be solved simultaneously.
toatextfileofabout30kB.Theextensivealgebraiccal-
culations necessary to transform the equations are done
by a CAS. 3.2 Newton’smethod inmultidimensions
The solutions of the arising nonlinear system are ob-
2.7 Boundary conditions
tainedbytheNewton-Raphsonmethodofsuccessiveap-
proximations(Press1993,chap.9.7),whichisamultidi-
The circular singularity of the twist disclination corre-
mensionalextensionofNewton’s rootfinding algorithm.
sponds to η = ∞. A toroidal core η > η sourround-
max
As in one dimension the algorithm starts from an ba-
ing the singularityη =∞ is removedfrom the material.
sicstateandusesfirstderivativestoiterativelycalculate
The chosen values of η = 2...3.625 correspond to
max
approximations of the solution.
core radii r = 0.23...0.049. Since η =ϑ=0 is the in-
The numerical computation of the 4n derivatives of
finitely far point, only a region η > η and ϑ > ϑ
min min
theEuler-Lagrangeequationsishowevernumericallyin-
is considered.
hibitive.Thisproblemissolvedbyanalytically differenti-
The distance R of the singularity line to the symme-
atingtheEuler-Lagrangeequationsateverylatticepoint
try axis η =0 (the size of the defect) is set to 1.
with respect to the displacements u,v and w and their
The cylindrical symmetry of the problem in the co-
derivatives using a CAS. Thus the 18 quantities λ , λ ,
ordinate ϕ reduces the 3D-problem to two dimensional η ϑ
u ,u ,v ,v ,w ,w ,u ,v ,w ,u ,v ,w ,w
one.Therefore,u,v andw dependonη andϑonly.Fur- η ϑ η ϑ η ϑ ηϑ ηϑ ηϑ ϑϑ ηη ηη ϑϑ
enteringtheEuler-Lagrangeequationsarediscretizedby
thermore, due to mirror symmetry ϑ is restricted to the
difference molecules (Bronstein 1985, p. 767f.).
region 0<ϑ<π.
Onlythesimplestdifferencemoleculesareimplemented
Theboundaryvaluesofu,vandwareleftfreewhere
to avoid instability of the Newton-Raphson method due
ever possible in orderto allow relaxing to the configura-
tounrealisticsmoothnessassumptions.Thepossiblegrid
tionofminimalenergy.Thisisdoneatthesurfaceofthe
dimensions do not allow to assume smoothness of the
torus sourrounding the singularity at η = η (torus
max
solution because the distorsions in the vicinity of the
shown in fig. 1), and at η =0. The boundary conditions
toroidalcoreregionusuallybecomeverylargeandhigher
at ϑ = 0 are enforced by the symmetry of the problem
orderdifferenceoperatorscontainnonext-neighborcou-
and at ϑ = π by the ’cut-and-glue’ condition with the
pling.
Frank angle Ω.
Inaddition,the differenceschemeallowsforaconve-
Thus at the disk defined by ϑ=π, v vanishes and u
nient formulation of the boundary conditions.
and w are fixed to a purely circular displacement corre-
sponding to the Frank angle Ω (see fig. 2). During en-
ergy minimization, the shape of the removed torus is 3.3 Global convergence
kept fixed, otherwise the material could overlap, which
Although the above implementation is numerically very
is physically impossible. The boundary values w(ϑ =
efficient,theglobalconvergenceofthedirectmultidimen-
7NotethatminmizationofW intheneo-Hookeancaseeqn.11 sional Newton-Raphson method is still not guaranteed.
stillleadstoalinearPDE.Onlytheconstraintofincompressibility
To obtain the global convergence of the algorithm one
(8)enforcesthesubstantialnonlinearityofthefinalPDE.
5
usually defines a functional V on the solution space as sons, only the component w of the diplacement is var-
the square of the l.h. sides of eqns. 16-19. If the New- ied. The components u and v are adjusted to meet the
ton step increases V with respect to the previous state, aboverequirement.Thisprocedurerequiresto solveonly
the algorithm used here ‘walks back’ along the Newton a onedimensional linear PDE for calculating the basic
direction looking for a one-dimensional local minimum states.
of V. (cfr. Press 1993, 9.7). The existence of such a lo-
cal minimum is guaranteed, since at the basic state, V 3.5 Interpolation and extrapolation
bydefinitionofthe gradientdecreasesalongthe Newton
direction. It turns out that convergence improves if the Forfinergridsandsmallercoreradii,itbecomesincreas-
‘walkback’alreadyisundertakenifVdecreasedslightly, ingly difficult to find solutions by means of the above
and the Newton step was ‘accepted’ only if V decreases method.
bylessthanafactorof10.Afactor10iseasilyreachedin It turns out that grid refinement by interpolating a
the vicinity of the solution, where the Newton-Raphson previouslyfoundsolutiononacoarsergridisbyfarmore
algorithm converges quadratically with distance. If the efficientthanstarting withanindependent new solution
Newton step is accepted, its result is used as the basic of the linearized system. The interpolation is performed
stateofthenextiterationsteparoundwhichthe PDEis bya2D-splinealgorithm(Press1993,chap.3.6)andthe
againlinearized.In contrastto the advicegiveninPress interpolatedfunction – evaluatedatthe refined grid– is
(1993), it is found here that it is useful to minimize V takenasnew basicstate.Startingfromthis state,anew
along the Newton direcetion precisely.Although this re- solution usually is obtained by Newton-Raphson. How-
quiresmorefunctionevaluations,itperformssuperiorto ever, this refined solution still has the same core radius
taking some premature value for a new linearization. as the previous coarse solution. For extending the grid
Fortheone-dimensionalminimizationalongtheNew- towardssmaller core radii,that is to largerηmax, a very
ton direction, the golden sectio algorithm (Press 1993, efficient method is to attach one lattice line to the old
chap. 10.1) is used. Where the function is sufficiently solutionbylinearextrapolation,andtousethisextended
flat8 (functionals V at three points vary by less than a solutionagainasanewbasicstate.Allsolutionsforvery
factor1.2),thegoldensectioisaccompaniedbyparabolic small core radii have been obtained in this way.
interpolation (Press 1993, chap. 10.2) of the minimum.
Thiscombination,isfoundtoperformquitewellinfind- 3.6 Matrix inversion and function eval-
ing solutions. Moreover,it is not sensitive to small vari- uation
ations of the above parameters.
Matrix size increases very rapidly with the number of
grid points. For example, the four functions (λ,u,v,w)
3.4 Basic states
on a 30×30 lattice lead to a 3600×3600 matrix with
The success of the Newton-Raphson method critically almost 13 × 106 coefficients, which has to be inverted
dependsonthebasicstateusedforthefirstlinearization. at every Newton step. Since the matrix is sparse, the
Inthemultidimensionalcasethereislittle hopethatthe inversion is more efficiently done by specialized sparse
algorithmdirectly convergesto the solution without the matrix algorithms. Here the corresponding packages of
above mentioned global-convergencemethods. But even MatLab(cid:13)c are applied over a data exchange interface
then the algorithm ends up in a local minimum if the with Mathematica(cid:13)c.
first basic state is too far from the solution. There is SinceCASareslowinevaluatingtrigonometricfunc-
for example, no chance to find a solution when starting tions, the evaluation is done by an external C routine.
with a nonspecialized basic state, e.g. with all variables Thisspeedsupthecomputationofthematrixcoefficients
settozero.Sometimes linearapproximations,e.g.states by a factor 10 as compared to the CAS.
obtained by assuming small deformations,are chosenas With this features, satisfactory results are obtained
basicstates.Inthepresentcase,thelinearizedconstraint by running the program on a usual PC.
of incompressibility div u= 0, still is too far from real-
ity to be used for constructing a first basic state. Thus,
4 Results
the incompressibilty condition is directly implemented
by using the symmetries of the toroidal coordinate sys-
4.1 Linear approximation
tem. In first approximation, it can be assumed that the
displacement takes place along concentric circles with In linear approximation the twist disclination can be
respect to the symmetry axis. Only for small displace- treatedanalytically.HuangandMura(1970)havefound
ments this corresponds to the tangential displacement that the total elastic energy is
of the linearized equation. Thus, no dilatation at all is
1 r r
allowedin the circularapproximation.For practicalrea- W =( µ Ω2R3){[2+(1− )2]K(1− )− (26)
3 R R
8Andifthepointinthemiddlewasnottooexcentric. r r
2
2[1+(1− ) ]E(1− )},
R R
6
W
10
8
6
4
2
Om
p p 3p p
4 2 4
Figure 3: Comparision of the numerical results for the
linear approximation (filled circles) with the analytical
results (solid line) of Huang and Mura (1970) in case of
a core radius of 0.09 (η = 3, η = ϑ = 0.05,
max min min
R = 1, µ = 3). The defect energy increases with the
amount of the Frank angle Ω (ranging from 0 to π). For
small Ω, the results agree.
where K and E denote the complete elliptic integral of
the first and second kind. Eqn. 26 can be approximated
by
1 1 4
2 3
W = µ Ω R ( log(8R/r)− ). (27)
3 2 3
This result is used to test the algorithm for calculating
thebasicstates(section3.4).Fig.3showsthatnumerical
and analytical results agree for small Ω. The predicted
Ω2-dependence in eqn. 26 is reproduced by the linear
algorithm.
For larger values of Ω, however, the energies calcu-
lated numerically lie below the parabola of Huang and
Mura (1970). This illustrates the two aspects of a linear
approximation:bothintheanalyticalandinthe numer-
ical treatment shear stresses are assumed to cause pure Figure 4: Solution of the nonlinear PDE. The respec-
shear,butthelatterapproachliftsthesmall-deformation tive undistorted state is shown in fig. 2, with singu-
assumption by taking the energy function (11). larity as dottet line. A very large deformation with a
Frankangle Ω=π is implemented by twisting the lower
4.2 Solutions of the complete nonlinear boundary ϑ = π by the amount of π/2. The lower part
PDE −π/2 < ϑ < 0 (not shown here) is twisted in the oppo-
sitedirection.The coreradiusisr about0.23,η =2,
max
The linear assumption of shear stresses causing shear η =ϑ =0.05.
min min
deformationsonlyisliftednow.Thesolutionsofthefully
nonlinear equation presented here are constraint to the
the case Ω = π where very large deformations occur.
Fig.4visualizesthesedeformationsforacoreradiusr =
0.23. The other parameters are µ = 3, R = 1,η =
min
ϑ = 0.05. The dottet line indicates the singularity.
min
Fig. 2 shows the corresponding undeformed state.
The components of the displacementvector u, v and
w, for the same parameters as in fig. 4, are shown in
fig. 5. Note that the grid shown in fig. 5 is equidistant
in the toroidal coordinates η and ϑ, but not in a carte-
sian frame, as it is evident from figs. 2 and 4. There,
7
Figure6:Distributionoftheelasticenergyforancircular
twist disclination with η = 2.5 (equivalent to r =
max
0.145), and η =ϑ =0.05.
min min
the lattice spacing is much finer in the region of large
deformations.
Even if it demands quite a lot of imagination, the
reader should verify that the displacement components
in fig. 5, attachedto fig. 2,yield indeed the deformation
fig. 4.
In fig. 6 (r = 0.145) the distribution of the elastic
energy plotted. At each grid point the energy density,
weightedby the correspondingtoroidalvolume element,
is shown. Most of the energy is stored in the region
nearby the cut surface, where its density is greater than
inthe vicinity of the core regionη =η .This effectis
max
even more pronounced for smaller core radiii r.
4.3 Elastic energy
Fig. 7 shows the total elastic energy of a circular twist
disclination in function of η which parametrizes the
max
core radius.
Using the global minimization techniques discussed
above, solutions over a wide range from η = 2.0
max
(which corresponds to r = 0.23) to η = 3.625 ( or
max
r =0.049)arefound.Fig.7showsthedependenceofthe
energy(filled circles)fromη.9 It is considerablyreduced
incomparisontothelinearapproximationofHuangand
Mura (1970) (solid line). Although from numerical re-
sults like those of fig. 7 it cannot be deduced that the
energy approachesa finite value in the limit r →0, fig.6
suggeststheenergydensitydecreasessufficientlyforcon-
vergence.
The total elastic energies shown in fig. 7 are calcu-
Figure 5:vector components u(η,ϑ), v(η,ϑ) andw(η,ϑ)
lated from a lattice of 16 points in ϑ- direction, whereas
for the solution fig. 4 (core radius r=0.23).
the number of points in η-direction ranged from 17 to
9η standsforηmax+ηmin here.
8
W ElasticEnergy W
14
12
10
8
6
4
2
eta
2.25 2.5 2.75 3 3.25 3.5
Figure 7: Dependence of the total elastic energy of cir-
cular edge disclinations for η = 2.0 (r = 0.23) to
max
η =3.625(r =0.049).The filled circles are obtained
max
from the nonlinear modelling, the solid line shows the
analytic resultof Huang andMura (1970).The cutoff of
the infinitely far point was η = ϑ = 0.05, which
min min
corresponds to a distance of about 20 from the origin.
There were 16 lattice points in ϑ-direction and a range Figure 8: Example of a small core radius r = 0.09. Vec-
from 17 to 30 for η. tor component v in ϑ-direction for a solution of the lin-
earized equations.
30. The outer limit of the modelled region is defined by
η = ϑ = 0.05. Other solutions with a different
min min
number of grid points, and smaller values for η and
min
ϑ (0.01, 0.02) showed only a very slight dependence
min
on these parameters.
4.4 Poynting effect
The solutionsofthenonlinearcaseshowacharacteristic
phenomenon which only occurs in the elasticity theory
of finite deformations, the so-called Poynting effect.
The Poynting effect in its classic form is (Truesdell
and Noll 1965, p. 193): ‘When an incompressible cylin-
der,freeonitsoutersurface,istwisted,itexperiencesan
elongation ultimately proportional to the square of the
twist.’
Toillustratethedifferencebetweenthelinearapprox-
imationandthenonlinearsolution,fig.8and9showthe
vector component v in ϑ-direction for r = 0.09. In the
regionclosetoη =0,i.e.inthevicinityofthesymmetry
axis,vhasnegativevaluesandthedisplacementvectoris
pointing upwardsin the regionwherez >0. (cfr.fig. 2).
By mirror symmetry it points downwards in the region
ofz <0 (−π <ϑ<0).As aconsequencethe materialis
stretched and incompressibility forces the central region
Figure 9: As fig.8, but for the nonlinear case. In the re-
close to the cut surface to contract.
gion nearby η = 0, i.e. in the vicinity of the symmetry
The elongationindicated for z >0 is observedinthe
axis, v has negative values. that means it is pointing
nonlinear case only. It is a displacement normal to the
upwards. The material is stretched and the cut surface
planeonwhichthestressisinduced(thecutsurface).In
contracts.
thelinearcasefromtheassumptionofpureshearcaused
v =0 at the symmetry axis (fig.8).
The Poynting contraction apparently approaches a
9
5 Conclusions
R Poynting Contraction
0.9
By combining analytical and numerical approaches, a
0.8 generalmethodforobtainingrigouroussolutionsofnon-
linear boundary value problems with nontrivial geome-
0.7 tries has been developed.
With this new method, for the first time the fully
0.6 nonlinear problem of a circular twist disclination in an
hyperelastic imcompressible material could be treated.
0.5 Whiletheextensivecalculationsnecessaryforformu-
lating Euler-Lagrange-equationsin a curvilinear coordi-
eta nate systemwere done by the computer algebrasystem,
2.25 2.5 2.75 3 3.25 3.5
thesolutionofthediscretizedequationswasobtainedby
theNewton-Raphsonmethodwithadditionalglobalcon-
Figure 10: Poynting contraction as fraction of the orig-
vergencefeatures.Theresultsgiveanonlinearcorrection
inal size of the defect. The contractions seems to reach
to the total elastic energy obtained in analytical treat-
a lower bound at .58, at a core radius of 0.049, which
ments (Huang and Mura 1970). Moreover, the method
corresponds to η =3.625.
max allowed for the first time to model the contraction of
the singularity line and to observe the Poynting effect,
positivelowerlimitforverysmallcoreradiir.Theresults whichischaracteristicforforlargedeformationsinfinite
showninfig.10suggestthatthecontractiondoesnotfall elasticity. An application of the outlined method to to
below 58% of the original radius. othernonlinearPDE’swith boundaryvalue problemsin
several dimensions seems possible.
4.5 Problems
Acknowledgement. We are grateful to Dr. Markus
Forbothηmax <2.0andηmax >3.75itbecomesincreas- Lazar for hints to the literature.
inglydifficulttofindsolutionsofthenonlinearequations.
In the first case the difficulty apears to be related to
References
the imposedboundary conditionswhichfix the torusdi-
ameter r (i.e. the core radius). In the second case the
Beatty,M.F.(1987).Topicsinfiniteelasticity:Hyper-
Poynting contraction leads to a net change in the torus
elasticity of rubber, elastomers, and biologicaltis-
volume, which apparently leads to numerical problems.
sues - with examples. Appl. Mech. Rev. 40, 1699–
On the other hand, for small core radii,the program
1733.
finds solutions which satisfy the numerical convergence
criterium, but show a physically unreasonable zigzag- Bilby, B. A., R. Bullough, and E. Smith (1955). Con-
patternin the displacements.Possiblythe occur inhuge tinous distributions ofdislocations:anew applica-
deformations (a volume element of the size of the core tion of the methods of non-riemannian geometry.
radiusbecomesstretchedtothelengthofthesingularity Proc. Roy. Soc. London, Ser. A 231, 263–273.
line) become numerically untractable at least with the
Bronstein, I. N. (1985). Taschenbuch der Mathematik
present method (‘relatively’ simple coordinate system).
(22 ed.). Leipzig: Teubner.
However,there still may be a deeper physical reasonfor
Cauchy (1827). De la pression ou tension dans un
this instability.
corps solide; sur la condensation et la dilatation
A simpler analogue to the circular edge disclination
des corps solides; sur les equations qui expriment
in case of large deformations is a twisted rubber band.
lesconditionsd’equilibreoulesloisdemouvement
Everyday experience shows that while the Poynting ef-
interieurd’un corpssolide.Exercises de mathema-
fect only slightly elongates the band, the band sponta-
tique.
neously looses its axial symmetry and coils up as soon
as the twist surpasses a critical value. deWit, R. (1973a). Theory of disclinations: II. conti-
In case of the edge disclination, the present model – nous and discrete disclinations in anisotropic elas-
whichessentiallydepends onrotationalandmirrorsym- ticity.J.Res.Nat.Bur.Standards 77A(1),49–100.
metry is not able to correctly model such a transition.
deWit, R. (1973b).Theory of disclinations: III. conti-
It is therefore not astonishing that meaningful solutions
nous and discrete disclinations in isotropicelastic-
are found up to a maximal η only.
max ity. J. Res. Nat. Bur. Standards 77A(3), 359–368.
Eshelby, J. D. (1949).Uniformly moving dislocations.
Proceedings of the Physical Society London A 62,
307–314.
10