Idea Transcript
WEIGHTED RESIDUAL METHOD
1
INTRODUCTION • Direct stiffness method is limited for simple 1D problems • PMPE is limited to potential problems • FEM can be applied to many engineering problems that are governed by a differential equation • Need systematic approaches to generate FE equations – Weighted residual method – Energy method
• Ordinary differential equation (second-order (second order or fourth-order) fourth order) can be solved using the weighted residual method, in particular using Galerkin method
2
EXACT VS. APPROXIMATE SOLUTION • Exact solution – Boundary value problem: differential equation + boundary conditions – Displacements in a uniaxial bar subject to a distributed force p(x) d 2u + p( x ) = 0, 0 ≤ x ≤ 1 dx 2 u(0) = 0 ⎫ ⎪ ⎪ ⎪ Boundary conditions ⎬ du (1) = 1⎪ ⎪ ⎪ dx ⎭ – Essential BC: The solution value at a point is prescribed (displacement or kinematic BC) – Natural BC: The derivative is given at a point (stress BC) – Exact E t solution l ti u(x): ( ) twice t i differential diff ti l function f ti – In general, it is difficult to find the exact solution when the domain and/or boundary conditions are complicated – Sometimes the solution may not exists even if the problem is well defined 3
EXACT VS. APPROXIMATE SOLUTION cont. • Approximate solution – It satisfies the essential BC, but not natural BC – The approximate solution may not satisfy the DE exactly – Residual: d 2u + p( x ) = R ( x ) dx 2 – Want to minimize the residual by y multiplying p y g with a weight g W and integrate over the domain 1
∫0 R( x )W ( x )dx = 0
Weight function
– If it satisfies for any W(x), then R(x) will approaches zero, and the approximate solution will approach the exact solution – Depending on choice of W(x): least square error method, method collocation method, Petrov-Galerkin method, and Galerkin method
4
GALERKIN METHOD • Approximate solution is a linear combination of trial functions N
u( x ) = ∑ ci φi ( x )
Trial function
i =1
– Accuracy depends on the choice of trial functions – The approximate solution must satisfy the essential BC
• Galerkin G l ki method th d – Use N trial functions for weight functions 1
0 ∫0 R( x )φi ( x )ddx = 0,
i =1 1,…, N
1⎛ d 2u
⎞⎟ ⎜⎜ + p ( x ) ⎟⎟φi ( x )dx = 0, i = 1,…, N ∫0 ⎜⎝ dx 2 ⎠ 1 d 2u
1
∫0 dx 2 φi ( x )dx = −∫0 p( x )φi ( x )dx,
i = 1, 1 …, N 5
GALERKIN METHOD cont. • Galerkin method cont cont. – Integration-by-parts: reduce the order of differentiation in u(x)
du φ dx i
1 0
−∫
1 du d dφ φi
0
dx dx
1
dx = −∫ p( x )φi ( x )dx, i = 1,…,N 0
– Apply Appl natural nat ral BC and rearrange 1 dφ i
∫0
du dx = dx d d dx
du
1
du
∫0 p( x )φi ( x )dx + ddx ((1))φi ((1)) − ddx ((0))φi ((0),),
i = 1,,…, N
– Same order of differentiation for both trial function and approx. solution – Substitute S b tit t the th approximate i t solution l ti 1 dφ i
∫0
N
c dx ∑ j j =1
dφ j dx = dx
1
du
du
∫0 p( x )φi ( x )dx + dx (1)φi (1) − dx (0)φi (0),
i = 1,…, N
6
GALERKIN METHOD cont. • Galerkin method cont cont. – Write in matrix form N
∑ Kij c j = Fi ,
[K] {c} = {F}
i = 1,…, N
( N×N )( N×1)
j =1
K ij =
Fi =
1
∫0
( N×1)
d φi d φ j dx dx dx
du
1
du
((1))φi ((1)) − ((0))φi (0) ( ) ∫0 p( x )φi ( x )dx + dx d d dx
– Coefficient matrix is symmetric; Kij = Kji – N equations with N unknown coefficients
7
EXAMPLE1 • Differential equation
Trial functions
2
d u + 1 = 0,0 ≤ x ≤ 1 dx 2 u(0) = 0 ⎪⎫ ⎪ ⎪ ⎬ Boundary conditions du (1) = 1⎪ ⎪ ⎪ dx ⎭
φ1( x ) = x
φ1′( x ) = 1
φ2 ( x ) = x 2
φ2′ ( x ) = 2 x
• Approximate solution (satisfies the essential BC) 2
u( x ) = ∑ ci φi ( x ) = c1x + c2 x 2 i=1
• Coefficient matrix and RHS vector 1
K11 = ∫ ( φ1′ ) dx = 1 2
0
1
K12 = K21 = ∫ ( φ1′φ2′ ) dx = 1 0
1
2 K22 = ∫ ( φ2′ ) dx = 0
4 3
1
du
3
F1 =
∫0 φ1( x )dx + φ1(1) − dx (0)φ1(0) = 2
F2 =
∫0 φ2 ( x )dx + φ2 (1) − dx (0)φ2 (0) = 3
1
du
4
8
EXAMPLE1 cont. • Matrix equation 1 ⎡3 3⎤ ⎥ [K ] = ⎢ 3 ⎢⎣ 3 4 ⎥⎦
{F} =
⎧ ⎪ 2 ⎫ ⎪ {c } = [K ]− 1 {F } = ⎪ ⎨ 1⎪ ⎬ ⎪ ⎪− 2 ⎪ ⎪ ⎩ ⎭
1 ⎪⎧ 9 ⎪⎫ ⎨ ⎬ 6 ⎪⎩⎪ 8 ⎪⎭⎪
• Approximate solution u( x ) = 2 x −
x2 2
– Approximate solution is also the exact solution because the linear combination of the trial functions can represent the exact solution
9
EXAMPLE2 • Differential equation
Trial functions
2
d u + x = 0, 0 ≤ x ≤1 dx 2 u(0) = 0 ⎪⎫ ⎪ ⎪ ⎬ Boundary conditions du (1) = 1⎪ ⎪ ⎪ dx ⎭
φ1( x ) = x
φ1′( x ) = 1
φ2 ( x ) = x 2
φ2′ ( x ) = 2 x
1 ⎧ 16 ⎫
⎪ ⎪ • Coefficient matrix is same, force vector: {F} = 12 ⎨⎪15 ⎬⎪ ⎩⎪
{c } = [K ]− 1 {F } =
⎧⎪ 19 ⎫⎪ ⎪⎨ 12 ⎪⎬ ⎪⎩⎪ − 41 ⎪⎭⎪
u( x ) =
⎭⎪
19 x2 x− 12 4
• Exact solution u( x ) =
3 x3 x− 2 6
– The trial functions cannot express p the exact solution;; thus,, approximate solution is different from the exact one 10
EXAMPLE2 cont. • Approximation is good for u(x), u(x) but not good for du/dx 1.6 1.4
u(x), du/dx u
1.2 1.0 0.8 06 0.6 0.4 0.2 0.0 0
0.2
u-exactt
u-approx.
0.4
x
0.6
d /d ((exact) du/dx t)
0.8
1
d /d ((approx.)) du/dx
11
HIGHER-ORDER DIFFERENTIAL EQUATIONS • Fourth-order Fourth order differential equation d 4w − p( x ) = 0, dx 4
0≤x≤L
– Beam bending under pressure load
• Approximate solution N
w (0) = 0 ⎫⎪ ⎪⎪ ⎬ Essential BC dw (0) = 0 ⎪⎪ ⎪⎭ dx 2 ⎫⎪ d w (L ) = M ⎪⎪ 2 ⎪⎪ dx ⎬ Natural BC ⎪⎪ d 3w L V ( ) = − ⎪ dx 3 ⎭⎪⎪
w ( x ) = ∑ ci φi ( x ) i =1
• Weighted W i ht d residual id l equation ti (G (Galerkin l ki method) th d) L ⎛ d 4w
⎞
∫0 ⎝⎜⎜⎜ dx 4 − p( x ) ⎠⎟⎟⎟φi ( x )dx = 0,
i = 1,…, N
– In order to make the order of differentiation same, integration-by-parts must be done twice
12
HIGHER-ORDER DE cont. • After integration-by-parts integration by parts twice d 3w φi dx 3
L
0
d 2w d φi − 2 dx dx
d 2φi dx = dx 2 dx 2
L d 2w
∫0
L
+∫
d 2φi dx = dx 2 dx 2
L d 2w
0
0
L
∫0
L
∫0
d 3w p( x )φi ( x )dx − 3 φi dx
L
0
p( x )φi ( x )dx, i = 1,…, N
d 2w d φi + 2 dx dx
L
, i = 1,…, N 0
• Substitute approximate solution L
∫0
d 2φ j d 2φi dx = ∑ c j dx d 2 d dx 2 j =1 N
L
∫0
d 3w p( x )φi ( x )dx − 3 φi d dx
L
0
d 2w d φi + 2 d d dx dx
L
, i = 1,…, N 0
– Do not substitute the approx. solution in the boundary terms
• Matrix form [K] {c} = {F}
N×N N×1
L d 2φ i 2 0
Kij = ∫
N×1
Fi = ∫
L
0
d 2φj
dx
dx 2
dx
d 3w p( x )φi ( x )dx − 3 φi dx
L
0
d 2w dφi + 2 dx dx
L
0
13
EXMAPLE dw (0) = 0 dx d 2w d 3w = (1) 2 (1) = −1 dx 2 dx 3
• Fourth-order Fourth order DE d 4w − 1 = 0, dx 4
w (0) = 0
0≤x≤L
• Two trial functions φ1 = x 2, φ2 = x 3
φ1′′ = 2, 2 φ2′′ = 6x
• Coefficient matrix K11 =
1
∫0 ( φ1′′)
K12 = K 21 = K 22 =
1
2
dx = 4 1
∫0 ( φ1′′φ2′′ ) dx = 6
∫0 ( φ2′′ )
2
⎡4 6 ⎤ ⎥ [K ] = ⎢ ⎢⎣ 6 12 ⎥⎦
dx = 12
14
EXAMPLE cont. • RHS F1 = F2 =
1
∫0
d 3w (0) d 2w (0) 16 ′ ′(0) = (0) M (1) φ + φ − φ 1 1 1 3 dx 3 dx 2
x 2dx + V φ1(1) +
1
∫0
x 3dx + V φ2 (1) +
d 3w (0) d 2w (0) 29 ′ (0) M (1) φ + φ − φ2′ (0) = 2 2 3 2 4 dx dx
• Approximate solution ⎧ 41 ⎪ ⎫ ⎪ {c} = [K ]−1 {F} = ⎪ ⎨ 241 ⎪ ⎬ ⎪⎩ ⎪ − 4 ⎪⎭ ⎪
w(x) =
41 2 1 3 x − x 24 4
• Exact solution w(x) =
1 4 1 3 7 2 x − x + x 24 3 4
15
EXAMPLE cont. 4 3
w'', w''
2 1 0 -1 -2 -3 3 0 w'' (exact)
0.2
0.4 w'' (approx.)
x
0.6 w''' (exact)
0.8
1 w''' (approx.)
16
FINITE ELEMENT APPROXIMATION • Domain Discretization – Weighted residual method is still difficult to obtain the trial functions that satisfy the essential BC – FEM is i to t divide di id th the entire ti d domain i iinto t a sett off simple i l sub-domains bd i (finite element) and share nodes with adjacent elements – Within a finite element, the solution is approximated in a simple polynomial l i l form f u(x)
Approximate solution x Finite elements
Analytical solution
– Wh When more number b off finite fi it elements l t are used, d th the approximated i t d piecewise linear solution may converge to the analytical solution 17
FINITE ELEMENT METHOD cont. • Types of finite elements
1D
2D
3D
• Variational equation is imposed on each element. 1
∫0
dx =
0.1
∫0
dx + ∫
0.2
01 0.1
dx +
+∫
1
0 0.9 9
dx
One element
18
TRIAL SOLUTION – Solution within an element is approximated using simple polynomials polynomials. 1 1
2 2
n
n−1 3
n−1
xi
n
n+1
xi+1 i
– i-th element is composed of two nodes: xi and xi+1. Since two unknowns are involved,, linear polynomial p y can be used: u ( x ) = a0 + a1x,
xi ≤ x ≤ xi +1
– The unknown coefficients, coefficients a0 and a1, will be expressed in terms of nodal solutions u(xi) and u(xi+1).
19
TRIAL SOLUTION cont. – Substitute two nodal values
⎧ u( xi ) = ui = a0 + a1xi ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u( xi +1) = ui +1 = a0 + a1xi +1 – Express a0 and a1 in terms of ui and ui+1. Then, the solution is approximated by
u( x ) =
xi +1 − x x−x ui + ( i ) i ui +1 (i ) L L Ni ( x )
Ni +1( x )
– Solution for i-th element: u ( x ) = Ni ( x )ui + Ni +1( x )ui +1,
– Ni(x) and Ni+1(x): Shape
xi ≤ x ≤ xi +1
Function or Interpolation Function 20
TRIAL SOLUTION cont. • Observations – Solution u(x) is interpolated using its nodal values ui and ui+1. – Ni(x) = 1 at node xi, and =0 at node xi+1. Ni(x)
Ni+1(x)
xi
xi+1
– The solution is approximated by piecewise linear polynomial and its gradient is constant within an element element. u
ui+2 ui
xi
ui+1 xi+1
du dx
xi
xi+2
xi+1
xi+2
– Stress and strain (derivative) are often averaged at the node. 21
GALERKIN METHOD • Relation between interpolation functions and trial functions – 1D problem with linear interpolation
ND
u( x ) = ∑ ui φi ( x ) i =1
⎧ 0, 0 ≤ x ≤ xi −1 ⎪ ⎪ ⎪ ⎪ ( i −1) x−x ⎪ Ni ( x ) = ( i −1)i −1 , xi −1 < x ≤ xi ⎪ ⎪ L φi ( x ) = ⎪ ⎨ ⎪ xi +1 − x (i ) ⎪ , xi < x ≤ xi +1 ⎪⎪ Ni ( x ) = L( i ) ⎪ ⎪ ⎪ xi +1 < x ≤ xND ⎪ ⎩ 0,
– Difference: the interpolation function does not exist in the entire domain but it exists only in elements connected to the node domain,
• Derivative ⎧⎪ 0, ⎪⎪ ⎪⎪ 1 , ⎪ d φi ( x ) ⎪⎪ L( i −1) =⎨ dx ⎪⎪ 1 ⎪⎪ − ( i ) , ⎪⎪ L ⎪⎪⎩ 0,
0 ≤ x ≤ xi −1
1 φi (x )
( i −1 )
xi −1 < x ≤ xi 1/ L
x i −2
xi < x ≤ xi +1 −1/ L( i ) xi +1 < x ≤ xND
x i −1
xi
x i +1
d φi (x ) dx 22
EXAMPLE • Solve using two equal equal-length length elements u (0) = 0 ⎫ ⎪ ⎪ ⎪ ⎬ Boundary conditions du (1) = 1⎪ ⎪ ⎪ dx ⎭
d 2u + 1 = 0,0 ≤ x ≤ 1 dx 2
• Three nodes at x = 0, 0.5, 1.0; displ at nodes = u1, u2, u3 • Approximate solution u( x ) = u1φ1( x ) + u2φ2 ( x ) + u3φ3 ( x ) ⎧1 − 2 x, 0 ≤ x ≤ 0.5 φ1( x ) = ⎪ ⎨ ⎪ 0.5 < x ≤ 1 ⎪ ⎩ 0, ⎧ 0 ≤ x ≤ 0.5 ⎪ 0, φ3 ( x ) = ⎨ 0.5 < x ≤ 1 ⎪ ⎩⎪ −1 + 2 x,
⎧ 2 x, 0 ≤ x ≤ 0.5 φ2 ( x ) = ⎪⎨ ⎪⎪⎩ 2 − 2 x, 0.5 < x ≤ 1 1
φ_1 φ_2 φ_3
φ 0.5
0 0
0.5
x
1
23
EXAMPLE cont. • Derivatives of interpolation functions 0 ≤ x ≤ 0.5 d φ1( x ) ⎧ ⎪ −2, =⎨ dx ⎪ ⎪⎩ 0, 0.5 < x ≤ 1 ⎧ 0, 0 ≤ x ≤ 0.5 d φ3 ( x ) ⎪ =⎨ dx ⎪ ⎪ ⎩ 2, 0.5 < x ≤ 1
d φ2 ( x ) ⎧⎪ 2, 0 ≤ x ≤ 0.5 =⎨ dx ⎪⎪⎩ −2, 0.5 < x ≤ 1
• Coefficient matrix 0.5 1 d φ1 d φ2 dx ( 2)(2) )( ) dx ((0)( )(−2))dx = −2 = − + ∫0 dx dx ∫0 ∫0.5 05 1 dφ dφ 0.5 1 2 2 =∫ dx = ∫ 4dx + ∫ 4dx = 4 0 dx dx 0 0.5 1
K12 = K 22
• RHS 0.5
F1 =
∫0
F2 =
∫0
1
1× (1− 2 x )dx + ∫ 1× (0)dx + 0.5
0.5
2 xdx + ∫
1
0.5
(2 − 2 x )dx +
du du du (1)φ1(1) − (0)φ1(0) = 0.25 − (0) dx dx dx
du du (1)φ2 (1) − (0)φ2 (0) = 0.5 dx dx 24
EXAMPLE cont. • Matrix equation ⎡ 2 −2 0 ⎤ ⎧⎪ u1 ⎫⎪ ⎧⎪ F1 ⎫⎪ ⎪⎪ ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎪ ⎢ −2 4 −2 ⎥ ⎪ ⎨ u2 ⎬ = ⎨ 0.5 ⎬ ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎢⎣ 0 −2 2 ⎥⎦ ⎩⎪ ⎪ u3 ⎭⎪ ⎪ ⎩⎪ ⎪1.25 ⎭⎪⎪
Consider it as unknown
• Striking the 1st row and striking the 1st column (BC) ⎡ 4 −2 ⎤ ⎧⎪ u2 ⎫⎪ ⎧⎪ 0.5 ⎫⎪ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ ⎢⎣ −2 2 ⎥⎦ ⎪⎩⎪ u3 ⎪⎪⎭ ⎩⎪⎪1.25 ⎭⎪⎪
• Solve for u2 = 0.875, u3 = 1.5 • Approximate solution ⎧1.75 1 75 x, ⎪ u( x ) = ⎨ ⎪ ⎪ ⎩ 0.25 + 1.25 x,
0≤x ≤0 0.5 5 0.5 ≤ x ≤ 1
– Piecewise linear solution
25
EXAMPLE cont. 16 1.6
u(x)
1.2
0.8 u-exact
0.4
u-approx.
0 0
0.2
0.4
x
0.6
0.8
1
2
1.5 du/dx x
• Solution comparison • Approx. solution has about 8% error • Derivative shows a large discrepancy • Approx. A derivative d i ti iis constant as the solution is piecewise linear p
1 du/dx (exact)
0.5
du/dx (approx.)
0 0
0.2
0.4
x
0.6
0.8
1
26
FORMAL PROCEDURE • Galerkin method is still not general enough for computer code • Apply Galerkin method to one element (e) at a time • Introduce a local coordinate ξ=
x = xi (1 − ξ ) + x j ξ
x − xi x−x = (e ) i x j − xi L
• Approximate solution within the element N 1(ξ) N 2 (ξ)
u( x ) = ui N1( x ) + u j N2 ( x )
N1(ξ ) = (1− ξ )
Element e
N2 (ξ ) = ξ
xi
⎛ x−x N1( x ) = ⎜⎜1 − ( e ) i ⎝ L N2 ( x ) =
⎞⎟ ⎠⎟⎟
x − xi L( e )
dN1 dN1 d ξ 1 = = − (e ) dx d ξ dx L dN2 dN2 d ξ 1 = = + (e ) dx d ξ dx L
ξ
xj
L(e )
27
FORMAL PROCEDURE cont. • Interpolation property N1( xi ) = 1, N1( x j ) = 0 N2 ( xi ) = 0,, N2 ( x j ) = 1
u( xi ) = ui u( x j ) = u j
• Derivative of approx. solution du dN1 dN2 + uj = ui dx dx dx ⎧ u1 ⎪⎫ 1 ⎢ dN du ⎢ dN1 dN2 ⎥ ⎪ =⎢ ⎥ ⎨ ⎬ = (e ) ⎢ 1 dx ⎣⎢ dx dx ⎦⎥ ⎪ ⎪ u2 ⎪⎭⎪ L ⎢⎣ d ξ ⎩
⎧ u1 ⎪ ⎫ dN2 ⎥ ⎪ ⎥⎨ ⎬ d ξ ⎦⎥ ⎪ ⎪ u2 ⎪ ⎪ ⎩ ⎭
• Apply Galerkin method in the element level xj
∫x
i
dNi du dx = dx dx
xj
∫x
i
p( x )Ni ( x )dx +
du du ( x j )Ni ( x j ) − ( xi )Ni ( xi ), i = 1,2 dx dx
28
FORMAL PROCEDURE cont. • Change variable from x to ξ ⎢ dN1 ⎢ ⎣⎢ d ξ
1 1 dNi L( e ) ∫0 d ξ +
1 ⎧ u1 ⎫ dN2 ⎥ ⎥ d ξ i⎪⎨ ⎪⎬ = L( e ) ∫ p( x )Ni (ξ )d ξ 0 d ξ ⎦⎥ ⎪⎩⎪ u2 ⎪⎭⎪
du du ( x j )Ni (1) − ( x )N (0), i = 1,2 dx dx i i
– Do not use approximate solution for boundary terms
• Element-level matrix equation ⎧ du ⎫ d ⎪ ⎪ ⎪ − ( xi ) ⎪ ⎪ ⎪ dx (e ) (e ) (e ) ⎪ ⎪ [k ]{u } = { f } + ⎨ ⎬ ⎪ ⎪ du ⎪ + ( x j )⎪ ⎪ ⎪ ⎪ dx ⎪ d ⎩ ⎭
1 ⎧ ⎪ N (ξ ) ⎫⎪ {f (e ) } = L(e ) ∫ p( x )⎨ 1 ⎬ d ξ 0 ⎪ ⎪ N2 (ξ ) ⎪⎭⎪ ⎩
⎡ ⎛ dN ⎞2 ⎢ ⎜ 1 ⎟⎟ 1⎢ ⎜ ⎝ d ξ ⎟⎠ 1 ⎡ k( e ) ⎤ = ⎢ ⎣ ⎦ L( e ) ∫0 ⎢ 2×2 ⎢ dN2 dN1 ⎢ dξ dξ ⎣
dN1 dN2 ⎤⎥ ⎡ 1 −1⎤ dξ dξ ⎥ ⎥ dξ = 1 ⎢ ⎥ (e ) ⎢ 2 ⎥ ⎥ 1 1 − L ⎣ ⎦ ⎛ dN2 ⎞⎟ ⎥ ⎜ ⎜⎝ d ξ ⎟⎟⎠ ⎥⎦ 29
FORMAL PROCEDURE cont. • Need to derive the element element-level level equation for all elements • Consider Elements 1 and 2 (connected at Node 2) ⎧ ⎫ du ⎪ ⎪ (1) ⎪ ⎪ − ( x ) 1 ⎡ k11 k12 ⎤ (1) ⎧ ⎪ ⎪ ⎫ ⎧ ⎫ u1 ⎪ ⎪ f1 ⎪ ⎪ dx ⎪⎬ ⎢ ⎥ ⎨ ⎬ = ⎨ ⎬ + ⎪⎨ ⎢⎣ k21 k22 ⎥⎦ ⎪ ⎪⎪ du ⎪ ⎪ u2 ⎪ ⎪ ⎪ f2 ⎪ ⎪ ⎩ ⎭ ⎪ ⎩ ⎭ ⎪⎪ + ( x2 ) ⎪⎪⎪ ⎩ dx ⎭ ⎧ ⎫ du ⎪ ⎪ (2) ⎪ ⎪ − ( x ) 2 ⎡ k11 k12 ⎤ (2) ⎧ ⎪ ⎪ ⎫ ⎧ ⎫ u2 ⎪ ⎪ f2 ⎪ ⎪ dx ⎪⎬ ⎢ ⎥ ⎨ ⎬ = ⎨ ⎬ + ⎪⎨ ⎢⎣ k21 k22 ⎦⎥ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ u f du ⎪ 3⎭ ⎪ ⎩ ⎪ 3⎭ ⎪ ⎪⎪ + ( x ) ⎪⎪⎪ ⎪ dx 3 ⎭ ⎪ ⎩
• Assembly ⎡ k (1) ⎢ 11 ⎢ k (1) ⎢ 21 ⎢ ⎢ ⎣ 0
(1) k12 (1) (2) k22 + k11 (2) k21
Vanished ⎧ ⎫ du ⎪ ⎪ (1) ⎪⎪⎫ ⎪⎪ − dx ( x1 ) ⎪⎪ unknown term 0 ⎤⎥ ⎪ ⎧ u1 ⎪⎫ ⎧⎪⎪ f1 ⎪⎪ ⎪ ⎪⎪ ⎪⎪ (1) ⎪ ⎪⎪⎪ (2) ⎥ ⎪ (2) ⎪ k12 u f f 0 = + + ⎨2 ⎨ ⎬⎪ 2⎬ ⎥⎨ 2 ⎬ ⎪ ⎪⎪ ⎪⎪ ⎪ ⎪ ⎪⎪⎪ (2) ⎥ ⎪ d u3 ⎪⎭ ⎪ f (2) ⎪⎪ ⎪⎪ du ⎥⎪ k22 ⎩ 3 ⎪⎭ ⎪⎪ ( x3 ) ⎪⎪ ⎦ ⎩⎪ ⎪⎩ dx ⎪⎭ 30
FORMAL PROCEDURE cont. • Assembly of NE elements (ND = NE + 1) ⎡ (1) ⎢ k11 ⎢ ⎢ k ((1)) ⎢ 21 ⎢ ⎢ 0 ⎢ ⎢ ⎢ ⎢ ⎢⎣ 0
(1) k12
0
(1) ( ) (2) ( ) k22 + k11
((2)) k12
(2) k221
(2) (2) + k11 k22
0
0
…
( NE ) k21
( ND×ND )
⎤ ⎧ du ⎫ ⎪ − ( x1 ) ⎪⎪ ⎪ ⎧⎪ f (1) ⎫ ⎪⎪ ⎪ 0 ⎥⎧ ⎫ 1 ⎪ ⎪ u ⎪ ⎪ dx ⎥⎪ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ f ((1)) + f ((2)) ⎪⎪ ⎪ ⎪ 0 ⎥⎥ ⎪ 0 u2 ⎪ ⎪ 2 2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎥ ⎪ ⎪ = ⎪ (2) ⎪ (3) ⎪+ ⎪ 0 ⎬ 0 ⎥ ⎨ u3 ⎬ ⎨ f3 + f3 ⎬ ⎨ ⎪ ⎪ ⎪ ⎪ ⎥⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎥⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎥⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ( NE ) ( NE ) ⎥ ⎪ du u ⎪ N⎭ ⎪ ⎪ ⎪ ⎪ f k22 ⎥ ⎩ ⎪ N ⎪ ⎪ + ( xN ) ⎪⎪ ⎩ ⎭ ⎪ dx ⎪ ( ND×1) ⎩ ⎭ ⎦ ( ND×1) ( ND×1)
[K ]{q} = {F}
• Coefficient matrix [K] is singular; it will become non-singular after applying boundary conditions
31
EXAMPLE • Use three equal equal-length length elements d 2u + x = 0, dx 2
0 ≤ x ≤1
u(0) = 0, u(1) = 0
• All elements have the same coefficient matrix 1 ⎡ 1 −1⎤ ⎡ 3 −3 ⎤ ⎡ k(e ) ⎤ = 1 2 3) ⎣ ⎦ 2×2 L( e ) ⎢⎢ −1 1 ⎥⎥ = ⎢⎢ −3 3 ⎥⎥ , (e = 1,2,3) ⎣ ⎦ ⎣ ⎦
• Change variable of p(x) = x to p(ξ): • RHS
p(ξ ) = xi (1 − ξ ) + x j ξ
1 1 ⎧⎪ N (ξ ) ⎫ ⎧ ⎪ ⎪1 − ξ ⎫ ⎪ {f (e ) } = L(e ) ∫ p( x ) ⎨ 1 ⎬ d ξ = L( e ) ∫ [ xi ( 1 − ξ ) + x j ξ ]⎨ ⎬dξ 0 0 ⎪⎪⎩ N2 (ξ ) ⎪ ⎪⎩ ⎪ ξ ⎪⎭ ⎪ ⎪⎭ ⎧ xi ⎫ xj ⎪ ⎪ ⎪ + ⎪ ⎪ 3 6⎪ ⎪ (e = 1,2,3) , , ) = L(e ) ⎨⎪ ⎬, ⎪⎪ xi ⎪ xj ⎪ ⎪⎪ + ⎪⎪ 3⎭ ⎩6 32
EXAMPLE cont. • RHS cont cont.
⎪⎪⎧ f1(1) ⎪⎪⎫ ⎫ 1 ⎪⎧ 1 ⎪ ⎨ (1) ⎬ = ⎨ ⎬, ⎪ ⎪ ⎪ f2 ⎪⎭ ⎪ 54 ⎪⎩⎪ 2 ⎪⎭ ⎩
• Assembly ⎡ 3 −3 ⎢ ⎢ −3 3 + 3 ⎢ ⎢ 0 −3 ⎢ ⎢⎣ 0 0
0 −3 3+3 −3
⎧ ⎪ f2(2) ⎪⎫ ⎪ ⎧4⎪ ⎫ 1⎪ ⎪ ⎨ (2) ⎬ = ⎨ ⎬, ⎪ 5 ⎪⎭ ⎪ ⎪⎪⎩ f3 ⎪⎪⎭ 54 ⎪⎩
⎪⎧ f3(3) ⎪⎪⎫ ⎧7⎪ ⎫ 1⎪ ⎪ ⎨ (3) ⎬ = ⎨ ⎬ ⎪ 8 ⎪⎭ ⎪ ⎪⎪⎩ f4 ⎪⎪⎭ 54 ⎪⎩
⎧⎪ 1 du ⎫⎪ ⎪⎪ − (0) ⎪⎪ ⎪⎪ 54 dx ⎪⎪ 0 ⎤ ⎪⎧ u1 ⎪⎫ ⎪ 2 ⎪ 4 ⎪ ⎥ ⎪⎪ ⎪⎪ ⎪ 0 ⎥ ⎪⎪ u2 ⎪⎪ ⎪⎪⎪ 54 + 54 ⎪⎪⎪ ⎥⎨ ⎬ = ⎨ ⎬ ⎪⎪ −3 ⎥ ⎪⎪ u3 ⎪⎪ ⎪⎪ 7 5 ⎥⎪ ⎪ ⎪ + ⎪ 3 ⎥⎦ ⎩⎪⎪ u4 ⎭⎪⎪ ⎪⎪ 54 54 ⎪⎪ ⎪⎪ 8 du ⎪⎪⎪ ⎪⎪ (1) + ⎩⎪ 54 dx ⎪⎭⎪
Element 1 Element 2 Element 3
• Apply boundary conditions
– Deleting 1st and 4th rows and columns ⎡ 6 −3 ⎤ ⎧ ⎪ u ⎫ 1 ⎧⎪ 1 ⎫ ⎢ ⎥⎨ 2 ⎪ ⎬= ⎨ ⎪ ⎬ ⎢⎣ −3 6 ⎥⎦ ⎩⎪ ⎪ u3 ⎭⎪ ⎪ 9 ⎪⎩⎪ 2 ⎪⎭⎪
u2 =
4 81
u3 =
5 81
33
EXAMPLE cont. • Approximate solution 1 3 1 2 ≤x≤ 3 3
0≤x≤
2 ≤ x ≤1 3
0.08
u-approx. u-exact u e act
0.06 u(x)
⎧⎪ 4 ⎪⎪ x, ⎪⎪ 27 ⎪⎪ 4 1⎛ 1⎞ u( x ) = ⎪⎨ + ⎜⎜ x − ⎟⎟⎟, 3⎠ ⎪⎪ 81 27 ⎝ ⎪⎪ 5 5⎛ 2⎞ ⎪⎪⎪ − ⎜⎜ x − ⎟⎟, 3⎠ ⎪⎩ 81 27 ⎝
0.04
0.02
0
• Exact solution u( x ) =
0
0.2
0.4
x
0.6
0.8
1
1 x (1 − x 2 ) 6
– Three element solutions are poor – Need more elements
34
CONVERGENCE • Weighted residual of 2m 2m-th th order DE has highest derivatives of order m • With exact arithmetic, the following is sufficient for convergence to true solution (φ) as mesh is refined: – Complete polynomials of at least order m inside element – Continuity C ti it across element l tb boundaries d i up tto d derivatives i ti of order m-1 – Element must be capable of representing exactly uniform φ and uniform derivatives up to order m-1. • Beam: 4-th order DE (m = 2) – – – –
Complete polynomials: v(x) = a0 + a1x + a2x2 + a3x3 Continuity on v(x) and dv(x)/dx across element boundaries Uniform v(x) ( ) = a0 Beam elements will converge Uniform derivative dv(x)/dx = a1 upon refinement 35
RIGID BODY MOTION • Rigid Ri id b body d motion ti ffor CST can llead d tto non zero strains! t i !
• Rigid body motion
⎧u⎫ ⎡a1 ⎨ ⎬=⎢ ⎩v⎭ ⎣a4
cosθ −1 sinθ
⎧1 ⎫ −sinθ ⎤⎪ ⎪ ⎨x⎬ cosθ −1⎦⎥⎪ ⎪ ⎩y⎭
∂u
• The normal strain ε x = ∂x = cosθ −1 ≠ 0 36
CONVERGENCE RATE Quadratic curve u=a+bx+cx2 modeled by linear FE ufe=a+(b+ch)x a (b ch)x
• Maximal error at mid-point D
uA + uC ch2 h2 '' eD = uD − uB = − uB = = u 2 4 8
• Maximal gradient error is maximal at ends
eA' =
uC − u A h − b = hc = u '' 2 h
• Error in function converges faster than in derivative! 37
QUADRATIC ELEMENT FOR CUBIC SOLUTION
• Exact solution
u = a + bx + cx 2 + dx 3 • Finite element approximation
1 3 ⎞ ⎛ ⎞ ⎛ u = a + ⎜ b − dh 2 ⎟ x + ⎜ c + dh ⎟ x 2 2 2 ⎠ ⎝ ⎠ ⎝ • Maximal errors 3dh3 h3 eD = − =− u' ' ' 64 128
and
dh2 h2 e =− = − u' ' ' 2 12 ' A
38
CONVERGENCE RATE • Useful to know convergence rate – Estimate how much to refine – Detect modeling crimes – Extrapolate
• Most studies just do series of refinements if anything
39