On the development of a triangular, multi-field user-element for Abaqus Dortmund, Technische ...

49 Pages·2010·2.27 MB·English

Dortmund, Technische Universit¨at
Institute of Mechanics
Dortmund, Germany
April 16, 2010
On the development of a triangular,
multi-field user-element for Abaqus
Author: M. van Dijk, 0625286
Supervisors: Dipl. Ing. C. Hortig
dr. ir. R.H.J. Peerlings
Document number: MT 10.13
Contents
1 Introduction 3
2 Weak form and linearization 4
2.1 Momentum balance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Heat balance equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Nonlocal balance equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Linearization of balance momentum . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.5 Linearization of heat balance equation and nonlocal balance equation. . . . . . . . 8
2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3 Finite element framework 12
3.1 Finite element approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2 Assembling of global system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Algorithmic tangent for thermo-elasticity . . . . . . . . . . . . . . . . . . . . . . . 17
4 Abaqus User-element 20
4.1 UEL Template . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2 Input file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.3 Obtaining Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5 Results 23
5.1 Test environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.2 Thermo-elastic element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.3 Thermo-viscoplastic element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
6 Conclusions and Recommendations 27
A Special Tensor operations 29
B Matlab Script 31
C UEL Template 38
D Input file 41
E Convention DOFs in ABAQUS 42
F Python script 43
Chapter 1
Introduction
In order to predict the structural behavior of products, such as the response on thermal and
mechanical loads, Finite Element Method (FEM) programs are used more and more frequently. A
wide variety of elements and material formulations have been developed to ensure reliable results.
Beyond the formulations implemented in commercial codes, advanced users desire to implement
their own elements and routines to solve non standard problems. Some of the commercial FE
codes, such as Abaqus, offer user interfaces for that purpose.
Due to the flexibility, the Abaqus subroutine User-element (UEL) will be employed in this
work. With this subroutine a maximum number of 20 additional degrees of freedom (DOF’s)1
becomes possible. The main purpose of a user-designed element is to provide the stiffness matrix
as well as the residual vector, as needed in a context of solving the nonlinear system of equations
using Newton-Rhapson. The User-element can be written in C-code or in Fortran. In this report
the User-element will partly be written in Matlab and then compiled to C-code, as developing
in Matlab is easier compared to C-code. In the following, we will derive a general algorithmic
background, allowing additional degrees of freedom such as temperature, damage etc. To keep
things simple, the final implementation will be restricted to a linear plane strain thermo-elastic
element. Thus, we work with the two standard DOF’s, displacement u , and one additional
1,2
degree of freedom, the temperature θ.
Themainpartofthisreportfocusesontheframeworktocalculatethestiffnessmatrixandthe
right hand side vector. In chapter 2 we start with deriving the weak form of momentum balance,
the balance of heat equation and the nonlocal balance equation. As will be shown, both, the
balanceofheatequationandthenonlocalbalanceequationcanbeexpressedwithacorresponding
general form. The next step is to express the derived equations in terms of a finite element
approximation. This will be discussed in chapter 3. Here, we also discuss the implementation in
Matlab. Because most of the formulas obtained are straight forward, they will not be discussed in
detail. In chapter 4 we discuss the implementation of the User-element into the user subroutine
UEL. We start with presenting the template file of such a UEL, as well as the input file, which
gives the model information to the Abaqus solver. Further, we will describe how to obtain the
results, calculated with the User-element. We end up with a demonstration of simulation results,
obtained with the developed User-elements. Some of these results will be compared with Abaqus
standard elements.
1Intotalthereare30DOF’savailable,butthefirst10cannotbeuserdefined. ForfurtherdetailsseeChapter4
Chapter 2
Weak form and linearization
In order to solve the linear system of equations, the stiffness matrix and right hand side vector are
needed. The first step, in order these quantities, is to calculation the weak form of the different
balance laws and to linearized them.
2.1 Momentum balance
The momentum balance is given by
∇(cid:126) ·σσσ+ρ(cid:126)b=ρ(cid:126)x¨ (2.1)
Here,σσσ is the stress tensor,(cid:126)b is the body force, and (cid:126)x¨ is the spatial acceleration. For deriving the
weak form of this equation the weighted residual form will be used. The weighted residual form
of equation (2.1) can be calculated by taking the scalar product with a vector-valued test function
δ(cid:126)x, followed by integration over the domain V
(cid:90) (cid:16) (cid:17)
∇(cid:126) ·σσσ+ρ(cid:126)b−ρ(cid:126)x¨ ·δ(cid:126)xdV =0 (2.2)
V
(cid:90) (cid:16) (cid:17) (cid:90) (cid:16) (cid:17)
∇(cid:126) ·σσσ ·δ(cid:126)xdV − ρ(cid:126)x¨−ρ(cid:126)b ·δ(cid:126)xdV =0 (2.3)
V V
By application of the chain rule, it can be shown that
(cid:16) (cid:17)
∇(cid:126) ·(σσσ·δ(cid:126)x)= ∇(cid:126) ·σσσ ·δ(cid:126)x+σσσ :∇(cid:126)δ(cid:126)x (2.4)
Using this rule in (2.3) gives
(cid:90) (cid:90) (cid:16) (cid:17)
{∇(cid:126) ·(σσσ·δ(cid:126)x)−σσσ :∇(cid:126)δ(cid:126)x}dV − ρ(cid:126)x¨−ρ(cid:126)b ·δ(cid:126)xdV =0 (2.5)
V V
Additionally, the divergence theorem says
(cid:90) (cid:90)
∇(cid:126) ·(cid:126)adV = (cid:126)n·(cid:126)adS (2.6)
V S
where (cid:126)a is a smooth vector function on the domain V and (cid:126)n is the unit outward normal to the
boundary S of the domain V. Using the divergence theorem in (2.5) gives
(cid:90) (cid:90) (cid:90)
(cid:126)n·(σσσ·δ(cid:126)x) dS− σσσ :∇(cid:126)δ(cid:126)xdV + (ρ(cid:126)b−ρ(cid:126)x¨)·δ(cid:126)xdV =0 (2.7)
S V V
Furthermore, the traction force(cid:126)t is defined as
(cid:126)t=(cid:126)n·σσσ (2.8)
2.2 Heat balance equation 5
Using this relation in (2.7) will give
(cid:90) (cid:90) (cid:90)
σσσ :∇(cid:126)δ(cid:126)xdV + (ρ(cid:126)x¨−ρ(cid:126)b)·δ(cid:126)xdV − (cid:126)t·δ(cid:126)xdS =0 (2.9)
V V S
2.2 Heat balance equation
The strong form of the heat balance equation is given by
ρcθ˙−κ ∇(cid:126) ·(∇(cid:126)θ)=r (2.10)
0
Here, ρ is the density, c the heat capacity at constant deformation, θ˙ the rate of temperature, κ
0
theheatfluxcoefficient,θthetemperature,andrtheradiation. Thederivationoftheweakformis
comparable with the derivation of the weak form for the balance momentum equation (2.9). First
thestrongformismultipliedwithascaler-valuedtestfunctionδθ, followedbyintegrationoverthe
domain V.
(cid:90) (cid:90) (cid:90)
ρcθ˙δθdV − κ ∇(cid:126) ·(∇(cid:126)θ)δθdV = rδθdV (2.11)
0
V V V
Applying the chain rule gives
(cid:16) (cid:17) (cid:16) (cid:17)
∇(cid:126) · δθ∇(cid:126)θ =∇(cid:126)θ·∇(cid:126) (δθ)+δθ∇(cid:126) ∇(cid:126)θ (2.12)
(cid:90) (cid:90) (cid:110) (cid:16) (cid:17)(cid:111) (cid:90)
ρcθ˙δθdV + κ ∇(cid:126)θ·∇(cid:126)(δθ)−∇(cid:126) · δθ∇(cid:126)θ dV = rδθdV (2.13)
0
V V V
(cid:90) (cid:90) (cid:90) (cid:16) (cid:17) (cid:90)
ρcθ˙δθdV + κ ∇(cid:126)θ·∇(cid:126)(δθ)dV − κ ∇(cid:126) · δθ∇(cid:126)θ dV = rδθdV (2.14)
0 0
V V V V
Using the divergence theorem (2.6) again gives
(cid:90) (cid:90) (cid:90) (cid:16) (cid:17) (cid:90)
ρcθ˙δθdV + κ ∇(cid:126)θ·∇(cid:126)(δθ)dV − κ (cid:126)n·(δθ∇(cid:126)θ) dS = rδθdV (2.15)
0 0
V V S V
Defining
κ ∇(cid:126)θ = −(cid:126)q (2.16)
0 c
(cid:126)q ·(cid:126)n = (cid:126)q (2.17)
c n
yields the final result
(cid:90) (cid:90) (cid:90) (cid:90)
ρcθ˙δθdV + κ ∇(cid:126)θ·∇(cid:126)(δθ)dV + (cid:126)q δθdS = rδθdV (2.18)
0 n
V V S V
2.3 Nonlocal balance equation
The nonlocal balance equation is formulated in strong form as
α¯−c(l)∇(cid:126) ·(∇(cid:126)α¯)=α (2.19)
With α being the local variable, α¯ being the nonlocal variable and c(l) being a constant term for
locality. Using the same strategy as before, multiplying the strong form with a test function δα¯
and then integrate over the volume, gives
(cid:90) (cid:90) (cid:90)
α¯δα¯dV − c(l)∇(cid:126) ·(∇(cid:126)α¯)δα¯dV = αδα¯dV (2.20)
V V V
and using the following relations
(cid:16) (cid:17)
∇(cid:126) · δα¯∇(cid:126)α¯ =∇(cid:126)α¯·∇(cid:126)(δα¯)+δα¯∇(cid:126) ·(∇(cid:126)α¯) (2.21)
6 2.4 Linearization of balance momentum
(cid:90) (cid:90)
∇(cid:126) ·(cid:126)adV = (cid:126)n·(cid:126)adS (2.22)
V S
c(l)∇(cid:126)α¯ =−(cid:126)qα (2.23)
c
(cid:126)qα·(cid:126)n=(cid:126)qα (2.24)
c n
Will lead to the weak formulation of the nonlocal balance equation
(cid:90) (cid:90) (cid:90) (cid:90)
α¯δα¯dV + c(l)∇(cid:126)α¯·∇(cid:126)(δα¯)+ (cid:126)qαδα¯dS = αδα¯dV (2.25)
n
V V S V
2.4 Linearization of balance momentum
In general, the above set of equation requires the application of an iterative solution scheme (e.g.
Newton-Rhapson). Thus, the system of equations have to be linearized. Recall formula (2.9), the
weak form of the balance momentum equation
(cid:90) (cid:90) (cid:90)
R = σσσ :∇(cid:126)(δ(cid:126)x)dV + (ρ(cid:126)x¨−ρ(cid:126)b)·δ(cid:126)xdV − (cid:126)t·δ(cid:126)xdS (2.26)
(cid:126)u
V V S
Assuming that (cid:126)b and (cid:126)t, the body force and traction vector at the boundary, respectively, are
independent from (cid:126)x and neglecting the inertia forces, the linearization of the residual (2.26) with
respect to the increment of displacement is given by
(cid:26)(cid:90) (cid:27)
∆ R =∆ σσσ :∇(cid:126)(δ(cid:126)x)dV (2.27)
(cid:126)u (cid:126)u (cid:126)u
V
which is best derived in quantities given in the referential configuration V . Here,
0
det(FFF) = J (2.28)
σσσ = J−1FFF ·SSS·FFFT (2.29)
∇(cid:126)(δ(cid:126)x) = ∇(cid:126) (δ(cid:126)x)·FFF−1 =FFF−T ·∇(cid:126) (δ(cid:126)x) (2.30)
0 0
dV = J dV (2.31)
0
Using these relations in (2.27) gives
(cid:26)(cid:90) (cid:27)
∆ R = ∆ J−1FFF ·SSS·FFFT :FFF−T ·∇(cid:126) (δ(cid:126)x)J dV
(cid:126)u (cid:126)u (cid:126)u 0 0
V0
(cid:26)(cid:90) (cid:27)
= ∆ FFF ·SSS :∇(cid:126) (δ(cid:126)x)dV (2.32)
(cid:126)u 0 0
V0
Because ∇(cid:126) (δ(cid:126)x) is independent of ∆ we get
0 (cid:126)u
(cid:90)
∆ R = (∆ FFF ·SSS+FFF ·∆ SSS):∇(cid:126) (δ(cid:126)x)dV (2.33)
(cid:126)u (cid:126)u (cid:126)u (cid:126)u 0 0
V0
Now, as
∆ FFF =∇(cid:126) (∆(cid:126)u) (2.34)
(cid:126)u 0
with ∆(cid:126)u the incremental displacement vector, relation (2.32) can be written as
(cid:90) (cid:16) (cid:17)
∆ R = ∇(cid:126) (∆(cid:126)u)·SSS+FFF ·∆ SSS :∇(cid:126) (δ(cid:126)x)dV (2.35)
(cid:126)u (cid:126)u 0 (cid:126)u 0 0
V0
2.4 Linearization of balance momentum 7
When now making use of the following relations
∇(cid:126) (∆(cid:126)u) = ∇(cid:126)(∆(cid:126)u)·FFF (2.36)
0
∇(cid:126) (δ(cid:126)x) = ∇(cid:126)(δ(cid:126)x)·FFF (2.37)
0
SSS = FFF−1·τττ ·FFF−T (2.38)
∆ FFF ·FFF−1 = ∇(cid:126)(∆(cid:126)u) (2.39)
(cid:126)u
whereτττ is the Kirchoff stress tensor and using
∆A−1 = −A−1·∆A·A−1 (2.40)
∆A−T = −A−T ·∆AT ·A−T (2.41)
we can calculateFFF ·∆ SSS as follows
(cid:126)u
FFF ·∆ SSS = FFF ·(cid:0)∆ FFF−1·τττ ·FFF−T +FFF−1·∆ τττ ·FFF−T +FFF−1·τττ ·∆ FFF−T(cid:1)
(cid:126)u (cid:126)u (cid:126)u (cid:126)u
= FFF ·∆ FFF−1·τττ ·FFF−T +∆ τττ ·FFF−T +τττ ·∆ FFF−T
(cid:126)u (cid:126)u (cid:126)u
= ∆ τττ ·FFF−T −FFF ·FFF−1·∆ FFF ·FFF−1·τττ ·FFF−T −τττ ·FFF−T ·∆ FFFT ·FFF−T
(cid:126)u (cid:126)u (cid:126)u
(cid:16) (cid:17)T
= ∆ τττ ·FFF−T −∇(cid:126)(∆(cid:126)u)·τττ ·FFF−T −τττ · ∇(cid:126)(∆(cid:126)u) ·FFF−T (2.42)
(cid:126)u
Whenpluggingtheabovesolutioninto(2.35)andcombiningitwithequations(2.30),(2.31),(2.36)
and (2.38), we obtain
(cid:90) (cid:18) (cid:16) (cid:17)T (cid:19)
∆ R = ∆ τττ ·FFF−T −τττ · ∇(cid:126)(∆(cid:126)u) ·FFF−T :∇(cid:126)(δ(cid:126)x)·FFF J−1dV (2.43)
(cid:126)u (cid:126)u (cid:126)u
V
Making use of the trace operation
(cid:16) (cid:17)
A:B =tr A·BT (2.44)
we finally find
(cid:90) (cid:18) (cid:16) (cid:17)T(cid:19)
∆ R = J−1∆ τττ −J−1τττ · ∇(cid:126)(∆(cid:126)u) :∇(cid:126)(δ(cid:126)x)dV
(cid:126)u (cid:126)u (cid:126)u
V
(cid:90) (cid:18) (cid:16) (cid:17)T(cid:19)
= J−1∆ τττ −σσσ· ∇(cid:126)(∆(cid:126)u) :∇(cid:126)(δ(cid:126)x)dV (2.45)
(cid:126)u
V
The linearization of the residual R with respect to an arbitrary scalar quantity α (e.g., the
(cid:126)u
temperature) that does not infect the deformation gradientFFF, is derived as follows. Starting from
equation (2.27):
(cid:26)(cid:90) (cid:27)
∆ R =∆ FFF ·SSS :∇(cid:126) (δ(cid:126)x)dV (2.46)
α (cid:126)u α 0 0
V0
∇(cid:126) (δ(cid:126)x) is again independent from ∆ , so we get
0 α
(cid:90)
∆ R = (∆ FFF ·SSS+FFF ·∆ SSS):∇(cid:126) (δ(cid:126)x)dV (2.47)
α (cid:126)u α α 0 0
V0
As said before, α does not infectFFF and will vanish from the formula
(cid:90)
∆ R = (FFF ·∆ SSS):∇(cid:126) (δ(cid:126)x)dV (2.48)
α (cid:126)u α 0 0
V0
Writing outFFF ·∆ SSS will give
α
FFF ·∆ SSS = ∆ τττ ·FFF−T −∆ FFF ·FFF−1·τττ ·FFF−T −τττ ·FFF−T ·∆ FFFT ·FFF−T
α α α α
= ∆ τττ ·FFF−T (2.49)
α
8 2.5 Linearization of heat balance equation and nonlocal balance equation
and combining these results with equation (2.47), will lead to
(cid:90)
∆ R = ∆ τττ ·FFF−T :∇(cid:126) (δ(cid:126)x)dV (2.50)
α (cid:126)u α 0 0
V0
Which leads to the final result when making use of (2.30), (2.31) and (2.36) as follows
(cid:90)
∆ R = J−1∆ τττ :∇(cid:126)(δ(cid:126)x)dV (2.51)
α (cid:126)u α
V
2.5 Linearization of heat balance equation and nonlocal bal-
ance equation
Tolinearizethebalanceofheatequationtogetherwiththenonlocalformulation,westartwiththe
residuals of the weak forms, given by equation (2.18) and (2.25)
(cid:90) (cid:90) (cid:90) (cid:90)
R = ρcθ˙δθdV − (cid:126)q ·∇(cid:126)(δθ)dV − rθδθdV + (cid:126)q δθdS
θ c c n
V V V S
(cid:90) (cid:90) (cid:90) (cid:90)
R = α¯δα¯dV + c(l)∇(cid:126)α¯·∇(cid:126)(δα¯)dV − αδα¯dV + (cid:126)qαδα¯dS
α¯ n
V V V S
The general form for both, the heat balance equation and the nonlocal formulation is given by
(cid:90) (cid:90) (cid:90) (cid:90)
R = f(α)δαdV − (cid:126)qα·∇(cid:126)(δα)dV − rαδαdV + qαδα¯dS , (2.52)
α c c n
V V V S
with
θ −θ
f(θ) = θ˙ = n+1 n (2.53)
∆t
(cid:126)qθ = −κ ∇(cid:126)(θ)α (2.54)
c 0
rθ = r (2.55)
c
f(α) = α (2.56)
(cid:126)qα = −c(l)∇(cid:126)(α) (2.57)
c
rα = r (2.58)
c
Assuming now, that (cid:126)qα is independent from α¯, the linearization of R is, in terms of the current
n α
configuration, given by
(cid:26)(cid:90) (cid:90) (cid:90) (cid:27)
∆R =∆ f(α)δαdV − (cid:126)qα·∇(cid:126)(δα)dV − rαδαdV (2.59)
α c c
V V V
Making use of
(cid:126)qα = J−1FFF ·(cid:126)qα (2.60)
c r
∇(cid:126)(δα) = FFF−T ·∇(cid:126) (δα) (2.61)
0
∇(cid:126)(α) = FFF−T ·∇(cid:126) (α) (2.62)
0
we can write equation (2.59) in terms of the reference configuration
(cid:26)(cid:90) (cid:90) (cid:90) (cid:27)
∆R =∆ f(α)δαdV − (cid:126)qα·∇(cid:126) (δα)dV − rαδαdV (2.63)
α 0 r 0 0 r 0
V0 V0 V0
Making use of (2.57) and of (2.62) we can write
(cid:126)qα =−cFFF−T ·∇(cid:126) (α) (2.64)
c 0
2.5 Linearization of heat balance equation and nonlocal balance equation 9
The flux can be given in terms of the reference configurations as follows:
(cid:126)qα = JFFF−1·(cid:126)qα
r c
= −cJFFF−1·FFF−T ·∇(cid:126) (α) (2.65)
0
or, applying the definition
(cid:126)qα :=FFF ·(cid:126)qα (2.66)
s r
as
(cid:126)qα = FFF−1·(−cJFFF−T ·∇(cid:126)α)
r
= FFF−1·(cid:126)qα (2.67)
s
Where the subscript s is a temporary state which will be used for simplicity reasons. Now, the
linearization of R can be reformulated as
α
(cid:26)(cid:90) (cid:90) (cid:90) (cid:27)
∆R =∆ f(α)δαdV − (FFF−1·(cid:126)qα)·∇(cid:126)(δα)dV − rαδαdV (2.68)
α 0 s 0 r 0
V0 V0 V0
The next step is calculating the linearization of R with respect to the increment of displacement
α
(cid:90) (cid:90)
∆ R =− (cid:0)∆ FFF−1·(cid:126)qα+FFF−1·∆ (cid:126)qα(cid:1)·∇(cid:126) (δα)dV − ∆ rαδαdV (2.69)
(cid:126)u α (cid:126)u s (cid:126)u s 0 0 (cid:126)u r 0
V0 V0
Because
∆ (FFF ·FFF−1)=∆ FFF ·FFF−1+FFF ·∆ FFF−1 =0 (2.70)
(cid:126)u (cid:126)u (cid:126)u
we find
∆ FFF−1 = −FFF−1·∆ FFF ·FFF−1
(cid:126)u (cid:126)u
= −FFF−1·∇(cid:126) (∆ )·FFF−1 (2.71)
0 (cid:126)u
Now, the linearization of (cid:126)qα with respect to the incremental displacement is given by
s
∆ (cid:126)qα =−c∆ (JFFF−T)·∇(cid:126) (α) (2.72)
(cid:126)u s (cid:126)u 0
When making use of
∆J = JFFF−T :∆FFFT
and using equation (2.41) we can write
∆ (JFFF−T) = ∆ JFFF−T +J ∆ FFF−T
(cid:126)u (cid:126)u (cid:126)u
= (JFFF−T :∆ FFFT)FFF−T −JFFF−T ·∆ FFFT ·FFF−T
(cid:126)u (cid:126)u
= (FFF−T :∆ FFFT)JFFF−T −FFF−T ·∆ FFFT ·JFFF−T
(cid:126)u (cid:126)u
= (cid:8)JFFF−T ⊗FFF−T −FFF−T ∆∆∆JFFF−T(cid:9)[∆ FFF] (2.73)
(cid:126)u
where we use the function∆∆∆ defined as follows (for proof see appendix A)
AAA·BBBT ·CCC ={AAA∆∆∆CCC} [BBB] (2.74)
and of
(AAA:BBB)CCC ={CCC⊗AAA} : [BBBT] (2.75)
Combined this will be
∆ (cid:126)qα = (cid:8)(cid:126)qα⊗FFF−T −FFF−T∆∆∆(cid:126)qα(cid:9)∆ FFF
(cid:126)u s s s (cid:126)u
= (cid:8)(cid:126)qα⊗FFF−T −FFF−T∆∆∆(cid:126)qα(cid:9)∇(cid:126) (∆(cid:126)u) (2.76)
s s 0
10 2.6 Summary
Finally, the linearization of R with respect to the incremental displacement in terms of the
α
reference configuration will be
∆ R = −(cid:90) (cid:16)FFF−1(cid:8)−III(cid:50)(cid:50)(cid:50)FFF−1·(cid:126)qα+(cid:126)qα⊗FFF−T −FFF−T∆∆∆(cid:126)qα(cid:9)∇(cid:126) (∆(cid:126)u)(cid:17)·∇(cid:126) (δα)dV
(cid:126)u α s s s 0 0 0
V0
(cid:90)
− ∆ rαδαdV (2.77)
(cid:126)u r 0
V0
Or in terms of the current configuration
(cid:90) (cid:16) (cid:17)
∆ R = − {−III(cid:50)(cid:50)(cid:50)(cid:126)qα+(cid:126)qα⊗III−III∆∆∆(cid:126)qα}∇(cid:126)(∆(cid:126)u) ·∇(cid:126)(δα)dV
(cid:126)u α c c c
V
(cid:90)
− ∆ rαδαdV (2.78)
(cid:126)u c
V
Where(cid:50)(cid:50)(cid:50) is a special function (see appendix A)
AAA·BBB·CCC = {AAA(cid:50)(cid:50)(cid:50)CCC} BBB (2.79)
AAA·BBB·(cid:126)c = {AAA(cid:50)(cid:50)(cid:50)(cid:126)c} BBB (2.80)
To derive the linearization of R with respect to α we start with equation (2.68)
α
(cid:90) (cid:90) (cid:90)
∆ R = ∆ f(α)δαdV − (FFF−1·∆ (cid:126)qα)·∇(cid:126) (δα)dV − ∆ rαδαdV (2.81)
α α α 0 α s 0 0 α r 0
V0 V0 V0
The linearization of (cid:126)qα to α is given by
s
∆ (cid:126)qα =−cJFFF−T ·∇(cid:126) (∆α) (2.82)
α s 0
When using this result in equation (2.80) we get the linearization of R with respect to α in the
α
reference configuration
(cid:90) (cid:90)
∆ R = ∆ f(α)δαdV − (−cJ)(FFF−1·FFF−T ·∇(cid:126) (∆α))·∇(cid:126) (δα)dV
α α α 0 0 0 0
V0 V0
(cid:90)
− ∆ rαδαdV (2.83)
α r 0
V0
or in terms of the current configuration as
(cid:90) (cid:90) (cid:90)
∆ R = ∆ f(α)δαdV + c∇(cid:126)(∆α)·∇(cid:126)(δα)dV − ∆ rαδαdV (2.84)
α α α α c
V V V
2.6 Summary
We started with the weak form of the momentum balance (2.9)
(cid:90) (cid:90) (cid:90)
R = σσσ :∇(cid:126)δ(cid:126)xdV + (ρ(cid:126)x¨−ρ(cid:126)b)·δ(cid:126)xdV − (cid:126)t·δ(cid:126)xdS =0 (2.85)
(cid:126)u
V V S
The linearization to (cid:126)u is given by
(cid:90) (cid:18) (cid:16) (cid:17)T(cid:19)
∆ R = J−1∆ τττ −σσσ· ∇(cid:126)(∆(cid:126)u) :∇(cid:126)(δ(cid:126)x)dV (2.86)
(cid:126)u (cid:126)u (cid:126)u
V
Now, making use of the function∆∆∆, we can write
(cid:40) (cid:41)
(cid:90) ∂τττ
∆ R = J−1 −σσσ∆∆∆III [∇(cid:126)(∆(cid:126)u)]:(∇(cid:126)(δ(cid:126)x)dV (2.87)
(cid:126)u (cid:126)u V ∂∇(cid:126)(∆(cid:126)u)