0
Before we start, some basic mathematic theorem.
Divergence
∇⋅v=vi,i
∇⋅T=Tij,iej
Gradient (dim+1)
∇ϕ=ϕ,iei
∇v=vi,jei⊗ej
∇T=Tij,kei⊗ej⊗ek
Curl (dim=)
∇×v=εijkvj,iek
∇×T=εijkTmj,iek⊗em
The product rule of differentiation
∵σ⋅ν=σijνjei
∴∇⋅(σ⋅ν)=(σijνj),i=σij,iνj+σijνj,i
also
(∇⋅σ)⋅ν=σij,iνj
∇ν:σ=νi,jσji
or as Prof. Zohdi suggested
∇ν:σ=νi,jσij
but since σ is symmetric
σij=σji
therefore
∇⋅(σ⋅ν)=(∇⋅σ)⋅ν+∇ν:σ(1)
Divergence Theorem
∫∂RϕndA=∫R∇⋅ϕdV
∫∂Rv⋅ndA=∫R∇⋅vdV(2)
∫∂RTndA=∫R∇⋅TdV
or
∫∂RϕnidA=∫Rϕ,idV
∫∂RvinidA=∫Rvi,idV
∫∂RTijnjdA=∫RTij,jdV
Stokes theorem
∫Cϕdx=∫Sn×∇ϕdA
∫Cv⋅dx=∫Sn⋅∇×vdA
∫CTdx=∫S(∇×T)TndA
or
∫Cϕdxi=∫Sεijknjϕ,kdA
∫Cvidxi=∫Sniεijkvk,jdA
∫CTijdxj=∫SεjpqTiq,pnjdA
1
The SMALL deformation of the body is governed by (strong form):
∇⋅(IE:∇u)+ρg=0
since
σ=IE:∇u
let
f=ρg
thus
∇⋅σ+f=0
so we can write
∫Ω(∇⋅σ+f)⋅νdΩ=0
using (1)
∫Ω(∇⋅(σ⋅ν)−∇ν:σ)dΩ+∫Ωf⋅νdΩ=0
using (2)
∫Ω∇ν:σdΩ=∫Ωf⋅νdΩ+∫∂Ωσ⋅ν⋅ndA
because
t=σ⋅n
therefore
∫Ω∇ν:σdΩ=∫Ωf⋅νdΩ+∫Γtt⋅νdA
so we have the weak form
Find u,u∣Γu=u∗, such that ∀ν,ν∣Γu=0∫Ω∇ν:IE:∇udΩ=∫Ωf⋅νdΩ+∫Γtt⋅νdA
where u∗ is the applied boundary displacement on Γu, t=t∗ on Γt
since the integrals must be finite, we improve the weak form as below
Find u∈H1(Ω),u∣Γu=u∗, such that ∀ν∈H1(Ω),ν∣Γu=0∫Ω∇ν:IE:∇udΩ=∫Ωf⋅νdΩ+∫Γtt⋅νdA
where
H1(Ω):=[H1(Ω)]3
explained as below
similar to the definition of u∈H1(Ω) which states
u∈H1(Ω) if ∣∣u∣∣H1(Ω)2:=∫Ωu,ju,jdΩ+∫ΩuudΩ<∞
the definition of u∈H1(Ω) states as below
u∈H1(Ω) if ∣∣u∣∣H1(Ω)2:=∫Ωui,jui,jdΩ+∫ΩuiuidΩ<∞
2
should u∗ be different from the applied boundary displacement on Γu, weak form can be stated as below
Find u∈H1(Ω), such that ∀ν∈H1(Ω)∫Ω∇ν:IE:∇udΩ=∫Ωf⋅νdΩ+∫Γtt⋅νdA+P⋆∫Γu(u∗−u)⋅νdA(3)
or
Find u∈H1(Ω), such that ∀ν∈H1(Ω)∫Ωνi,jIEijkluk,ldΩ=∫ΩfiνidΩ+∫ΓttiνidA+P⋆∫Γuui∗νi−uiνidA
where the (penalty) parameter P⋆ is a large positive number.
3
To have
{ϕ^i(ζ1,ζ2,ζ3)=1ϕ^i(ζ1,ζ2,ζ3)=0 , ζ1=ζi1,ζ2=ζi2,ζ3=ζi3 , others
for each ϕ^i 8 functions (for each of the 8 points) with 8 unknowns, aka. ζ1mζ2nζ3k where m,n,k=0,1 can be derived, leading to the only solution as following
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧ϕ^1=81(1−ζ1)(1−ζ2)(1−ζ3)ϕ^2=81(1+ζ1)(1−ζ2)(1−ζ3)ϕ^3=81(1+ζ1)(1+ζ2)(1−ζ3)ϕ^4=81(1−ζ1)(1+ζ2)(1−ζ3)ϕ^5=81(1−ζ1)(1−ζ2)(1+ζ3)ϕ^6=81(1+ζ1)(1−ζ2)(1+ζ3)ϕ^7=81(1+ζ1)(1+ζ2)(1+ζ3)ϕ^8=81(1−ζ1)(1+ζ2)(1+ζ3)
4
IE has the following symmetries
IEijkl=IEklij=IEjikl=IEijlk
allowing us to rewrite (3) in a more compact matrix form
∫Ω([D]{ν})T[IE]([D]{u})dΩ=∫Ω{ν}T{f}dΩ+∫Γt{ν}T{t∗}dA+P⋆∫Γu{ν}T{u∗−u}dA(4)
where
[D]:=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡∂x1∂00∂x2∂0∂x3∂0∂x2∂0∂x1∂∂x3∂000∂x3∂0∂x2∂∂x1∂⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,{u}:=⎩⎪⎨⎪⎧u1u2u3⎭⎪⎬⎪⎫,{f}:=⎩⎪⎨⎪⎧f1f2f3⎭⎪⎬⎪⎫,{t∗}:=⎩⎪⎨⎪⎧t1∗t2∗t3∗⎭⎪⎬⎪⎫,(5)
and [IE] is the elastic stiffness matrix of the material.
[IE]:=⎣⎢⎢⎢⎢⎢⎢⎢⎡IE1111IE2211IE3311IE1211IE2311IE1311IE1122IE2222IE3322IE1222IE2322IE1322IE1133IE2233IE3333IE1233IE2333IE1333IE1112IE2212IE3312IE1212IE2312IE1312IE1123IE2223IE3323IE1223IE2323IE1323IE1113IE2213IE3313IE1213IE2313IE1313⎦⎥⎥⎥⎥⎥⎥⎥⎤
Further rewrite {uh}
{uh}=[ϕ]{a}(6)
where
{a}:=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1a2a3...a3N⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫,[ϕ]:=⎣⎢⎡ϕ1ϕ2...ϕNϕ1ϕ2...ϕNϕ1ϕ2...ϕN⎦⎥⎤,
choose ν with the same basis, but a different linear combination
{νh}=[ϕ]{b}(7)
with (6),(7) we can rewrite (4)
∫Ω([D][ϕ]{b})T[IE]([D][ϕ]{a})dΩ=∫Ω([ϕ]{b})T{f}dΩ+∫Γt([ϕ]{b})T{t∗}dA+P⋆∫Γu([ϕ]{b})T{u∗−[ϕ]{a}}dA(8)
We can rewrite (8)
{b}T{[K]{a}−{R}}=0
where
[K]:=∫Ω([D][ϕ])T[IE]([D][ϕ])dΩ+P⋆∫Γu[ϕ]T[ϕ]dA(9)
{R}:=∫Ω[ϕ]T{f}dΩ+∫Γt[ϕ]T{t∗}dA+P⋆∫Γu[ϕ]T{u∗}dA(10)
Consider
[K]=e=1∑Ne[K]e
we can rewrite (9),(10) for e=1,2..,Ne
[K]e:=∫Ωe([D][ϕ])T[IE]([D][ϕ])dΩe+P⋆∫Γu,e[ϕ]T[ϕ]dAe(11)
{R}e:=∫Ωe[ϕ]T{f}dΩe+∫Γt,e[ϕ]T{t∗}dAe+P⋆∫Γu,e[ϕ]T{u∗}dAe(12)
where
Γu,e=Γu∩∂Ωe
Γt,e=Γt∩∂Ωe
5
In 3D,
[ϕ^]:=⎣⎢⎡ϕ^1ϕ^2...ϕ^8ϕ^1ϕ^2...ϕ^8ϕ^1ϕ^2...ϕ^8⎦⎥⎤
express [D] in terms ζ1,ζ2,ζ3
[D(ϕ(x1,x2,x3))]=[D^ϕ(Mx1(ζ1,ζ2,ζ3),Mx2(ζ1,ζ2,ζ3),Mx3(ζ1,ζ2,ζ3))]
so that
[D^]=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡ζi∂∂x1ζi00ζi∂∂x2ζi0ζi∂∂x3ζi0ζi∂∂x2ζi0ζi∂∂x1ζiζi∂∂x3ζi000ζi∂∂x3ζi0ζi∂∂x2ζiζi∂∂x1ζi⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
therefore
[D^][ϕ^]=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡ζi∂ϕ^1∂x1ζiζi∂ϕ^1∂x2ζiζi∂ϕ^1∂x3ζiζi∂ϕ^2∂x1ζiζi∂ϕ^2∂x2ζiζi∂ϕ^2∂x3ζi.........ζi∂ϕ^8∂x1ζiζi∂ϕ^8∂x2ζiζi∂ϕ^8∂x3ζiζi∂ϕ^1∂x2ζiζi∂ϕ^1∂x1ζiζi∂ϕ^1∂x3ζiζi∂ϕ^2∂x2ζiζi∂ϕ^2∂x1ζiζi∂ϕ^2∂x3ζi.........ζi∂ϕ^8∂x2ζiζi∂ϕ^8∂x1ζiζi∂ϕ^8∂x3ζiζi∂ϕ^1∂x3ζiζi∂ϕ^1∂x2ζiζi∂ϕ^1∂x1ζiζi∂ϕ^2∂x3ζiζi∂ϕ^2∂x2ζiζi∂ϕ^2∂x1ζi.........ζi∂ϕ^8∂x3ζiζi∂ϕ^8∂x2ζiζi∂ϕ^8∂x1ζi⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
with Gaussian Quadrature
∫0LF(x)dx=∫−11F(x(ζ))J(ζ)dζ=1∑GwiF(ζi)J(ζi)
we can rewrite (11),(12)
[K]e=1∑G1∑G1∑Gwqwrws([D^][ϕ^])T[IE^][D^][ϕ^]J+P⋆1∑G1∑Gwqwr[ϕ]T[ϕ]Js
{R}e=1∑G1∑G1∑Gwqwrws[ϕ]T{f}J+1∑G1∑Gwqwr[ϕ]T{t∗}Js+P⋆1∑G1∑Gwqwr[ϕ]T{u∗}Js
where
F:=∇x(ζ1,ζ2,ζ3)
J:=∣F∣
and
Js=JF−T⋅N⋅n
derived as below
Nanson’s formula
ndAe=JF−T⋅NdAe^
where Ae is an area of a region in the current configuration, Ae^ is the same area in the reference configuration, and n is the outward normal to the area element in the current configuration while N is the outward normal in the reference configuration.
multiply n on both sides
dAe=JF−T⋅N⋅ndAe^
therefore
Js=JF−T⋅N⋅n
6
6.1
∇⋅(IK⋅∇θ)+ρs=0
let the heat flux density q
q=IK⋅∇θ
and let
f=ρs
therefore
∇⋅q+f=0
so we can write
∫Ω(∇⋅q+f)⋅νdΩ=0(13)
now consider
(∇⋅q)⋅ν=qi,iν
∇⋅(qν)=(qiν),i=qi,iν+qiν,i
(∇ν)⋅q=qiν,i
therefore we can write (13) as
∫Ω(∇⋅(qν)−(∇ν)⋅q)dΩ+∫ΩfνdΩ=0
∫Ω(∇ν)⋅qdΩ=∫ΩfνdΩ+∫Ω∇⋅(qν)
using (2)
∫Ω(∇ν)⋅qdΩ=∫ΩfνdΩ+∫∂Ωq⋅nνdA
so we have the weak form
Find θ∈H1(Ω),θ∣Γθ=θ∗, such that ∀ν∈H1(Ω),ν∣Γθ=0∫Ω(∇ν)⋅IK⋅∇θdΩ=∫ΩfνdΩ+∫∂Ωq⋅nνdA
6.2
should θ∗ be different from the applied boundary displacement on Γθ, weak form can be stated as below
Find θ∈H1(Ω),θ∣Γθ=θ∗, such that ∀ν∈H1(Ω),ν∣Γθ=0∫Ω(∇ν)⋅IK⋅∇θdΩ=∫ΩfνdΩ+∫∂Ωq⋅nνdA+P⋆∫Γθ(θ∗−θ)νdA(14)
where the (penalty) parameter P⋆ is a large positive number.
6.3
same as 3
6.4
rewrite θh
θh=[ϕ]{a}=ϕiai(15)
where
{a}:=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1a2a3...aN⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫,[ϕ]:=[ϕ1ϕ2...ϕN],
choose ν with the same basis, but a different linear combination
νh=[ϕ]{b}=ϕjbj(16)
with (15),(16) we can rewrite (14)
∫Ω(∇[ϕ]{b})⋅IK⋅∇[ϕ]{a}dΩ=∫Ωf[ϕ]{b}dΩ+∫∂Ωq⋅n[ϕ]{b}dA+P⋆∫Γθ(θ∗−[ϕ]{a})⋅[ϕ]{b}dA(17)
We can rewrite (17)
{b}T{K{a}−{R}}=0
where
K:=∫Ω(∇[ϕ])⋅IK⋅∇[ϕ]dΩ+P⋆∫Γθ[ϕ]T[ϕ]dA(18)
{R}:=∫Ω[ϕ]TfdΩ+∫∂Ω[ϕ]Tq⋅ndA+P⋆∫Γθ[ϕ]Tθ∗dA(19)
Consider
K=e=1∑NeKe
we can rewrite (18),(19) for e=1,2..,Ne
Ke:=∫Ωe(∇[ϕ])⋅IK⋅∇[ϕ]dΩe+P⋆∫Γu,e[ϕ]T[ϕ]dAe
{R}e:=∫Ωe[ϕ]TfdΩe+∫∂Ω,e[ϕ]Tq⋅ndAe+P⋆∫Γu,e[ϕ]Tθ∗dAe
where
Γθ,e=Γθ∩∂Ωe
7
Key difference is that the diemsion of the equations derived for thermodynamics are generally lower than those for linear elasticity.
This is because θ and s are scalers rather than vectors as u and g, and IK is a second order tensor rather than a forth order tensor as IE.
Reference