User:Eml5526.s11.team7/HW4

Problem Description
Re do Problem 2.9, fixing $$\mathbf{\phi=\pi/2}$$ $$\Rightarrow \mathbf{F} = \left \{ 1, sinx, sin2x, sin3x, \cdot \cdot\cdot\cdot \right \}$$

For n=2
PART 1 From G1D1M.1.0/D2 $$P(u)=0 = \frac{\partial }{\partial x}\left[ 2\frac{\partial u}{\partial x}\right] + 3$$ therefore the exact solution is $$u(x)=\frac{-3}{4}x^2-4x+\frac{19}{4} $$ $$ d=\left \{ d_0, d_1, d_2 \right \}$$ $$ b=\left \{ 1,sinx,sin2x \right \}$$ $$ u^h(x) = \sum_{j=0}^2{d_jb_j(x)}=d_0 + d_{2}sin(2x) + d_{1}sin(x)$$ $$ (u^h)' =\frac{d{u}^{h}}{dx}(x) = \sum_{j=0}^2{d_jb_j'(x)}=d_{1}cos(x) + 2d_{2}cos(2x)$$ PART 2 The boundary conditions are enforced to generate two equations as follows: $$\displaystyle u^h(x=1) =0=d_{0} + 0.8414d_{1} + 0.9092d_{2}$$ $$\displaystyle (u^h)'(x=0) =-4=d_{1} + 2d_{2}$$  PART 3  An additional equation is generated by projected the residue onto a basis function as follows: $$\int_{0}^{1}b_{k}(x)P(u^h)dx=-0.9193d_1 - 5.664d_2=-3$$ PART 4  The equations are assembled into stiffness coefficient and force matrices in the form $$\mathbf{Kd}=\mathbf{f}$$:

$$ \begin{bmatrix} 1.0000 &   0.8415 &   0.9093\\         0   & 1.0000  &  2.0000\\         0  & -0.9194  & -5.6646   \end{bmatrix} \begin{bmatrix} d_0\\ d_1\\ d_2 \end{bmatrix} = \begin{bmatrix} 0\\ -4\\ -3 \end{bmatrix} $$ Note that the stiffness matrix, $$\mathbf{K}$$, is not symmetric PART 5  The coefficients are generated by solving the matrices above yielding:

$$ d=\begin{pmatrix} 4.7162\\  -7.4908\\    1.7454 \end{pmatrix} $$

For n=4
$$ d=\begin{pmatrix} 4.7491\\ -13.1144\\    8.2087\\   -3.6116\\    0.8829 \end{pmatrix} $$

For n=6
$$ d=\begin{pmatrix} 4.7500\\ -257.2135\\ 341.5393\\ -261.6484\\  128.7013\\  -38.5824\\    5.5312 \end{pmatrix} $$

Author
Eml5526.s11.team7.Srilalithkumar 05:05, 2 March 2011 (UTC)

Problem Description
“A First Course in Finite Elements” by Fish & Belytschko P73 Problem 3.5

Obtain the weak form of the equations of heat conduction with the following boundary conditions: The essential boundary condition defined at $$x=0$$ $$T(0) = 100$$ and the natural boundary condition defined at $$x=10$$ $$\bar q(10) = hT$$. The natural boundary condition is a convection condition.

Solution
The strong form of the heat problem is:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\frac{d}(AK\frac) + s = 0

$$     (4.2.1)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

T(0) = 100

$$     (4.2.2)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\bar q(10) = hT

$$     (4.2.3) So the weak form is obtained as follows:
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\int_0^{10} {w(\frac{d}(AK\frac) + s)} dx = 0

$$     (4.2.4)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\Rightarrow wAK\frac|_0^{10} - \int_0^{10} {\fracAK\fracdx + \int_0^{10} {wsdx} } = 0

$$     (4.2.5)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\Rightarrow w(10)AK\frac - w(0)AK\frac - \int_0^{10} {\fracAK\fracdx + \int_0^{10} {wsdx} } = 0

$$     (4.2.6)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\Rightarrow w(10)AhT - 0 - \int_0^{10} {\fracAK\fracdx + \int_0^{10} {wsdx} } = 0\;\;at\;\;w(0) = 0

$$     (4.2.7)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\Rightarrow w(10)AhT + \int_0^{10} {wsdx} = \int_0^{10} {\fracAK\fracdx}\;\;\;\; \forall w\;|\;w(0) = 0

$$     (4.2.8 )
 * 
 * }

To sum up, the weak form of the heat problem is:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

w(10)AhT + \int_0^{10} {wsdx} = \int_0^{10} {\fracAK\fracdx}\;\;\;\; \forall w\;|\;w(0) = 0

$$     (4.2.9)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

T(0) = 100

$$     (4.2.10)
 * <p style="text-align:right">
 * }

Author
Zongyi Yang

Brandon Hua

Problem Description

 * [From Fish & Belytschko, Chapter 3, Problem 3.6]

Given the strong form for the heat conduction problem in a circular plate:


 * $$k \frac{d}{dr}\left( r \frac{dT}{dr}\right)+rs=0 \qquad \text{on} \qquad 0 < r \leq R.$$


 * natural boundary condition: $$\displaystyle \frac{dT}{dr}(r=0) = 0,$$


 * essential boundary condition: $$\displaystyle T(r=R) = 0,$$


 * where R is the total radius of the plate, s is the heat source per unit length along the plate radius, T is the temperature and k is the conductivity. Assume that k, s, and R are given.

Find

 * 1) Construct the weak form for the above strong form.
 * 2) Use quadratic trial (candidate) solutions of the form $$\displaystyle T = \alpha_0 + \alpha_1 r + \alpha_2 r^2$$ and weight functions of the same form to obtain a solution of the weak form.
 * 3) Solve the differential equation with the boundary conditions and show that the temperature distribution along the radius is given by
 * $$ \displaystyle T = \frac{s}{4k} (R^2 - r^2)$$.

