We will now derive the finite element equations for the updated Lagrangian formulation for three-dimensional problems in solid mechanics.
Updated Lagrangian governing equations
edit
The updated Lagrangian equations are
ρ
0
(
X
)
=
J
(
X
,
t
)
ρ
(
X
,
t
)
ρ
v
˙
=
∇
⋅
σ
T
+
ρ
b
σ
=
σ
T
ρ
e
˙
=
σ
:
d
−
∇
⋅
q
+
ρ
s
L
φ
σ
=
G
(
σ
,
d
,
T
,
…
)
d
=
1
2
(
l
+
l
T
)
{\displaystyle {\begin{aligned}\rho _{0}(\mathbf {X} )&=J(\mathbf {X} ,t)~\rho (\mathbf {X} ,t)\\\rho ~{\dot {\mathbf {v} }}&={\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}^{T}+\rho ~\mathbf {b} \\{\boldsymbol {\sigma }}&={\boldsymbol {\sigma }}^{T}\\\rho ~{\dot {e}}&={\boldsymbol {\sigma }}:{\boldsymbol {d}}-{\boldsymbol {\nabla }}\cdot \mathbf {q} +\rho ~s\\{\mathcal {L}}_{\varphi }~{\boldsymbol {\sigma }}&={\boldsymbol {G}}({\boldsymbol {\sigma }},{\boldsymbol {d}},T,\dots )\\{\boldsymbol {d}}&={\frac {1}{2}}({\boldsymbol {l}}+{\boldsymbol {l}}^{T})\end{aligned}}}
The boundary conditions are
σ
T
⋅
n
=
t
on
∂
Ω
t
v
=
v
¯
on
∂
Ω
v
{\displaystyle {\begin{aligned}{\boldsymbol {\sigma }}^{T}\cdot \mathbf {n} &=\mathbf {t} \qquad {\text{on}}~~\partial \Omega _{t}\\\mathbf {v} &={\overline {\mathbf {v} }}\qquad {\text{on}}~~\partial \Omega _{v}\end{aligned}}}
The initial conditions are
v
(
X
,
0
)
=
v
0
(
X
)
σ
(
X
,
0
)
=
σ
0
(
X
)
or
u
(
X
,
0
)
=
u
0
(
X
)
{\displaystyle {\begin{aligned}\mathbf {v} (\mathbf {X} ,0)&=\mathbf {v} _{0}(\mathbf {X} )\\{\boldsymbol {\sigma }}(\mathbf {X} ,0)&={\boldsymbol {\sigma }}_{0}(\mathbf {X} )\qquad {\text{or}}\qquad \mathbf {u} (\mathbf {X} ,0)&=\mathbf {u} _{0}(\mathbf {X} )\end{aligned}}}
Note that the dependent variables in this case are all expressed as functions of the Lagrangian (material) coordinates, i.e.,
Velocity (
v
(
X
,
t
)
{\displaystyle \mathbf {v} (\mathbf {X} ,t)}
)
Cauchy stress (
σ
(
X
,
t
)
{\displaystyle {\boldsymbol {\sigma }}(\mathbf {X} ,t)}
)
Rate of deformation (
d
(
v
,
t
)
{\displaystyle {\boldsymbol {d}}(\mathbf {v} ,t)}
)
Also note that instead of
σ
{\displaystyle {\boldsymbol {\sigma }}}
and
d
{\displaystyle {\boldsymbol {d}}}
, we may use
S
{\displaystyle {\boldsymbol {S}}}
and
E
{\displaystyle {\boldsymbol {E}}}
in our calculations and then transform to
σ
{\displaystyle {\boldsymbol {\sigma }}}
.
Another name for the weak form is the Principle of Virtual Power (or Principle of Virtual Work if rates are not involved).
The strong form consists of the balance of linear momentum and the traction boundary conditions. That is,
∇
⋅
σ
T
+
ρ
b
=
ρ
a
in
Ω
σ
T
⋅
n
=
t
on
∂
Ω
t
{\displaystyle {\begin{aligned}{\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}^{T}+\rho ~\mathbf {b} &=\rho ~\mathbf {a} \qquad {\text{in}}~~\Omega \\{\boldsymbol {\sigma }}^{T}\cdot \mathbf {n} &=\mathbf {t} \qquad {\text{on}}~~\partial \Omega _{t}\end{aligned}}}
where
a
=
v
˙
=
u
¨
{\displaystyle \mathbf {a} ={\dot {\mathbf {v} }}={\ddot {\mathbf {u} }}}
is the acceleration.
The weak form of this equation can be written as
∫
Ω
ρ
a
⋅
w
dV
+
∫
Ω
σ
:
(
∇
w
)
dV
=
∫
Ω
ρ
b
⋅
w
dV
+
∫
∂
Ω
t
t
⋅
w
dA
{\displaystyle {\int _{\Omega }\rho ~\mathbf {a} \cdot \mathbf {w} ~{\text{dV}}+\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\mathbf {w} )~{\text{dV}}=\int _{\Omega }\rho ~\mathbf {b} \cdot \mathbf {w} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}\mathbf {t} \cdot \mathbf {w} ~{\text{dA}}}}
where
w
{\displaystyle \mathbf {w} }
is a vector valued weighting function (also called a test function ) that is zero at the points on the boundary where the essential boundary conditions on
v
{\displaystyle \mathbf {v} }
are applied. Therefore we can think of
w
{\displaystyle \mathbf {w} }
as a variation of
v
{\displaystyle \mathbf {v} }
.
We proceed in the usual manner to derive the weak form. We multiply the differential equation with the weighting function and integrate over the volume to get
∫
Ω
w
⋅
(
∇
⋅
σ
T
+
ρ
b
−
ρ
a
)
dV
=
0
{\displaystyle \int _{\Omega }\mathbf {w} \cdot \left({\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}^{T}+\rho ~\mathbf {b} -\rho ~\mathbf {a} \right)~{\text{dV}}=0}
Using the identity
∇
⋅
(
A
T
⋅
b
)
=
b
⋅
(
∇
⋅
A
)
+
A
:
(
∇
b
)
{\displaystyle {\boldsymbol {\nabla }}\cdot ({\boldsymbol {A}}^{T}\cdot \mathbf {b} )=\mathbf {b} \cdot ({\boldsymbol {\nabla }}\cdot {\boldsymbol {A}})+{\boldsymbol {A}}:({\boldsymbol {\nabla }}\mathbf {b} )}
and noting that
σ
=
σ
T
{\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {\sigma }}^{T}}
(conservation of angular momentum), we
have
∫
Ω
w
⋅
(
∇
⋅
σ
T
)
dV
=
∫
Ω
∇
⋅
(
σ
⋅
w
)
dV
−
∫
Ω
σ
:
(
∇
w
)
dV
{\displaystyle \int _{\Omega }\mathbf {w} \cdot ({\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}^{T})~{\text{dV}}=\int _{\Omega }{\boldsymbol {\nabla }}\cdot ({\boldsymbol {\sigma }}\cdot \mathbf {w} )~{\text{dV}}-\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\mathbf {w} )~{\text{dV}}}
From the divergence theorem we have
∫
Ω
∇
⋅
(
σ
⋅
w
)
dV
=
∫
∂
Ω
n
⋅
(
σ
⋅
w
)
dA
{\displaystyle \int _{\Omega }{\boldsymbol {\nabla }}\cdot ({\boldsymbol {\sigma }}\cdot \mathbf {w} )~{\text{dV}}=\int _{\partial \Omega }\mathbf {n} \cdot ({\boldsymbol {\sigma }}\cdot \mathbf {w} )~{\text{dA}}}
where
n
{\displaystyle \mathbf {n} }
is the outward unit normal to the boundary. So,
∫
Ω
w
⋅
(
∇
⋅
σ
)
dV
=
∫
∂
Ω
n
⋅
(
σ
⋅
w
)
dA
−
∫
Ω
σ
:
(
∇
w
)
dV
{\displaystyle \int _{\Omega }\mathbf {w} \cdot ({\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }})~{\text{dV}}=\int _{\partial \Omega }\mathbf {n} \cdot ({\boldsymbol {\sigma }}\cdot \mathbf {w} )~{\text{dA}}-\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\mathbf {w} )~{\text{dV}}}
Now, on the boundary
n
⋅
(
σ
⋅
w
)
≡
n
i
σ
i
j
w
j
=
σ
j
i
n
i
w
j
≡
(
σ
T
⋅
n
)
⋅
w
=
t
⋅
w
{\displaystyle \mathbf {n} \cdot ({\boldsymbol {\sigma }}\cdot \mathbf {w} )\equiv n_{i}~\sigma _{ij}~w_{j}=\sigma _{ji}~n_{i}~w_{j}\equiv ({\boldsymbol {\sigma }}^{T}\cdot \mathbf {n} )\cdot \mathbf {w} =\mathbf {t} \cdot \mathbf {w} }
On the part of the boundary where no tractions are applied, we have
t
=
0
{\displaystyle \mathbf {t} =\mathbf {0} }
and hence
n
⋅
(
σ
⋅
w
)
=
0
{\displaystyle \mathbf {n} \cdot ({\boldsymbol {\sigma }}\cdot \mathbf {w} )=\mathbf {0} }
. Therefore we can write
∫
Ω
w
⋅
(
∇
⋅
σ
)
dV
=
∫
∂
Ω
t
t
⋅
w
dA
−
∫
Ω
σ
:
(
∇
w
)
dV
{\displaystyle \int _{\Omega }\mathbf {w} \cdot ({\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }})~{\text{dV}}={\int _{\partial \Omega }}_{t}\mathbf {t} \cdot \mathbf {w} ~{\text{dA}}-\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\mathbf {w} )~{\text{dV}}}
where
∂
Ω
t
{\displaystyle \partial \Omega _{t}}
is the part of the boundary where the tractions are applied.
Plugging this expression into the weighted average form of the balance of linear momentum we get
∫
∂
Ω
t
t
⋅
w
dA
−
∫
Ω
σ
:
(
∇
w
)
dV
+
∫
Ω
w
⋅
(
ρ
b
−
ρ
a
)
dV
=
0
{\displaystyle {\int _{\partial \Omega }}_{t}\mathbf {t} \cdot \mathbf {w} ~{\text{dA}}-\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\mathbf {w} )~{\text{dV}}+\int _{\Omega }\mathbf {w} \cdot \left(\rho ~\mathbf {b} -\rho ~\mathbf {a} \right)~{\text{dV}}=0}
Rearranging terms we get the required weak form
∫
Ω
ρ
a
⋅
w
dV
+
∫
Ω
σ
:
(
∇
w
)
dV
=
∫
Ω
ρ
b
⋅
w
dV
+
∫
∂
Ω
t
t
⋅
w
dA
{\displaystyle \int _{\Omega }\rho ~\mathbf {a} \cdot \mathbf {w} ~{\text{dV}}+\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\mathbf {w} )~{\text{dV}}=\int _{\Omega }\rho ~\mathbf {b} \cdot \mathbf {w} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}\mathbf {t} \cdot \mathbf {w} ~{\text{dA}}}
Physical interpretation of terms
edit
To provide a physical interpretation of the four terms in the weak form, we write the weighting function as a variation of the velocity, i.e.,
w
=
δ
v
{\displaystyle \mathbf {w} =\delta \mathbf {v} }
. Then we have the familiar principle of virtual power:
∫
Ω
ρ
a
⋅
δ
v
dV
+
∫
Ω
σ
:
(
∇
δ
v
)
dV
=
∫
Ω
ρ
b
⋅
δ
v
dV
+
∫
∂
Ω
t
t
⋅
δ
v
dA
{\displaystyle \int _{\Omega }\rho ~\mathbf {a} \cdot \delta \mathbf {v} ~{\text{dV}}+\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\delta \mathbf {v} )~{\text{dV}}=\int _{\Omega }\rho ~\mathbf {b} \cdot \delta \mathbf {v} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}\mathbf {t} \cdot \delta \mathbf {v} ~{\text{dA}}}
The first term on the left can then be identified as the virtual kinetic power, the second term on the left as the virtual internal power, and the sum of the two terms on the right as the virtual external power. In other words, we have
δ
P
kin
:=
∫
Ω
ρ
a
⋅
δ
v
dV
δ
P
int
:=
∫
Ω
σ
:
(
∇
δ
v
)
dV
δ
P
ext
:=
∫
Ω
ρ
b
⋅
δ
v
dV
+
∫
∂
Ω
t
t
⋅
δ
v
dA
{\displaystyle {\begin{aligned}\delta P_{\text{kin}}&:=\int _{\Omega }\rho ~\mathbf {a} \cdot \delta \mathbf {v} ~{\text{dV}}\\\delta P_{\text{int}}&:=\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\delta \mathbf {v} )~{\text{dV}}\\\delta P_{\text{ext}}&:=\int _{\Omega }\rho ~\mathbf {b} \cdot \delta \mathbf {v} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}\mathbf {t} \cdot \delta \mathbf {v} ~{\text{dA}}\end{aligned}}}
and
δ
P
kin
+
δ
P
int
=
δ
P
ext
.
{\displaystyle \delta P_{\text{kin}}+\delta P_{\text{int}}=\delta P_{\text{ext}}~.}
We can see why the second terms is identified as the virtual internal power
as follows. Note that
∇
δ
v
=
δ
l
=
δ
(
d
+
w
)
{\displaystyle {\boldsymbol {\nabla }}\delta \mathbf {v} =\delta {\boldsymbol {l}}=\delta ({\boldsymbol {d}}+{\boldsymbol {w}})}
where
l
{\displaystyle {\boldsymbol {l}}}
is the velocity gradient,
d
{\displaystyle {\boldsymbol {d}}}
is symmetric rate of deformation tensor, and
w
{\displaystyle {\boldsymbol {w}}}
is the skew symmetric spin tensor. Therefore, since the product of a symmetric and a skew symmetric tensor is zero, we have
σ
:
∇
δ
v
=
σ
:
δ
d
+
σ
:
δ
w
=
σ
:
δ
d
{\displaystyle {\boldsymbol {\sigma }}:{\boldsymbol {\nabla }}\delta \mathbf {v} ={\boldsymbol {\sigma }}:\delta {\boldsymbol {d}}+{\boldsymbol {\sigma }}:\delta {\boldsymbol {w}}={\boldsymbol {\sigma }}:\delta {\boldsymbol {d}}}
Hence,
δ
P
int
=
∫
Ω
σ
:
δ
d
dV
{\displaystyle \delta P_{\text{int}}=\int _{\Omega }{\boldsymbol {\sigma }}:\delta {\boldsymbol {d}}~{\text{dV}}}
This has the same form as the expression for internal power that we discussed
earlier.
Finite element discretization
edit
Let us go back to the original weak form that we derived, i.e.,
∫
Ω
ρ
a
⋅
w
dV
+
∫
Ω
σ
:
(
∇
w
)
dV
=
∫
Ω
ρ
b
⋅
w
dV
+
∫
∂
Ω
t
t
⋅
w
dA
{\displaystyle \int _{\Omega }\rho ~\mathbf {a} \cdot \mathbf {w} ~{\text{dV}}+\int _{\Omega }{\boldsymbol {\sigma }}:({\boldsymbol {\nabla }}\mathbf {w} )~{\text{dV}}=\int _{\Omega }\rho ~\mathbf {b} \cdot \mathbf {w} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}\mathbf {t} \cdot \mathbf {w} ~{\text{dA}}}
To start the discretization process we have to first choose an approximate solution
v
h
{\displaystyle \mathbf {v} _{h}}
from the space
S
h
{\displaystyle {\mathcal {S}}_{h}}
of admissible finite element trial functions . These functions must satisfy at least
C
0
{\displaystyle C^{0}}
continuity and the essential boundary condition
v
=
v
¯
{\displaystyle \mathbf {v} ={\overline {\mathbf {v} }}}
on
∂
Ω
v
{\displaystyle {\partial \Omega }_{v}}
.
More formally, we choose
v
h
(
X
,
t
)
∈
S
h
{\displaystyle \mathbf {v} _{h}(\mathbf {X} ,t)\in {\mathcal {S}}_{h}}
where
S
h
=
{
v
|
v
∈
C
0
(
X
)
,
v
=
v
¯
on
∂
Ω
v
}
{\displaystyle {\mathcal {S}}_{h}=\{\mathbf {v} ~|~\mathbf {v} \in C^{0}(\mathbf {X} ),~~\mathbf {v} ={\overline {\mathbf {v} }}~{\text{on}}~\partial \Omega _{v}\}}
If we discretize the body into a number of elements with
n
{\displaystyle n}
nodes,
then we can write the approximate solution as
v
h
(
X
,
t
)
=
∑
I
=
1
n
N
I
(
X
)
v
I
(
t
)
{\displaystyle \mathbf {v} _{h}(\mathbf {X} ,t)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {v} _{I}(t)}
The quantities
N
I
(
X
)
{\displaystyle N_{I}(\mathbf {X} )}
are the shape functions which provide the required level of continuity and are chosen so that
v
h
∈
S
h
{\displaystyle \mathbf {v} _{h}\in {\mathcal {S}}_{h}}
.
We also have to choose appropriate weighting functions (
w
h
{\displaystyle \mathbf {w} _{h}}
) which are also known as test functions .
In Galerkin finite elements we choose these functions from the same space as the trial
functions with the additional restriction that these functions go to zero
at points on the boundary where essential boundary conditions are applied.
Formally, choose
w
h
(
X
)
∈
S
t
{\displaystyle \mathbf {w} _{h}(\mathbf {X} )\in {\mathcal {S}}_{t}}
where
S
t
=
{
w
|
w
∈
C
0
(
X
)
,
w
=
0
on
∂
Ω
v
}
{\displaystyle {\mathcal {S}}_{t}=\{\mathbf {w} ~|~\mathbf {w} \in C^{0}(\mathbf {X} ),~~\mathbf {w} =\mathbf {0} ~{\text{on}}~\partial \Omega _{v}\}}
For finite element analysis, these weighting functions have the same form
as the trial function, i.e.,
w
h
(
X
)
=
∑
I
=
1
n
N
I
(
X
)
w
I
{\displaystyle \mathbf {w} _{h}(\mathbf {X} )=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {w} _{I}}
Instead of a trial function for the velocity, we can start with one for the
motion and then take time derivatives. In this case, we choose the motion
to be approximated by
x
h
(
X
,
t
)
=
∑
I
=
1
n
N
I
(
X
)
x
I
(
t
)
{\displaystyle \mathbf {x} _{h}(\mathbf {X} ,t)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {x} _{I}(t)}
where
N
I
(
X
)
{\displaystyle N_{I}(\mathbf {X} )}
are the nodal shape functions (note that these are functions of
X
{\displaystyle \mathbf {X} }
) and
x
I
(
t
)
{\displaystyle \mathbf {x} _{I}(t)}
are the positions of nodes
i
{\displaystyle i}
in the current configuration. Therefore, nodes remain coincident with material point labels at all times.
Therefore, in the reference configuration, we have
x
h
(
X
,
0
)
=
∑
I
=
1
n
N
I
(
X
)
x
I
(
0
)
=
∑
I
=
1
n
N
I
(
X
)
X
I
{\displaystyle \mathbf {x} _{h}(\mathbf {X} ,0)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {x} _{I}(0)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {X} _{I}}
The shape functions are also chosen such that
N
I
(
X
J
)
=
δ
I
J
{\displaystyle N_{I}(\mathbf {X} _{J})=\delta _{IJ}}
Then the displacement
u
=
x
−
X
{\displaystyle \mathbf {u} =\mathbf {x} -\mathbf {X} }
can be approximated by
u
h
(
X
,
t
)
=
x
h
(
X
,
t
)
−
X
h
=
∑
I
=
1
n
N
I
(
X
)
[
x
I
(
t
)
−
X
I
]
=
∑
I
=
1
n
N
I
(
X
)
u
I
(
t
)
{\displaystyle \mathbf {u} _{h}(\mathbf {X} ,t)=\mathbf {x} _{h}(\mathbf {X} ,t)-\mathbf {X} _{h}=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\left[\mathbf {x} _{I}(t)-\mathbf {X} _{I}\right]=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {u} _{I}(t)}
The velocity
v
{\displaystyle \mathbf {v} }
is given by the time derivative of
u
{\displaystyle \mathbf {u} }
keeping
X
{\displaystyle \mathbf {X} }
fixed. So we have
v
h
(
X
,
t
)
=
∑
I
=
1
n
N
I
(
X
)
u
˙
I
(
t
)
=
∑
I
=
1
n
N
I
(
X
)
v
I
(
t
)
{\displaystyle \mathbf {v} _{h}(\mathbf {X} ,t)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~{\dot {\mathbf {u} }}_{I}(t)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {v} _{I}(t)}
Note that this approximation has the same form as the trial function for the velocity that we started with.
The acceleration
a
{\displaystyle \mathbf {a} }
is given by the time derivative of
v
{\displaystyle \mathbf {v} }
. Hence
a
h
(
X
,
t
)
=
∑
I
=
1
n
N
I
(
X
)
u
¨
I
(
t
)
=
∑
I
=
1
n
N
I
(
X
)
v
˙
I
(
t
)
=
∑
I
=
1
n
N
I
(
X
)
a
I
(
t
)
{\displaystyle \mathbf {a} _{h}(\mathbf {X} ,t)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~{\ddot {\mathbf {u} }}_{I}(t)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~{\dot {\mathbf {v} }}_{I}(t)=\sum _{I=1}^{n}N_{I}(\mathbf {X} )~\mathbf {a} _{I}(t)}
The spatial velocity gradient
l
{\displaystyle {\boldsymbol {l}}}
is given by
l
(
x
,
t
)
=
∇
v
(
x
,
t
)
{\displaystyle {\boldsymbol {l}}(\mathbf {x} ,t)={\boldsymbol {\nabla }}\mathbf {v} (\mathbf {x} ,t)}
Note that the derivatives are with respect to
x
{\displaystyle \mathbf {x} }
and not
X
{\displaystyle \mathbf {X} }
.
Therefore, the approximation that we seek has the form
l
h
(
x
,
t
)
=
∑
I
=
1
n
∇
[
N
I
(
X
)
v
I
(
t
)
]
{\displaystyle {\boldsymbol {l}}_{h}(\mathbf {x} ,t)=\sum _{I=1}^{n}{\boldsymbol {\nabla }}[N_{I}(\mathbf {X} )~\mathbf {v} _{I}(t)]}
Using the identity
∇
(
φ
a
)
=
a
⊗
(
∇
φ
)
+
φ
∇
a
{\displaystyle {\boldsymbol {\nabla }}(\varphi ~\mathbf {a} )=\mathbf {a} \otimes ({\boldsymbol {\nabla }}\varphi )+\varphi ~{\boldsymbol {\nabla }}\mathbf {a} }
we get
l
h
(
x
,
t
)
=
∑
I
=
1
n
v
I
(
t
)
⊗
∇
N
I
(
X
)
{\displaystyle {\boldsymbol {l}}_{h}(\mathbf {x} ,t)=\sum _{I=1}^{n}\mathbf {v} _{I}(t)\otimes {\boldsymbol {\nabla }}N_{I}(\mathbf {X} )}
In index notation,
[
l
i
j
]
h
=
∑
I
=
1
n
v
I
i
∂
N
I
(
X
,
t
)
∂
x
j
=
v
I
i
N
I
,
j
{\displaystyle [l_{ij}]_{h}=\sum _{I=1}^{n}v_{Ii}~{\frac {\partial N_{I}(\mathbf {X} ,t)}{\partial x_{j}}}=v_{Ii}~N_{I,j}}
The rate of deformation is then given by
d
h
(
x
,
t
)
=
1
2
(
l
h
+
l
h
T
)
=
∑
I
=
1
n
1
2
[
v
I
(
t
)
⊗
∇
N
I
(
X
)
+
∇
N
I
(
X
)
⊗
v
I
(
t
)
]
{\displaystyle {\boldsymbol {d}}_{h}(\mathbf {x} ,t)={\frac {1}{2}}~\left({\boldsymbol {l}}_{h}+{\boldsymbol {l}}_{h}^{T}\right)=\sum _{I=1}^{n}{\frac {1}{2}}\left[\mathbf {v} _{I}(t)\otimes {\boldsymbol {\nabla }}N_{I}(\mathbf {X} )+{\boldsymbol {\nabla }}N_{I}(\mathbf {X} )\otimes \mathbf {v} _{I}(t)\right]}
In index notation,
[
d
i
j
]
h
=
∑
I
=
1
n
1
2
[
v
I
i
N
I
,
j
+
v
I
j
N
I
,
i
]
{\displaystyle [d_{ij}]_{h}=\sum _{I=1}^{n}{\frac {1}{2}}[v_{Ii}~N_{I,j}+v_{Ij}~N_{I,i}]}
We can now substitute the trial and test functions into the weak form of the
momentum equation to get (using the same procedure as in one-dimension)
∫
Ω
ρ
N
I
a
h
dV
+
∫
Ω
σ
⋅
(
∇
N
I
)
dV
=
∫
Ω
ρ
N
I
b
dV
+
∫
∂
Ω
t
N
I
t
dA
{\displaystyle {\int _{\Omega }\rho ~N_{I}~\mathbf {a} _{h}~{\text{dV}}+\int _{\Omega }{\boldsymbol {\sigma }}\cdot ({\boldsymbol {\nabla }}N_{I})~{\text{dV}}=\int _{\Omega }\rho ~N_{I}~\mathbf {b} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}N_{I}~\mathbf {t} ~{\text{dA}}}}
We can write this approximate weak form as
f
I
kin
+
f
I
int
=
f
I
ext
{\displaystyle \mathbf {f} _{I}^{\text{kin}}+\mathbf {f} _{I}^{\text{int}}=\mathbf {f} _{I}^{\text{ext}}}
where the first term represents the inertial force, the second term the
internal force, and the third term the external force. Thus,
f
I
kin
:=
∫
Ω
ρ
N
I
a
h
dV
f
I
int
:=
∫
Ω
σ
⋅
(
∇
N
I
)
dV
f
I
ext
:=
∫
Ω
ρ
N
I
b
dV
+
∫
∂
Ω
t
N
I
t
dA
{\displaystyle {\begin{aligned}\mathbf {f} _{I}^{\text{kin}}&:=\int _{\Omega }\rho ~N_{I}~\mathbf {a} _{h}~{\text{dV}}\\\mathbf {f} _{I}^{\text{int}}&:=\int _{\Omega }{\boldsymbol {\sigma }}\cdot ({\boldsymbol {\nabla }}N_{I})~{\text{dV}}\\\mathbf {f} _{I}^{\text{ext}}&:=\int _{\Omega }\rho ~N_{I}~\mathbf {b} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}N_{I}~\mathbf {t} ~{\text{dA}}\end{aligned}}}
The term
f
I
kin
:=
∫
Ω
ρ
N
I
a
h
dV
{\displaystyle \mathbf {f} _{I}^{\text{kin}}:=\int _{\Omega }\rho ~N_{I}~\mathbf {a} _{h}~{\text{dV}}}
is the inertial force term.
If we substitute
a
(
X
,
t
)
=
∑
J
=
1
n
N
J
(
X
)
a
J
(
t
)
{\displaystyle \mathbf {a} (\mathbf {X} ,t)=\sum _{J=1}^{n}N_{J}(\mathbf {X} )~\mathbf {a} _{J}(t)}
into the inertial force terms, we get
f
I
kin
=
∑
J
=
1
n
[
∫
Ω
ρ
N
I
N
J
dV
]
a
J
{\displaystyle \mathbf {f} _{I}^{\text{kin}}=\sum _{J=1}^{n}\left[\int _{\Omega }\rho ~N_{I}~N_{J}~{\text{dV}}\right]~\mathbf {a} _{J}}
We define the mass matrix as the quantity
M
I
J
:=
∫
Ω
ρ
N
I
N
J
dV
=
∫
Ω
0
ρ
0
N
I
N
J
dV
0
{\displaystyle {M_{IJ}:=\int _{\Omega }\rho ~N_{I}~N_{J}~{\text{dV}}=\int _{\Omega _{0}}\rho _{0}~N_{I}~N_{J}~{\text{dV}}_{0}}}
Therefore, the inertial force can be expressed, in analogy with Newton's
second law, as
f
I
kin
=
∑
J
=
1
n
M
I
J
a
J
{\displaystyle \mathbf {f} _{I}^{\text{kin}}=\sum _{J=1}^{n}M_{IJ}~\mathbf {a} _{J}}
Semi-discrete finite element equations
edit
We can then write the semi-discrete version of the weak form in the shape
of a matrix equation as
[
M
]
[
a
]
=
[
f
ext
]
−
[
f
int
]
{\displaystyle [\mathbf {M} ][\mathbf {a} ]=[\mathbf {f} ^{\text{ext}}]-[\mathbf {f} ^{\text{int}}]}
These equations are ordinary differential equations because we still have
time derivatives on the left hand side - hence the system of equations is
semi-discrete.
There are a number of ways in which the discretization in time can be
effected. Some common approaches are the generalized midpoint rule, the
Newmark
β
{\displaystyle \beta }
method, or Runge-Kutta explicit integration schemes. There is a large literature on the best way of doing the discretization in time
with the goal of conserving both momentum and energy as accurately as
possible.
Taking derivatives and calculating integrals
edit
At this stage we are faced with the problem of calculating derivatives and
integrals with respect to
x
{\displaystyle \mathbf {x} }
of quantities that depend on
X
{\displaystyle \mathbf {X} }
. How
can we proceed?
The standard way of resolving this issue is to do all our calculations with
respect to a parent element and then map the results to the reference
or current configurations as necessary.
To make things more concrete let us consider a single element and label
the three domains as
Parent element (
◻
{\displaystyle \square }
).
Reference element (
Ω
0
{\displaystyle \Omega _{0}}
).
Current element (
Ω
{\displaystyle \Omega }
).
Let the maps that we need be
Parent
→
{\displaystyle \rightarrow }
Reference :
X
=
X
(
ξ
)
=
x
(
ξ
,
0
)
{\displaystyle \mathbf {X} =\mathbf {X} ({\boldsymbol {\xi }})=\mathbf {x} ({\boldsymbol {\xi }},0)}
.
Parent
→
{\displaystyle \rightarrow }
Current :
x
=
x
(
ξ
,
t
)
{\displaystyle \mathbf {x} =\mathbf {x} ({\boldsymbol {\xi }},t)}
.
Reference
→
{\displaystyle \rightarrow }
Current :
x
=
x
(
X
,
t
)
{\displaystyle \mathbf {x} =\mathbf {x} (\mathbf {X} ,t)}
.
This situation is illustrated (for a two dimensional situation) in
the following figure.
Maps from parent element to reference and current configurations
We define the map between the parent element and the
element in the current configuration as
x
(
ξ
,
t
)
=
∑
I
N
I
(
ξ
)
x
I
(
t
)
{\displaystyle \mathbf {x} ({\boldsymbol {\xi }},t)=\sum _{I}N_{I}({\boldsymbol {\xi }})~\mathbf {x} _{I}(t)}
Then,
v
(
ξ
,
t
)
=
∑
I
N
I
(
ξ
)
v
I
(
t
)
{\displaystyle \mathbf {v} ({\boldsymbol {\xi }},t)=\sum _{I}N_{I}({\boldsymbol {\xi }})~\mathbf {v} _{I}(t)}
Let us now try to compute the spatial velocity gradient. We have
∂
v
(
x
,
t
)
∂
ξ
=
∂
v
∂
x
⋅
∂
x
∂
ξ
=
∂
v
∂
x
⋅
F
ξ
{\displaystyle {\frac {\partial \mathbf {v} (\mathbf {x} ,t)}{\partial {\boldsymbol {\xi }}}}={\frac {\partial \mathbf {v} }{\partial \mathbf {x} }}\cdot {\frac {\partial \mathbf {x} }{\partial {\boldsymbol {\xi }}}}={\frac {\partial \mathbf {v} }{\partial \mathbf {x} }}\cdot {\boldsymbol {F}}_{\xi }}
The quantity
F
ξ
:=
∂
x
∂
ξ
=
∑
I
x
I
⊗
∂
N
I
∂
ξ
{\displaystyle {\boldsymbol {F}}_{\xi }:={\frac {\partial \mathbf {x} }{\partial {\boldsymbol {\xi }}}}=\sum _{I}\mathbf {x} _{I}\otimes {\frac {\partial N_{I}}{\partial {\boldsymbol {\xi }}}}}
is called the Jacobian of the motion and is a function of time.
Inverting, we get
∂
v
∂
x
=
∂
v
∂
ξ
⋅
F
ξ
−
1
{\displaystyle {{\frac {\partial \mathbf {v} }{\partial \mathbf {x} }}={\frac {\partial \mathbf {v} }{\partial {\boldsymbol {\xi }}}}\cdot {\boldsymbol {F}}_{\xi }^{-1}}}
But the spatial velocity gradient is
l
=
∂
v
∂
x
{\displaystyle {\boldsymbol {l}}={\frac {\partial \mathbf {v} }{\partial \mathbf {x} }}}
Therefore,
l
=
∂
v
∂
ξ
⋅
F
ξ
−
1
{\displaystyle {\boldsymbol {l}}={\frac {\partial \mathbf {v} }{\partial {\boldsymbol {\xi }}}}\cdot {\boldsymbol {F}}_{\xi }^{-1}}
If we substitute the approximation
v
(
ξ
,
t
)
=
∑
I
N
I
(
ξ
)
v
I
(
t
)
{\displaystyle \mathbf {v} ({\boldsymbol {\xi }},t)=\sum _{I}N_{I}({\boldsymbol {\xi }})~\mathbf {v} _{I}(t)}
we get a formula for the velocity gradient
l
=
[
∑
I
v
I
⊗
∂
N
I
∂
ξ
]
⋅
F
ξ
−
1
{\displaystyle {{\boldsymbol {l}}=\left[\sum _{I}\mathbf {v} _{I}\otimes {\frac {\partial N_{I}}{\partial {\boldsymbol {\xi }}}}\right]\cdot {\boldsymbol {F}}_{\xi }^{-1}}}
A similar procedure can be followed when other spatial gradients need to be
calculated. The most frequently encountered situation is the calculation of
the gradient of the shape functions. In that case, we have
∇
N
I
(
X
)
=
∂
N
I
(
X
)
∂
x
=
∂
N
I
∂
ξ
⋅
F
ξ
−
1
{\displaystyle {{\boldsymbol {\nabla }}N_{I}(\mathbf {X} )={\frac {\partial N_{I}(\mathbf {X} )}{\partial \mathbf {x} }}={\frac {\partial N_{I}}{\partial {\boldsymbol {\xi }}}}\cdot {\boldsymbol {F}}_{\xi }^{-1}}}
For computing integrals, we observe that for a function
f
(
x
)
{\displaystyle f(\mathbf {x} )}
on
Ω
{\displaystyle \Omega }
we can write
∫
Ω
f
(
x
)
dV
=
∫
Ω
0
f
^
(
X
)
J
dV
0
=
∫
◻
f
~
(
ξ
)
J
ξ
d
◻
{\displaystyle \int _{\Omega }f(\mathbf {x} )~{\text{dV}}=\int _{\Omega _{0}}{\hat {f}}(\mathbf {X} )~J~{\text{dV}}_{0}=\int _{\square }{\tilde {f}}({\boldsymbol {\xi }})~J_{\xi }~{\text{d}}\square }
where
J
ξ
=
det
(
∂
x
∂
ξ
)
=
det
(
F
ξ
)
{\displaystyle J_{\xi }=\det \left({\frac {\partial \mathbf {x} }{\partial {\boldsymbol {\xi }}}}\right)=\det({\boldsymbol {F}}_{\xi })}
Similarly for a function
g
(
X
)
{\displaystyle g(\mathbf {X} )}
over
Ω
0
{\displaystyle \Omega _{0}}
, we have
∫
Ω
0
g
(
X
)
dV
0
=
∫
◻
g
~
(
ξ
)
J
ξ
0
d
◻
{\displaystyle \int _{\Omega _{0}}g(\mathbf {X} )~{\text{dV}}_{0}=\int _{\square }{\tilde {g}}({\boldsymbol {\xi }})~J_{\xi }^{0}~{\text{d}}\square }
where
J
ξ
0
=
det
(
∂
X
∂
ξ
)
=
det
(
F
ξ
0
)
{\displaystyle J_{\xi }^{0}=\det \left({\frac {\partial \mathbf {X} }{\partial {\boldsymbol {\xi }}}}\right)=\det({\boldsymbol {F}}_{\xi }^{0})}
Recall that the internal force terms have the form
f
I
int
=
∫
Ω
σ
⋅
(
∇
N
I
)
dV
{\displaystyle \mathbf {f} _{I}^{\text{int}}=\int _{\Omega }{\boldsymbol {\sigma }}\cdot ({\boldsymbol {\nabla }}N_{I})~{\text{dV}}}
Expressed in terms of an integral over the parent element, we then have
f
I
int
=
∫
◻
σ
⋅
(
∇
N
I
)
J
ξ
d
◻
{\displaystyle {\mathbf {f} _{I}^{\text{int}}=\int _{\square }{\boldsymbol {\sigma }}\cdot ({\boldsymbol {\nabla }}N_{I})~J_{\xi }{\text{d}}\square }}
where
∇
N
I
=
∂
N
I
∂
ξ
⋅
F
ξ
−
1
.
{\displaystyle {\boldsymbol {\nabla }}N_{I}={\frac {\partial N_{I}}{\partial {\boldsymbol {\xi }}}}\cdot {\boldsymbol {F}}_{\xi }^{-1}~.}
Finite element implementation
edit
Recall that the finite element system of equations can be written as
[
M
]
[
a
]
=
[
f
ext
]
−
[
f
int
]
{\displaystyle [\mathbf {M} ][\mathbf {a} ]=[\mathbf {f} ^{\text{ext}}]-[\mathbf {f} ^{\text{int}}]}
where
M
I
J
=
∫
Ω
0
ρ
0
N
I
N
J
dV
0
;
f
I
int
=
∫
Ω
σ
⋅
(
∇
N
I
)
dV
;
and
f
I
ext
=
∫
Ω
ρ
N
I
b
dV
+
∫
∂
Ω
t
N
I
t
dA
{\displaystyle M_{IJ}=\int _{\Omega _{0}}\rho _{0}~N_{I}~N_{J}~{\text{dV}}_{0}~;~~\mathbf {f} _{I}^{\text{int}}=\int _{\Omega }{\boldsymbol {\sigma }}\cdot ({\boldsymbol {\nabla }}N_{I})~{\text{dV}}~;~~{\text{and}}~~\mathbf {f} _{I}^{\text{ext}}=\int _{\Omega }\rho ~N_{I}~\mathbf {b} ~{\text{dV}}+{\int _{\partial \Omega }}_{t}N_{I}~\mathbf {t} ~{\text{dA}}}
To get a feel for the structure of the equation, let us consider a four-noded
plane element. If the components of the velocity in the two coordinate
directions are
(
u
,
v
)
{\displaystyle (u,v)}
, the acceleration vector has the form
a
=
[
a
1
a
2
a
3
a
4
]
=
[
u
˙
1
v
˙
1
u
˙
2
v
˙
2
u
˙
3
v
˙
3
u
˙
4
v
˙
4
]
{\displaystyle \mathbf {a} ={\begin{bmatrix}\mathbf {a} _{1}\\\mathbf {a} _{2}\\\mathbf {a} _{3}\\\mathbf {a} _{4}\end{bmatrix}}={\begin{bmatrix}{\dot {u}}_{1}&{\dot {v}}_{1}\\{\dot {u}}_{2}&{\dot {v}}_{2}\\{\dot {u}}_{3}&{\dot {v}}_{3}\\{\dot {u}}_{4}&{\dot {v}}_{4}\end{bmatrix}}}
The mass matrix has the form
M
=
[
M
11
M
12
M
13
M
14
M
21
M
22
M
23
M
24
M
31
M
32
M
33
M
34
M
41
M
42
M
43
M
44
]
{\displaystyle \mathbf {M} ={\begin{bmatrix}\mathbf {M} _{11}&\mathbf {M} _{12}&\mathbf {M} _{13}&\mathbf {M} _{14}\\\mathbf {M} _{21}&\mathbf {M} _{22}&\mathbf {M} _{23}&\mathbf {M} _{24}\\\mathbf {M} _{31}&\mathbf {M} _{32}&\mathbf {M} _{33}&\mathbf {M} _{34}\\\mathbf {M} _{41}&\mathbf {M} _{42}&\mathbf {M} _{43}&\mathbf {M} _{44}\end{bmatrix}}}
In expanded form
M
=
∫
Ω
0
ρ
0
[
N
1
N
1
N
1
N
2
N
1
N
3
N
1
N
4
N
2
N
1
N
2
N
2
N
2
N
3
N
2
N
4
N
3
N
1
N
3
N
2
N
3
N
3
N
3
N
4
N
4
N
1
N
4
N
2
N
4
N
3
N
4
N
4
]
dV
{\displaystyle \mathbf {M} =\int _{\Omega _{0}}\rho _{0}~{\begin{bmatrix}N_{1}~N_{1}&N_{1}~N_{2}&N_{1}~N_{3}&N_{1}~N_{4}\\N_{2}~N_{1}&N_{2}~N_{2}&N_{2}~N_{3}&N_{2}~N_{4}\\N_{3}~N_{1}&N_{3}~N_{2}&N_{3}~N_{3}&N_{3}~N_{4}\\N_{4}~N_{1}&N_{4}~N_{2}&N_{4}~N_{3}&N_{4}~N_{4}\end{bmatrix}}~{\text{dV}}}
Next, we consider the internal force term. Let the components of the internal
force in the two coordinate directions be (
f
,
g
{\displaystyle f,g}
). Then,
f
int
=
[
f
1
int
f
2
int
f
3
int
f
4
int
]
=
[
f
1
g
1
f
2
g
2
f
3
g
3
f
4
g
4
]
{\displaystyle \mathbf {f} ^{\text{int}}={\begin{bmatrix}\mathbf {f} _{1}^{\text{int}}\\\mathbf {f} _{2}^{\text{int}}\\\mathbf {f} _{3}^{\text{int}}\\\mathbf {f} _{4}^{\text{int}}\end{bmatrix}}={\begin{bmatrix}f_{1}&g_{1}\\f_{2}&g_{2}\\f_{3}&g_{3}\\f_{4}&g_{4}\end{bmatrix}}}
At this stage, recall that the velocity gradient is given by
l
=
∇
v
=
∑
I
v
I
⊗
∇
N
I
=
∑
I
v
I
⊗
∂
N
I
∂
x
{\displaystyle {\boldsymbol {l}}={\boldsymbol {\nabla }}\mathbf {v} =\sum _{I}\mathbf {v} _{I}\otimes {\boldsymbol {\nabla }}N_{I}=\sum _{I}\mathbf {v} _{I}\otimes {\frac {\partial N_{I}}{\partial \mathbf {x} }}}
In index notation,
l
i
j
=
∑
I
v
I
i
∂
N
I
∂
x
j
.
{\displaystyle l_{ij}=\sum _{I}v_{Ii}~{\frac {\partial N_{I}}{\partial x_{j}}}~.}
We define
B
I
j
:=
∂
N
I
∂
x
j
{\displaystyle B_{Ij}:={\frac {\partial N_{I}}{\partial x_{j}}}}
Then,
l
i
j
=
∑
I
v
I
i
B
I
j
{\displaystyle l_{ij}=\sum _{I}v_{Ii}~B_{Ij}}
In matrix form
l
=
[
v
1
v
2
v
3
v
4
]
[
B
1
B
2
B
3
B
4
]
{\displaystyle \mathbf {l} ={\begin{bmatrix}\mathbf {v} _{1}&\mathbf {v} _{2}&\mathbf {v} _{3}&\mathbf {v} _{4}\end{bmatrix}}{\begin{bmatrix}\mathbf {B} _{1}\\\mathbf {B} _{2}\\\mathbf {B} _{3}\\\mathbf {B} _{4}\end{bmatrix}}}
Expanded out
[
l
11
l
12
l
21
l
22
]
=
[
u
1
u
2
u
3
u
4
v
1
v
2
v
3
v
4
]
[
∂
N
1
∂
x
∂
N
1
∂
y
∂
N
2
∂
x
∂
N
2
∂
y
∂
N
3
∂
x
∂
N
3
∂
y
∂
N
4
∂
x
∂
N
4
∂
y
]
{\displaystyle {\begin{bmatrix}l_{11}&l_{12}\\l_{21}&l_{22}\end{bmatrix}}={\begin{bmatrix}u_{1}&u_{2}&u_{3}&u_{4}\\v_{1}&v_{2}&v_{3}&v_{4}\end{bmatrix}}{\begin{bmatrix}{\frac {\partial N_{1}}{\partial x}}&{\frac {\partial N_{1}}{\partial y}}\\{\frac {\partial N_{2}}{\partial x}}&{\frac {\partial N_{2}}{\partial y}}\\{\frac {\partial N_{3}}{\partial x}}&{\frac {\partial N_{3}}{\partial y}}\\{\frac {\partial N_{4}}{\partial x}}&{\frac {\partial N_{4}}{\partial y}}\end{bmatrix}}}
In more compact form, we can just write
[
l
]
=
[
v
]
T
[
B
]
{\displaystyle {[\mathbf {l} ]=[\mathbf {v} ]^{T}~[\mathbf {B} ]}}
Similarly, for the internal force term,
f
I
int
=
∫
Ω
σ
⋅
∇
N
I
dV
{\displaystyle \mathbf {f} _{I}^{\text{int}}=\int _{\Omega }{\boldsymbol {\sigma }}\cdot {\boldsymbol {\nabla }}N_{I}~{\text{dV}}}
we have
f
I
i
int
=
∫
Ω
σ
i
j
∂
N
I
∂
x
j
dV
=
∫
Ω
σ
i
j
B
I
j
dV
{\displaystyle f_{Ii}^{\text{int}}=\int _{\Omega }\sigma _{ij}~{\frac {\partial N_{I}}{\partial x_{j}}}~{\text{dV}}=\int _{\Omega }\sigma _{ij}~B_{Ij}~{\text{dV}}}
In matrix form,
[
f
1
int
f
2
int
f
3
int
f
4
int
]
=
∫
Ω
[
B
1
B
2
B
3
B
4
]
σ
T
{\displaystyle {\begin{bmatrix}\mathbf {f} _{1}^{\text{int}}\\\mathbf {f} _{2}^{\text{int}}\\\mathbf {f} _{3}^{\text{int}}\\\mathbf {f} _{4}^{\text{int}}\end{bmatrix}}=\int _{\Omega }{\begin{bmatrix}\mathbf {B} _{1}\\\mathbf {B} _{2}\\\mathbf {B} _{3}\\\mathbf {B} _{4}\end{bmatrix}}~{\boldsymbol {\sigma }}^{T}}
Expanded out,
[
f
1
g
1
f
2
g
2
f
3
g
3
f
4
g
4
]
=
∫
Ω
[
∂
N
1
∂
x
∂
N
1
∂
y
∂
N
2
∂
x
∂
N
2
∂
y
∂
N
3
∂
x
∂
N
3
∂
y
∂
N
4
∂
x
∂
N
4
∂
y
]
[
σ
11
σ
21
σ
12
σ
22
]
dV
{\displaystyle {\begin{bmatrix}f_{1}&g_{1}\\f_{2}&g_{2}\\f_{3}&g_{3}\\f_{4}&g_{4}\end{bmatrix}}=\int _{\Omega }{\begin{bmatrix}{\frac {\partial N_{1}}{\partial x}}&{\frac {\partial N_{1}}{\partial y}}\\{\frac {\partial N_{2}}{\partial x}}&{\frac {\partial N_{2}}{\partial y}}\\{\frac {\partial N_{3}}{\partial x}}&{\frac {\partial N_{3}}{\partial y}}\\{\frac {\partial N_{4}}{\partial x}}&{\frac {\partial N_{4}}{\partial y}}\end{bmatrix}}{\begin{bmatrix}\sigma _{11}&\sigma _{21}\\\sigma _{12}&\sigma _{22}\end{bmatrix}}~{\text{dV}}}
In compact form, we have
[
f
int
]
=
∫
Ω
[
B
]
[
σ
]
T
dV
{\displaystyle {[\mathbf {f} ^{\text{int}}]=\int _{\Omega }[\mathbf {B} ][{\boldsymbol {\sigma }}]^{T}~{\text{dV}}}}
We can express the external force term in a similar manner. Notice that
we are not taking adavantage of the symmetry of the stress tensor in the above
expressions.
A widely used alternative way of expressing the finite element system of
equations is the Voigt notation. This notation can be found in discussions
of introductory finite element analysis.
Algorithm for computing internal forces
edit
The major complications and computational effort in the finite element
implementation are encountered during the computation of internal forces.
The basic procedure for computing the internal force at the nodes of
an element is given below.
Set
[
f
int
]
=
[
0
]
{\displaystyle [\mathbf {f} ^{\text{int}}]=[\mathbf {0} ]}
.
For all quadrature points in the parent element (
ξ
Q
{\displaystyle {\boldsymbol {\xi }}_{Q}}
)
Compute the gradient of the shape functions (
[
B
I
]
{\displaystyle [\mathbf {B} _{I}]}
) for all nodes
I
{\displaystyle I}
.
Compute the velocity gradient
[
l
]
=
[
v
I
]
[
B
]
{\displaystyle [\mathbf {l} ]=[\mathbf {v} _{I}][\mathbf {B} ]}
.
Compute the rate of deformation tensor.
Compute the deformation gradient/ the Lagrangian Green tensor.
Compute the Cauchy stress or the 2nd P-K stress using the constitutive equation.
Update the nodal internal force (
[
f
I
int
]
{\displaystyle [\mathbf {f} _{I}^{\text{int}}]}
) using
f
I
int
=
f
I
int
+
[
B
I
]
[
σ
]
T
J
ξ
W
Q
{\displaystyle \mathbf {f} _{I}^{\text{int}}=\mathbf {f} _{I}^{\text{int}}+[\mathbf {B} _{I}][{\boldsymbol {\sigma }}]^{T}~J_{\xi }~W_{Q}}
where
W
Q
{\displaystyle W_{Q}}
are the weights for Gaussian integration.
The mass matrix and integration in time will be discussed later.