Part 1
First, the governing equation and natural boundary condition are multiplied by the weight function, $$\displaystyle w(r) $$, and integrated over the domain ($$ 0 < r \leq R$$). The restriction on the weight function is that $$\displaystyle w(R)=0$$.


 * {| style="width:100%" border="0"

$$ \int_{0}^{R}w(r) \left [ k \frac{d}{dr} \left ( r \frac{dT}{dr} \right ) + rs \right ] dr = 0 \qquad \text{and} \qquad w(r) \frac{dT}{dr} = 0.$$ (4.3.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The first equation can be expanded as


 * {| style="width:100%" border="0"

$$ \int_{0}^{R}w(r) k \frac{d}{dr} \left ( r \frac{dT}{dr} \right ) dr + \int_{0}^{R}w(r) r \ s \ dr = 0 \qquad \forall w.$$ (4.3.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The following integration by parts is simplified if the substitution $$\displaystyle f = r \frac{dT}{dr}$$ is made. Also note that the weight function will subsequently be referred to as $$\displaystyle w$$ with the dependence on the spatial domain, $$\displaystyle r $$, being implied. Equation 4.3.2 can therefore be written as,


 * {| style="width:100%" border="0"

$$ \int_{0}^{R}w k \frac{df}{dr} dr + \int_{0}^{R}wrs \ dr = 0.$$ (4.3.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Integration by parts yields,


 * {| style="width:100%" border="0"

$$ \int_{0}^{R} \frac{d}{dr} (wkf)dr - \int_{0}^{R} fk \frac{dw}{dr} dr + \int_{0}^{R}wrs \ dr = 0$$ (4.3.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ wkf|_{0}^{R} - \int_{0}^{R} kf \frac{dw}{dr} dr + \int_{0}^{R}wrs \ dr = 0.$$ (4.3.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The substitution in equation 4.3.3, $$\displaystyle f = r \frac{dT}{dr}$$, can now be reintroduced, and equation 4.3.5 becomes


 * {| style="width:100%" border="0"

$$ \left( wkr \frac{dT}{dr} \right)_{r=R} - \left( wkr \frac{dT}{dr} \right)_{r=0} - \int_{0}^{R} kr \frac {dT}{dr} \frac{dw}{dr} dr + \int_{0}^{R}wrs \ dr = 0.$$ (4.3.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The first term reduces to zero due to the prescribed weight function at $$\displaystyle r=R $$. In other words, $$\displaystyle w(R)=0 $$. The second term also reduces to zero since the natural boundary condition defines the heat flux at $$\displaystyle r=0 $$ to equal zero. Therefore, the weak form of the heat conduction problem becomes


 * {| style="width:100%" border="0"

$$ \int_{0}^{R} \frac{dw}{dr} kr \frac {dT}{dr} dr = \int_{0}^{R}wrs \ dr \qquad \forall w \quad \text{with} \quad w(R)=0.$$ (4.3.7)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Part 2
Consider the candidate and weight solutions of the following forms:


 * Candidate: $$\displaystyle T(r) = \alpha_0 + \alpha_1 r + \alpha_2 r^2$$,
 * Weight: $$\displaystyle w(r) = \beta_0 +\beta_1 r + \beta_2 r^2 $$.

For the candidate solution to be admissible, it must satisfy the essential boundary condition T(R) = 0. Likewise, the weight function must vanish at r = R. Implementing these conditions yields expressions with respect to the current and total radii of the plate.


 * {| style="width:100%" border="0"

$$ \displaystyle T(R) = \alpha_0 + \alpha_1 R + \alpha_2 R^2 = 0 \qquad \Rightarrow \qquad \alpha_0 = - \alpha_1 R - \alpha_2 R^2 $$     (4.3.8)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle w(R) = \beta_0 +\beta_1 R + \beta_2 R^2 = 0 \qquad \Rightarrow \qquad \beta_0 = - \beta_1 R - \beta_2 R^2 $$     (4.3.9)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Substituting these expressions for $$\displaystyle \alpha_0$$ and $$\displaystyle \beta_0$$ back into the candidate solution and weight function,


 * {| style="width:100%" border="0"

$$ \displaystyle T(r) = \alpha_1 \left( r-R \right) + \alpha_2 \left( r^2 - R^2 \right) $$     (4.3.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle w(r) = \beta_1 \left( r-R \right) + \beta_2 \left( r^2 - R^2 \right) $$     (4.3.11)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The derivatives of equations 4.3.10 and 4.3.11


 * {| style="width:100%" border="0"

$$ \frac{dT}{dr} = \alpha_1 + 2 \alpha_2 r,$$ (4.3.12)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \frac{dw}{dr} = \beta_1 + 2 \beta_2 r. $$ (4.3.13)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Equations 4.3.10 through 4.3.13 are substituted into the weak form, equation 4.3.7:


 * {| style="width:100%" border="0"

$$ \int_{0}^{R} \left( \beta_1 + 2 \beta_2 r \right) kr \left( \alpha_1 + 2 \alpha_2 r \right) dr = \int_{0}^{R} \left[ \beta_1 (r-R) + \beta_2 (r^2 - R^2) \right] rs \ dr $$ (4.3.14)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Integrating the both sides of the equation over the domain yields


 * {| style="width:100%" border="0"

$$ \text{LHS:} \quad \int_{0}^{R} \left( \beta_1 + 2 \beta_2 r \right) kr \left( \alpha_1 + 2 \alpha_2 r \right) dr \quad = \quad \frac{2}{3} \alpha_1 \beta_2 k R^3 \ + \ \frac{1}{2} \alpha_1 \beta_1 k R^2 \ + \ \frac{2}{3} \alpha_2 \beta_1 k R^3 \ + \ \alpha_2 \beta_2 k R^4 $$     (4.3.15)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \text{RHS:} \quad \int_{0}^{R} \left[ \beta_1 (r-R) + \beta_2 (r^2 - R^2) \right] rs \ dr \quad = \quad - \frac{1}{4} \beta_2 s R^4 \ - \ \frac{1}{6} \beta_1 s R^3 $$     (4.3.16)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Equations 4.3.15 and 4.3.16 are substituted into 4.3.14, and the entire equation is rewritten to collect the terms $$\displaystyle \beta_1$$ and $$\displaystyle \beta_2$$.


 * {| style="width:100%" border="0"

$$ \beta_1 \left( \frac{1}{2} \alpha_1 k R^2 + \frac{2}{3} \alpha_2 k R^3 + \frac{1}{6} s R^3 \right ) + \beta_2 \left( \frac{2}{3} \alpha_1 k R^3 + \alpha_2 k R^4 + \frac{1}{4} s R^4 \right) =0 $$     (4.3.17)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The coefficients $$\displaystyle \alpha_1$$ and $$\displaystyle \alpha_2$$ can be determined in terms of the known components s, k, and R.


 * {| style="width:100%" border="0"

$$ \displaystyle \begin{bmatrix} {\frac{1}{2}R} & { \frac{2}{3}R^3} \\ \\ {\frac{2}{3}R^3} & {R^4} \end{bmatrix} \cdot \begin{bmatrix} \alpha_1\\ \\ \alpha_2 \end{bmatrix} = \begin{bmatrix} {- \frac{1}{6}sR^3} \\ \\ {- \frac{1}{4}sR^4} \end{bmatrix}
 * style="width:95%" |
 * style="width:95%" |

$$     (4.3.18)
 * <p style="text-align:right">
 * }


 * Therefore,


 * {| style="width:100%" border="0"

$$ \alpha_1 = - \frac{ks}{12} \left( 2R^7 + R^5 \right) \qquad \text{and} \qquad \alpha_2 = - \frac{ks}{36} \left( 9R^8 + 4R^6 \right) $$     (4.3.19)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Substituting the values of $$\displaystyle \alpha_1$$ and $$\displaystyle \alpha_2$$ into equation 4.3.10 yields the candidate solution in terms of known values. There is some concern regarding the validity of this equation, since the system units do not correspond to an appropriate unit measurement of temperature.


 * {| style="width:100%" border="0"

$$ T(r) = - \frac{ks}{12} \left( 2R^7 + R^5 \right)(r-R) - \frac{ks}{36} \left( 9R^8 + 4R^6 \right)(r^2 - R^2) $$     (4.3.20)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

The candidate (trial) solution above is substituted into the governing equation to confirm equilibrium.


 * {| style="width:100%" border="0"

$$ k \frac{d}{dr}\left( r \frac{dT}{dr}\right)+rs=0 \qquad \Rightarrow \qquad -k \left[ \frac{ks(2R^7 + R^5)}{12} + \frac{rks(9R^8 + 4R^6)}{9} \right] +rs = 0 $$     (4.3.21)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * Solving for $$\displaystyle r$$ and simplifying, we are given an expression for when $$\displaystyle r$$ is valid for given values of k and R. Note that the validity of $$\displaystyle r$$ is independent of the heat source per unit length, s.


 * {| style="width:100%" border="0"

$$ r = - \frac{6R^7 k^2 + 3R^5k^2}{36R^8k^2 + 16R^6k^2 - 36} $$     (4.3.22)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Part 3
The differential - or governing - equation below will now be solved for $$\displaystyle T(r)$$.


 * {| style="width:100%" border="0"

$$ k \frac{d}{dr}\left( r \frac{dT}{dr}\right)+rs=0 \qquad \text{on} \qquad 0 < r \leq R $$ (4.3.23)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

The equation is rewritten with $$\displaystyle rs$$ and $$\displaystyle k$$ moved to the right hand side. Then both sides can be integrated over $$\displaystyle r$$.


 * {| style="width:100%" border="0"

$$ \int \frac{d}{dr} \left( r \frac{dT}{dr} \right) dr = - \int \frac{s}{k} r dr $$ (4.3.24)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

This reduces to


 * {| style="width:100%" border="0"

$$ r \frac{dT}{dr} = - \frac{s}{2k} r^2 + constant $$     (4.3.25)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Recalling the natural boundary condition and applying it to equation 4.3.25,


 * {| style="width:100%" border="0"

$$ \frac{dT}{dr}(0) = 0 \quad \Rightarrow \quad r \frac{dT}{dr} = - \frac{s}{2k} r^2 + constant $$     (4.3.26)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ 0 = 0 + constant \qquad \text{therefore} \qquad r \frac{dT}{dr}  = - \frac{s}{2k} r^2. $$     (4.3.27)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Rewrite as


 * {| style="width:100%" border="0"

$$ dT = - \frac{s}{2k} r dr \qquad \Rightarrow \qquad \int dT  = - \int \frac{s}{2k} r dr. $$ (4.3.28)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Integrate recalling the essential boundary condition, $$\displaystyle T(R) = 0$$.


 * {| style="width:100%" border="0"

$$ T = - \frac{s}{4k} r^2 + constant \quad \Rightarrow \quad T(R) = - \frac{s}{4k} R^2 + constant = 0 $$     (4.3.29)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * $$\quad \Rightarrow constant = \frac{s}{4k} R^2$$

This confirms that the temperature distribution along the radius is given by


 * {| style="width:100%" border="0"

$$ T(r) = \frac{s}{4k} \left( R^2 - r^2 \right). $$     (4.3.30)
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

Authors
Philip Flater Braden Snook

Problem Description
HW 4.4 p.21.1 Problem 3.7 from "A First Course in Finite Elements" by Fish and Belytschko

Given the strong form for the circular bar in torsion:
 * {| style="width:100%" border="0"

$$ \displaystyle \frac{d}{dx}\left ( JG\frac{d\phi }{dx} \right )+m=0, 0\leq x\leq l $$ (4.4.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Essential Boundary Condition:
 * {| style="width:100%" border="0"

$$ \displaystyle \phi \left ( x=0 \right )=\overline{\phi } $$       (4.4.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Natural Boundary Condition:
 * {| style="width:100%" border="0"

$$ \displaystyle M\left ( x=l \right )=\left ( JG\frac{d\phi }{dx} \right )_{l}=\overline{M} $$       (4.4.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Where:

$$ \displaystyle m(x) $$ is a distributed moment per unit length $$ \displaystyle M $$ is the torsion moment $$ \displaystyle \phi $$ is the angle of rotation $$ \displaystyle G $$ is the shear modulus $$ \displaystyle C $$ is the radius of the circular shaft $$ \displaystyle J $$ is the polar moment of inertia given by the equation


 * {| style="width:100%" border="0"

$$ \displaystyle J=\frac{\pi C^{4}}{2}$$ (4.4.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

a) Contruct the weak form for the circular bar in torsion. b) Assume that (m=0) and integrate the differential equation given above (4.4.1). Find the integration constants using the boundary conditions (4.4.2) and (4.4.3).

Hint: "I suggest to use this opportunity to solve the G1DM1.0 / D2 (at least a slightly modified version of it) as stated on p.9-2 of Mtg 9 (e); you just need to swap the boundary locations around to conform to the statement for HW 4.4, i.e., let x=0 be the location for the ess. b.c., and x=1 the location for the nat. b.c. So use the same data as stated in G1DM1.0 / D2, with the boundary locations swapped around; we can call this problem G1DM1.0 / D2b (data set 2b). This way, you can plot the exact solution, which can then be used to verify the results for future problems." -Prof Vu Quoc

Part a
Mulitply (4.4.1) by an arbitrary function $$ \displaystyle w(x) $$ and integrate over the domain $$ \displaystyle \Omega = ]0,l[ $$


 * {| style="width:100%" border="0"

$$ \displaystyle \int_{0}^{l}w\left [ \frac{d}{dx}JG\frac{d\phi}{dx} \right ]dx+\int_{0}^{l}wMdx=0, \forall w $$ (4.4.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle \frac{d}{dx}\left ( wf \right )=\frac{d(wf)}{dx}-f\frac{dw}{dx} $$       (4.4.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle \int_{0}^{l}w\frac{dw}{dx}=\int_{0}^{l}\frac{d}{dx}(wf)dx-\int_{0}^{l}f\frac{dw}{dx}dx $$       (4.4.7)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle \int_{0}^{l}w\frac{dw}{dx}=(wf)_{x=l}-(wf)_{x=0}-\int_{0}^{l}f\frac{dw}{dx}dx $$       (4.4.8)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Let
 * {| style="width:100%" border="0"

$$ \displaystyle f=JG\frac{d\phi}{dx} $$       (4.4.9)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

(4.4.5) becomes
 * {| style="width:100%" border="0"

$$ \displaystyle \left (wJG\frac{d\phi}{dx} \right )_{x=l}-\left (wJG\frac{d\phi}{dx}  \right )_{x=0}-\int_{0}^{l}\frac{dw}{dx}JG\frac{d\phi}{dx}dx-\int_{0}^{l}wmdx=0,\forall w, w(l)=0 $$       (4.4.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle G\frac{d\phi}{dx}=\tau $$       (4.4.11)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

(4.4.10) becomes
 * {| style="width:100%" border="0"


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * {| style="width:100%" border="0"
 * {| style="width:100%" border="0"

$$ \displaystyle \left (wJ \tau  \right )_{x=l}-\left (wJ \tau   \right )_{x=0}-\int_{0}^{l}\frac{dw}{dx}JG\frac{d\phi}{dx}dx-\int_{0}^{l}wmdx=0,\forall w $$ (4.4.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * style = |
 * }
 * }

Recall:


 * {| style="width:100%" border="0"

$$ \displaystyle w(0)=0 $$       (4.4.12)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Hint: "The problem stated on p.9-2 of Mtg 9 (e) corresponds to the G1DM1.0 / D1 (data set 1, in which the coefficient a_2 is a linear function of x), and NOT the G1DM1.0 / D2 (data set 2, in which the coefficient a_2 is a constant). To be precise, the problem that I referred to on p.9-2 of Mtg 9 (e) is stated in Eq.(3) on that page.  Hence, let's define G1DM1.0 / D1b (data set 1b, NOT 2b) as the problem when you swap the boundary conditions around, i.e., let x=0 be the location for the ess. b.c., and x=1 the location for the nat. b.c." -Prof Vu-Quoc From (3) 9-2


 * {| style="width:100%" border="0"

$$ \displaystyle \frac{d}{dx}\left [ \left ( 2+3x \right ) \frac{du}{dx} \right ]+5x=0,\forall x\varepsilon ]0,1[ $$ $$ \displaystyle \Rightarrow JG=2+3x $$ $$ \displaystyle \Rightarrow m(x)=5x $$       (4.4.13)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle u(0)=4 $$ $$ \displaystyle \Rightarrow \phi(0)=4 $$       (4.4.14)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle \frac{du}{dx}(x=1)=-6  $$ $$ \displaystyle \Rightarrow \frac{d\phi}{dx}(x=1)=-6 $$       (4.4.15)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Let
 * {| style="width:100%" border="0"

$$ \displaystyle w(0)=0 $$ (4.4.10) becomes:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \left (w(2+3x)\frac{d\phi}{dx} \right )_{x=l}-\int_{0}^{1}\frac{dw}{dx}(2+3x)\frac{d\phi}{dx}dx-\int_{0}^{1}w(5x)dx=0,\forall w, w(0)=0 $$     (4.4.16)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle w(2+3*1)(-6) -\int_{0}^{1}\frac{dw}{dx}(2+3x)\frac{d\phi}{dx}dx-\int_{0}^{1}w(5x)dx=0 $$     (4.4.17)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle 30w -\int_{0}^{1}\frac{dw}{dx}(2+3x)\frac{d\phi}{dx}dx-\int_{0}^{1}w(5x)dx=0 $$     (4.4.18)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle \int_{0}^{1}\frac{dw}{dx}(2+3x)\frac{d\phi}{dx}dx=30w -\int_{0}^{1}w(5x)dx $$     (4.4.19)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * {| style="width:100%" border="0"
 * {| style="width:100%" border="0"

$$ \displaystyle \int_{0}^{1}\frac{dw}{dx}(2+3x)\frac{d\phi}{dx}dx=30w(x) -\int_{0}^{1}w(x)(5x)dx $$       (4.4.20)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * style = |
 * }
 * }

Part b
From (4.4.23)
 * {| style="width:100%" border="0"

$$ \displaystyle \frac{d}{dx}\left ( (2+3x)\frac{d\phi}{dx} \right )+5x=0 $$     (4.4.21)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle \int d(2+3x)\frac{d\phi}{dx}+\int5xdx=0 $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle (2+3x)\frac{d\phi}{dx}+ \frac{5x^{2}}{2}+C_{1}=0 $$     (4.4.22)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \int (2+3x)d\phi= \frac{- \frac{5x^{2}}{2}+C_{1}}{(2+3x)} $$     (4.4.23)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \phi= \int - \frac{5x^{2}}{2(2+3x)}+\int \frac{C_{1}}{(2+3x)} $$     (4.4.24)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

From 2.4.4
 * {| style="width:100%" border="0"

$$ \displaystyle \int \frac{x^2}{2+3x}dx=\frac{9x^{2}-12x+8log(2+3x)}{54} $$     (4.4.25)  (4.4.24) becomes:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * {| style="width:100%" border="0"
 * {| style="width:100%" border="0"

$$ \displaystyle \phi= -\frac{5}{2}\left (\frac{9x^{2}-12x+8log(2+3x)}{54} \right ) + \frac{C_{1}}{3} {log(2+3x)}+C_{2} $$       (4.4.26) Apply (4.4.15) to (4.4.22)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * style = |
 * }
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle (2+3*1)*-6+ \frac{5*1^{2}}{2}+C_{1}=0$$ $$ \displaystyle C_{1}=-27.5 $$     (4.4.27)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Apply (4.4.27) to (4.4.26)
 * {| style="width:100%" border="0"

$$ \displaystyle \phi= -\frac{5}{2}\left (\frac{9x^{2}-12x+8log(2+3x)}{54} \right ) + \frac{-27.5}{3} {log(2+3x)}+C_{2} $$     (4.4.28)  Apply (4.4.14) to (4.4.28)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \phi(x=0)=4= -\frac{5}{2}\left (\frac{9*0^{2}-12*0+8log(2+3*0)}{54} \right ) + \frac{-27.5}{3} {log(2+3*0)}+C_{2} $$ $$ \displaystyle C_{2}=6.8709342 $$     (4.4.29)  The final solution is:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"


 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * {| style="width:100%" border="0"
 * {| style="width:100%" border="0"

$$ \displaystyle \phi= -\frac{5}{2}\left (\frac{9x^{2}-12x+8log(2+3x)}{54} \right ) + \frac{-27.5}{3} {log(2+3x)}+ 6.8709342 $$       (4.4.26)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * style = |
 * }
 * }

Author
Brandonhua 04:06, 2 March 2011 (UTC)

Problem Description
Possible basis functions are given as:

Polynomial: $${{F}_{P}}=\left\{ {{x}^{j}},\text{ }j=0,1,2,\ldots \right\}$$

Cosine: $${{F}_{C}}=\left\{ \cos jx,\text{ }j=0,1,2,\ldots \right\}$$

Sine: $${{F}_{S}}=\left\{ 1,\text{ }\sin jx,\text{ }j=1,2,\ldots \right\}$$

Fourier: $${{F}_{F}}=\left\{ \cos jx,\text{ }j=0,1,2,\ldots ;\text{ }\sin kx,\text{ }k=1,2,\ldots \right\}$$

Exponential: $${{F}_{E}}=\left\{ {{e}^{jx}},\text{ }j=0,1,2,\ldots \right\}$$

For $${{\Gamma }_{g}}=\left\{ \beta \right\}$$, and each of the above families of basis functions, i.e., $${{\Gamma }_{I}}=\left\{ {{b}_{j}}(x),\text{ }j=0,1,\ldots \right\},\text{ }I\in \left\{ P,C,S,F,E \right\}$$

Find
Find corresponding family $${{\bar{F}}_{I}}=\left\{ {{{\bar{b}}}_{j}} \right\}$$ satisfying constraint breaking solution (CBS), i.e., $${{b}_{0}}(\beta )\ne 0$$ and $${{b}_{i}}(\beta )=0,\text{ }i=1,2,\ldots ,n$$. Let $$\Omega =\left] \alpha ,\beta \right[=\left] -2,4 \right[$$.

1) $${{\bar{b}}_{0}}(x)=\text{const}$$, plot $${{b}_{j}}$$ and $${{\bar{b}}_{j}}$$ for $$j=1,2,3$$.

2) Show $$\left\{ {{e}^{jx}},\text{ }j=0,1,2 \right\}$$ linear independent on $$\Omega $$.

Part 1
The corresponding family $${{\bar{F}}_{I}}=\left\{ {{{\bar{b}}}_{j}} \right\}$$ satisfying CBS:

Polynomial:

$${{F}_{P}}=\left\{ {{\left( x-4 \right)}^{j}},\text{ }j=0,1,2,\ldots \right\}$$

Cosine:

$${{F}_{C}}=\left\{ 1,\text{ }\cos \left( j\left( x-4 \right)+\frac{\pi }{2} \right),\text{ }j=1,2,\ldots \right\}$$

Sine:

$${{F}_{S}}=\left\{ 1,\text{ }\sin \left( j\left( x-4 \right) \right),\text{ }j=1,2,\ldots \right\}$$

Fourier:

$${{F}_{F}}=\left\{ 1,\text{ }\cos \left( j\left( x-4 \right)+\frac{\pi }{2} \right),\text{ }j=1,2,\ldots ;\text{ }\sin \left( k\left( x-4 \right) \right),\text{ }k=1,2,\ldots \right\}$$

Exponential:

$${{F}_{E}}=\left\{ 1,\text{ }{{e}^{j\left( x-4 \right)}}-1,\text{ }j=1,2,\ldots \right\}$$

The following figures show the polynomial basis functions (j=1,2,3):





The following figures show the cosine basis functions (j=1,2,3):





The following figures show the sine basis functions (j=1,2,3):





The following figures show the fourier basis functions (j=1,2,3):





The following figures show the exponential basis functions (j=1,2,3):





Part 2
The exponential basis function is $${{b}_{j}}=\left\{ {{e}^{jx}},\text{ }j=0,1,2 \right\}$$,

Thus the Gram matrix is:

$$\Gamma =\left[ \begin{matrix} 6 & {{e}^{4}}-{{e}^{-2}} & \frac{1}{2}\left( {{e}^{8}}-{{e}^{-4}} \right) \\ {{e}^{4}}-{{e}^{-2}} & \frac{1}{2}\left( {{e}^{8}}-{{e}^{-4}} \right) & \frac{1}{3}\left( {{e}^{12}}-{{e}^{-6}} \right) \\ \frac{1}{2}\left( {{e}^{8}}-{{e}^{-4}} \right) & \frac{1}{3}\left( {{e}^{12}}-{{e}^{-6}} \right) & \frac{1}{4}\left( {{e}^{16}}-{{e}^{-8}} \right) \\ \end{matrix} \right]$$

And the determinant of Gram matrix is:

$$\det \left( \Gamma \right)=1.1145\times {{10}^{9}}\ne 0$$

Thus, the exponential basis function is linear independent on $$\Omega $$.

Author
Jiang Jin

Problem Description
An elastic bar is subjected along its length to a variable distributed spring p(x) as seen in the figure below. An axial force proportional to the displacement is imposed by the spring on the bar. Adapted from "A First Course in Finite Elements" by Fish and Belytschko p. 74 pb.3.9

a. Construct the strong form.

b. Construct the weak form.



Given
This bar is of length l, cross-sectional area A(x), Young's modulus E(x) with body force b(x) and with boundary conditions:


 * {| style="width:100%" border="0"

$$\displaystyle u(0)=0$$
 * style="width:10%; |
 * style="width:10%; |

$$\displaystyle E(l)A(l)\frac{du(l)}{dx}=t$$
 * <p style="text-align:right">
 * }

As seen in the figure above.

Solution
Assumptions:
 * The bar must be in equilibrium
 * It satisfies Hooke's law
 * Compatible Displacement Field
 * The Strain-Displacement Equation must be Satisfied

a. Construct the Strong Form
Following a similar solution method as Team 2 used here with Eq. 1.1-1.10. Using the following diagram from Team 7 homework 1 to start creating the strong form:



Performing a force balance on the rod using the free body diagram above, the net force is equal to the internal forces plus the external, or applied, forces. Note, in this specific problem we will designate our body force as b(x+dx/2) (instead of f(x+dx/2 as seen in the diagram from Team 7) and have the additional body force in the opposite direction caused by the springs (as seen by the figure in the problem description) designated as p(x+dx/2)


 * {| style="width:100%" border="0"

$$ \displaystyle \Sigma \text{ Internal Force } + \text{ Applied Force } = \text{ Net Force}$$ (4.6.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \Sigma \text{ Internal Force }=N(x) + N(x+\Delta x)$$ (4.6.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \text{ Applied Force }=b(x+\frac{dx}{2})dx-p(x+\frac{dx}{2})dx$$ (4.6.3)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Since we are looking at the static case:


 * {| style="width:100%" border="0"

$$ \displaystyle \text{ Net Force}= m*a =0$$ (4.6.4)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Substituting equations 4.6.2-4.6.4 into the force balance equation, 4.6.1:


 * {| style="width:100%" border="0"

$$ \displaystyle N(x) + N(x+dx) + b(x+\frac{dx}{2})dx-p(x+\frac{dx}{2})dx =0$$ (4.6.5)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Plugging in for the direction of the surface normal, the balance of forces at position x can be expressed as:


 * {| style="width:100%" border="0"

$$ \displaystyle N(x)=n(x)\sigma (x)A(x)=-\sigma (x)A(x)$$ (4.6.6)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Also plugging in for the surface normal, the corresponding force at x + dx:


 * {| style="width:100%" border="0"

$$ \displaystyle N(x+dx)=n(x+dx)\sigma (x+dx)A(x+dx)=\sigma (x+dx)A(x+dx)$$ (4.6.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Using a Taylor Series:


 * {| style="width:100%" border="0"

$$ \displaystyle \sigma (x+dx)=\sigma (x)+\frac{d\sigma (x)}{dx}dx+...$$ (4.6.8)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle A(x+dx)=A(x)+\frac{dA(x)}{dx}dx+...$$ (4.6.9)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Additionally, in Eq. 4.6.3 the body force $$b(x+\frac{dx}{2})$$ has been averaged for the differential volume we are considering to formulate the governing differential equation. The same has been done for the force caused by the variable distributed spring, given as $$p(x+\frac{dx}{2})$$. Taylor expanding for these as well:


 * {| style="width:100%" border="0"

$$\displaystyle b(x+\frac{dx}{2})=b(x)+\frac{db(x)}{dx}\frac{dx}{2}....$$. (4.6.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle p(x+\frac{dx}{2})=p(x)+\frac{dp(x)}{dx}\frac{dx}{2}....$$. (4.6.11)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Plugging in Eq.'s 4.6.8 and 4.6.9 into Eq. 4.6.7, and then plugging this and Eq.'s 4.6.6, 4.6.10 and 4.6.11 into Eq. 4.6.5 gives:
 * {| style="width:100%" border="0"

$$\displaystyle \sigma (x)A(x) + \sigma (x)\frac{dA(x)}{dx}dx+A(x)\frac{d\sigma (x)}{dx}-\sigma (x)A(x)+b(x)dx-p(x)dx$$. (4.6.12)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Using the product rule, canceling out terms, and dividing through by dx gives:
 * {| style="width:100%" border="0"

$$\displaystyle \frac{d}{dx}\sigma (x)A(x)+b(x)-p(x)=0$$. (4.6.13)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Noting that the strain-displacement equation is
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\varepsilon \left( x \right) = \frac

$$     (4.6.14) And by Hooke ‘ s Law
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\sigma \left( x \right) = E\left( x \right)\varepsilon \left( x \right) = E\left( x \right)\frac

$$     (4.6.15)
 * <p style="text-align:right">
 * }

Gives the following strong form:
 * {| style="width:92%; padding:10px; border:2px solid #8888aa"

$$\displaystyle \frac{d}{dx}E(x)A(x)\frac{du}{dx}+b(x)-p(x)=0$$. (4.6.16)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

b. Construct the Weak Form
Following the same solutions as Team 7's previous work here.

In order to reduce the requirements or constraints on what trial solutions could be substituted in for u, we will generate the weak form which requires the basis solutions only be once differentiable. This occurs from the natural boundary condition being absorbed in the weak form, as will be seen shortly. Also, this creates a symmetric stiffness matrix.

As with Team 7's solution, the strong form and natural boundary condition are used in an inner product with an arbitrary weighting function giving the following equations:


 * {| style="width:100%" border="0"

$$\displaystyle \int_{0}^{l}w(x)\big [\frac{d}{dx}E(x)A(x)\frac{du}{dx}+b(x)-p(x)\big ]dx=0$$. (4.6.17)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle w(l)E(l)A(l)\frac{du(l)}{dx}-w(l)t=0$$. (4.6.18)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Distributing the weighting function w(x) and the integral of Eq. 4.6.17 and looking at the first term:


 * {| style="width:100%" border="0"

$$\displaystyle \int_{0}^{l}w\frac{d}{dx}E(x)A(x)\frac{du}{dx}dx=w(x)E(x)A(x)|_{0}^{l}-int_{0}^{l}\frac{dw}{dx}E(x)A(x)\frac{du}{dx}dx$$. (4.6.19)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

In Eq. 4.6.19 integration by parts (derived and example of use here) was used. Now to remove the unknown flux at the essential boundary location (x=0), we will make the arbitrary weighting function be zero at x equal to zero, or


 * {| style="width:100%" border="0"

$$\displaystyle w(0)=0$$. (4.6.20) Using Eq. 4.6.7 and 4.6.9, Eq. 4.6.8 becomes:
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle w(x)E(x)A(x)|_{0}^{l}-\int_{0}^{l}\frac{dw}{dx}E(x)A(x)\frac{du}{dx}dx=w(l)t-\int_{0}^{l}\frac{dw}{dx}E(x)A(x)\frac{du}{dx}dx$$. (4.6.21)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Combining Equations 4.6.18 and 4.6.21 gives the following weak form:


 * {| style="width:92%; padding:10px; border:2px solid #8888aa"

$$\displaystyle w(l)t-\int_{0}^{l}\frac{dw}{dx}E(x)A(x)\frac{du}{dx}dx+\int_{0}^{l}w(x)p(x)dx+\int_{0}^{l}w(x)b(x)dx$$. (4.6.22)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Note in equation 4.6.22 the weighting function must be zero at the essential boundary condition (w(0)=0). Also, the weighting function and the displacement are both first order derivatives and now a symmetric stiffness matrix can be formed.

Authors
Braden Snook, Zongyi Yang

Problem Description
A calculix manual for novice/dummies

Author
Eml5526.s11.team7.Srilalithkumar

Problem Description
Consider a trial solution of the form $$u(x) = {\alpha _0} + {\alpha _1}(x - 3) + {\alpha _2}{(x - 3)^2} $$and a weight function of the same form. Obtain a solution to the weak form in Problem 3.1, "A First Course in Finite Elements", by Fish & Belytschko P72. Find $$\tilde M$$, assuming $$A = 1,E = 2,\bar m = 3$$, Do for $$n = 3$$. For dynamics, with $${u^h}(\beta ,t) = g(t) = \sin 2t$$, find $$F(t)$$ The weak form of Problem 3.1 is
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\int\limits_1^3 {\fracAE\frac} dx = - 0.1{(wA)_{x = 1}} + \int\limits_1^3 {2xwdx\forall w(3) = 0}

$$     (4.8.1)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

u(3) = 0.001

$$
 * <p style="text-align:right">
 * }

Solution
Part 1: The basis functions $${b_i} = \{ 1,(x - 3),{(x - 3)^2},{(x - 3)^3}\} i = 0,1,2,3$$
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde M = \left[ {\begin{array}{ccccccccccccccc} & \\  & \end{array}} \right]

$$     (4.8.2)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{M_{EE}} = {\tilde M_{00}} = \tilde m({b_0},{b_0}) = \int\limits_1^3 {{b_0}\bar m{b_0}dx} = \int\limits_1^3 {3dx}  = 6

$$     (4.8.3)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{M_{EF}} = \left[ {{M_{0j}};j = 1,2,3} \right] = \tilde m({b_0},{b_j})

$$     (4.8.4)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_0},{b_1}) = \int\limits_1^3 {{b_0}\bar m{b_1}dx} = \int\limits_1^3 {3(x - 3)dx}  =  - 6

$$     (4.8.5)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_0},{b_2}) = \int\limits_1^3 {{b_0}\bar m{b_2}dx} = \int\limits_1^3 {3{{(x - 3)}^2}dx}  = 8

$$     (4.8.6)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_0},{b_3}) = \int\limits_1^3 {{b_0}\bar m{b_3}dx} = \int\limits_1^3 {3{{(x - 3)}^3}dx}  =  - 12

$$     (4.8.7)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{M_{EF}} = \left[ {\begin{array}{ccccccccccccccc} { - 6}&8&{ - 12} \end{array}} \right]

$$     (4.8.8)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{M_{FE}} = {M_{EF}}^T = \left[ {\begin{array}{ccccccccccccccc} { - 6} \\  8 \\   { - 12} \end{array}} \right]

$$     (4.8.9)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{M_{FF}} = \left[ {{M_{ij}};i,j = 1,2,3} \right] = \tilde m({b_i},{b_j})

$$     (4.8.10)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_1},{b_1}) = \int\limits_1^3 {{b_1}\bar m{b_1}dx} = 3\int\limits_1^3 {(x - 3)(x - 3)dx}  = 8

$$     (4.8.11)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_1},{b_2}) = \int\limits_1^3 {{b_1}\bar m{b_2}dx} = 3\int\limits_1^3 {(x - 3){{(x - 3)}^2}dx}  =  - 12

$$     (4.8.12)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_1},{b_3}) = \int\limits_1^3 {{b_1}\bar m{b_3}dx} = 3\int\limits_1^3 {(x - 3){{(x - 3)}^3}dx}  = \frac{5}

$$     (4.8.13)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_2},{b_1}) = \int\limits_1^3 {{b_2}\bar m{b_1}dx} = 3\int\limits_1^3 {{{(x - 3)}^2}(x - 3)dx}  =  - 12

$$     (4.8.14)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_2},{b_2}) = \int\limits_1^3 {{b_2}\bar m{b_2}dx} = 3\int\limits_1^3 {{{(x - 3)}^2}{{(x - 3)}^2}dx}  = \frac{5}

$$     (4.8.15)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_2},{b_3}) = \int\limits_1^3 {{b_2}\bar m{b_3}dx} = 3\int\limits_1^3 {{{(x - 3)}^2}{{(x - 3)}^3}dx}  =  - 32

$$     (4.8.16)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_3},{b_1}) = \int\limits_1^3 {{b_3}\bar m{b_1}dx} = 3\int\limits_1^3 {{{(x - 3)}^3}(x - 3)dx}  = \frac{5}

$$     (4.8.17)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_3},{b_2}) = \int\limits_1^3 {{b_3}\bar m{b_2}dx} = 3\int\limits_1^3 {{{(x - 3)}^3}{{(x - 3)}^2}dx}  =  - 32

$$     (4.8.18)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde m({b_3},{b_3}) = \int\limits_1^3 {{b_3}\bar m{b_3}dx} = 3\int\limits_1^3 {{{(x - 3)}^3}{{(x - 3)}^3}dx}  = \frac{7}

$$     (4.8.19)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\therefore {M_{FF}} = \left[ {\begin{array}{ccccccccccccccc} 8&{ - 12}&{\frac{5}} \\ { - 12}&{\frac{5}}&{ - 32} \\ {\frac{5}}&{ - 32}&{\frac{7}} \end{array}} \right]

$$     (4.8.20)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

\therefore {M_{FF}} = \left[ {\begin{array}{ccccccccccccccc} 6&{ - 6}&8&{ - 12} \\  { - 6}&8&{ - 12}&{\frac{5}} \\ 8&{ - 12}&{\frac{5}}&{ - 32} \\ { - 12}&{\frac{5}}&{ - 32}&{\frac{7}} \end{array}} \right]

$$     (4.8.21) Part 2:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

F = {F_F} - \left( {{M_{FE}}{g^{(2)}} + {K_{FE}}g} \right)

$$     (4.8.22)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{F_F} = \left\{ {{F_i},i = 1,...,n} \right\} = \left\{ {\tilde f\left( \right),i = 1,...,n} \right\}

$$     (4.8.23)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\because \tilde f({b_i}) = {b_i}(\alpha )h + \int\limits_\alpha ^\beta {{b_i}fdx}

$$     (4.8.24)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\because f = 0

$$     (4.8.25)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\therefore \tilde f({b_i}) = {b_i}(1)h

$$     (4.8.26)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

h = {a_2}(\alpha )\frac = AE\frac = 0.1

$$     (4.8.27)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde f({b_1}) = {b_1}(1)h = - 2 \times 0.1 =  - 0.2

$$     (4.8.28)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde f({b_2}) = {b_2}(1)h = {\left( { - 2} \right)^2} \times 0.1 = 0.4

$$     (4.8.29)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde f({b_3}) = {b_3}(1)h = {\left( { - 2} \right)^3} \times 0.1 = - 0.8

$$     (4.8.30)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\therefore {F_F} = \left[ {\begin{array}{ccccccccccccccc} { - 0.2} \\  {0.4} \\   { - 0.8} \end{array}} \right]

$$     (4.8.31)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{g^{(2)}} = - 4\sin (2t)

$$     (4.8.32)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{K_{FE}} = \left\{ {{K_{0j}};j = 1,...,n} \right\} = \left\{ {\tilde k\left( {{b_0},{b_j}} \right),j = 1,...,n} \right\}

$$     (4.8.33)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\because \tilde k({b_i},{b_j}) = \int\limits_1^3 {b{'_i}{a_2}b{'_j}dx} = \int\limits_1^3 {b{'_i}AEb{'_j}dx}

$$     (4.8.34) So
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\tilde k({b_0},{b_1}) = \int\limits_1^3 {b{'_0}{a_2}b{'_1}dx} = \int\limits_1^3 {b{'_0}AEb{'_1}dx}  = 0 = \tilde k({b_0},{b_2}) = \tilde k({b_0},{b_3}) = 0

$$     (4.8.35)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\therefore {K_{EF}} = \left[ {\begin{array}{ccccccccccccccc} 0 \\  0 \\   0 \end{array}} \right]

$$     (4.8.36) And：
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

\because {M_{FE}} = {M_{EF}}^T = \left[ {\begin{array}{ccccccccccccccc} { - 6} \\  8 \\   { - 12} \end{array}} \right]

$$     (4.8.37)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

F = {F_F} - \left( {{M_{FE}}{g^{(2)}} + {K_{FE}}g} \right) = \left[ {\begin{array}{ccccccccccccccc} { - 0.2 - 24\sin (2t)} \\ {0.4 + 32\sin (2t)} \\ { - 0.8 - 48\sin (2t)} \end{array}} \right]

$$     (4.8.38)
 * <p style="text-align:right">
 * }

Authors
Zongyi Yang

Jiang Jin