User:Eml5526.s11.team3/Homework 3

 Homework 3 =Problem3.1 Solve for trial solution of the PDE using Polynomial basis function with k=1=

Problem Statement
Do Problem 2.9 with polynomial basis functions with $$k=1$$: $$\{{{(x+k)}^{j}},j=0,1,2,\cdots \}$$.

For more information about the problem statement, please refer to Homework2.9.

Solution
1)

For $$n=2$$ :

The basis function is $$\displaystyle

\left\{ {{b}_{i}}\left( x \right)={{(x+1)}^{i}},i=0,1,2 \right\}

$$ then,$$b_{0}(x)=1,\quad b_{1}(x)=(x+1),\quad b_{2}(x)=(x+1)^{2}$$

the main aim of this problem is to solve the PDE given by Eq(1) by discrete Weighted Residual form method, in the method we assume the solution $$u(x)$$ as $$u^{h}(x)= \sum_{j=0}^{n}d_{j}b_{j}(x)$$

for $$n=2$$, $$u^{h}(x)= \sum_{j=0}^{2}d_{j}b_{j}(x)$$

i.e.

given essential boundary condition as $$u(1)=0$$

hence the first equation which enforces boundary condition is $$ d_{0}1+d_{1}(1+1)+d_{2}(1+1)^{2}=0 $$

Essential boundary condition is

given natural boundary condition as $$\frac{du}{dx}(0)=-4$$

natual boundary condition written as $$d_{0}\frac{db_{0}}{dx}+d_{1}\frac{db_{1}}{dx}+d_{2}\frac{db_{2}}{dx}=-4$$ at $$x=0$$

i.e $$ d_{0}0+d_{1}(1)+d_{2}2(0+1)^{2}=0 $$

Natural boundary condition is

2) We have to find one more equation to solve for $$d=\{{{d}_{j}}\}$$ by projecting the residue $$P({{u}^{h}})$$ on a basis function $${{b}_{k}}(x)$$ with k = 0, 1, 2 such that the additional equation is linearly independent from the above two equations in 1.)

hence we project the residue $$P({{u}^{h}})$$ on a basis function $${{b}_{1}}(x)$$ to obtain the third equation which is linearly independent to the two boundary conditions obtained.

we get,

given $$a_{2}=2$$,$$f(x)=3$$ and since we are projecting on first basis function $$b_{0}(x)=1$$

hence Eq (3.4) becomes $$2d_{0}\int_{0}^{1}\frac{\mathrm{d}^{2} (1)}{\mathrm{d} x^{2}}dx+2d_{1}\int_{0}^{1}\frac{\mathrm{d}^{2} (x+1)}{\mathrm{d} x^{2}}dx+2d_{2}\int_{0}^{1}\frac{\mathrm{d}^{2} (x+1)^{2}}{\mathrm{d} x^{2}}dx=-3\int_{0}^{1}dx$$

the third equation is

3) Eq.(3.2), (3.3) and (3.5) can be written in Matrix form as

Hence Clearly $$\mathbf{K}$$ is non-Symmetric matrix

4) Using Matlab, we can easily get

5) Hence $$u^{h}(x)=8-\frac{5}{2}(x+1)-\frac{3}{4}(x+1)^{2}$$

which on simplification gives

the solution of PDE Eq(3.1) gives the exact solution $$u(x)$$

Eq(3.1)is $$\frac{\mathrm{d}^{2}u }{\mathrm{d} x^{2}}=\frac{-3}{2}$$

integrating the above equation once we get

but from essential BC, $$c=-4$$

hence integrating again Eq(3.8) and applying Natural B.C we get exact solution

Matlab Code:

6) To compute $$u^{h}(x=0.5)$$ for n = 2 Substituting $$x=0.5$$ in Eq(3.7)

we get,

$$u(0.5)=2.5625$$

Hence error $${{e}_{n}}(0.5)=2.5625-2.5625=0$$

Plot $${{e}_{n}}(0.5)$$ vs n (convergence)

For $$n=4$$ and $$n=6$$ the computation becomes difficult hence Matlab is used for computational purpose.

For $$n=4$$

1) The basis function is $$\displaystyle

\left\{ {{b}_{i}}\left( x \right)={{(x+1)}^{i}},i=0,1,2,3,4 \right\}

$$ then,$$b_{0}(x)=1,b_{1}(x)=(x+1),b_{2}(x)=(x+1)^{2},b_{3}(x)=(x+1)^{3},b_{4}(x)=(x+1)^{4}$$

Two equations that enforce boundary conditions are

Essential boundary condition

and Natural boundary condition is

2)

Additional equations to solve for $$d=\{{{d}_{j}}\}$$ by projecting the residue $$P({{u}^{h}})$$ on the basis function $${{b}_{k}}(x)$$ with k = 0, 1, 2 such that the additional equations are linearly independent from the above two equations in 1.)

3) Equations in matrix form

Hence Clearly $$\mathbf{K}$$ is non-Symmetric matrix

4) Using Matlab, we can easily get

5) Hence $$u^{h}(x)=8-\frac{5}{2}(x+1)-\frac{3}{4}(x+1)^{2}$$

which on simplification gives

and exact solution

Matlab Code:

Note: It is observed that on solving for $$n=4$$ we got $$d_{3}=0$$ and $$d_{4}=0$$ hence the plot for $$n=2$$ and $$n=4$$ are same

6) Since $$d_{3}=0$$ and $$d_{4}=0$$ for $$n=4$$ hence $$u^{h}(x)$$ and $$u(x)$$ will be same for $$n=2$$ and $$n=4$$ hence

and $$u(0.5)=2.5625$$

Hence error $${{e}_{n}}(0.5)=2.5625-2.5625=0$$

For $$n=6$$

1) The basis function is $$\displaystyle

\left\{ {{b}_{i}}\left( x \right)={{(x+1)}^{i}},i=0,1,2,3,4,5,6 \right\}

$$ then,$$b_{0}(x)=1,b_{1}(x)=(x+1),b_{2}(x)=(x+1)^{2},b_{3}(x)=(x+1)^{3},b_{4}(x)=(x+1)^{4},b_{5}(x)=(x+1)^{5},b_{6}(x)=(x+1)^{6}$$

Two equations that enforce boundary conditions are

Essential boundary condition

and Natural boundary condition is

2)

Additional equations to solve for $$d=\{{{d}_{j}}\}$$ by projecting the residue $$P({{u}^{h}})$$ on the basis function $${{b}_{k}}(x)$$ with k = 0, 1, 2,3,4 such that the additional equations are linearly independent from the above two equations in 1.)

3) Equations in matrix form

Hence Clearly $$\mathbf{K}$$ is non-Symmetric matrix

4) Using Matlab, we can easily get

5) Hence $$u^{h}(x)=8-\frac{5}{2}(x+1)-\frac{3}{4}(x+1)^{2}$$

which on simplification gives

and exact solution

Plot of $$u^{h}(x) vs u(x)$$ for n=2,n=4,n=6

the approximate $$u^{h}(x)$$and exact solution $$ u(x) $$ for n=2 ,n=4 and n=6 are found to be same, which is also evident from the plot where all the plots overlap on one another

Matlab Code:

Note: It is observed that on solving for $$n=6$$ we got $$d_{3}=0$$, $$d_{4}=0$$,$$d_{5}=0$$ and $$d_{6}=0$$ hence the plot for $$n=4$$ and $$n=6$$ are same

6) Since $$d_{3}=0$$ ,$$d_{4}=0$$,$$d_{5}=0$$ and $$d_{6}=0$$ for $$n=6$$ hence $$u^{h}(x)$$ and $$u(x)$$ will be same for $$n=4$$ and $$n=6$$ hence

and $$u(0.5)=2.5625$$

Hence error $${{e}_{n}}(0.5)=2.5625-2.5625=0$$

Plot showing error between approximate and actual solution at x=0.5

since we got both approximate and actual solution to be the same hence the error is zero through out

Problem solved by Avinash.

Comment
The exact solution of this problem is a second order polynomial. And we are using the polynomial basis to approximate it. From the view of interpolation, we can get the exact solution by interpolating it using polynomial basis. And since it's the second order polynomial, and when n=2 we can get the exact solution. So, for n=4 and n=6, the approximation result is the same as n=2. That's the same as what we get for this problem.

=Problem3.2 Find global stiffness matrix, nodal displacements and reactions of the given spring element=

Problem Statement
Fish J., and Belytschko T., John Wiley & Sons Ltd, 2007, p37, problem 2.1

For the spring system given below,

a. Number the elements and nodes

b. Assemble the global stiffness and force matrix

c. Partition the system and solve for the nodal displacements

d. Compute the reaction forces



3.2(a): Elements and Nodes Numbering


The numbers within parenthesis are element numbers and the rest are node numbers. The stiffness of the elements are given. They are as follows

3.2(b) : Assembled Global Stiffness matrix and Force matrix
The stiffness matrices of individual elements are assembled according to the connectivity table to get the global stiffness matrix. Connectivity table is a one which shows the nodes associated with each element.

Element stiffness matrix is given by,

Connectivity table :

$$\begin{array}{|c|c||c|}  Element No. & Local Node 1 & Local Node 2 \\ \hline 1&1&4\\ 2&2&4\\ 3&1&3\\ 4&3&4\\ \end{array}$$

Assembled global stiffness matrix,

Since the nodes 1 and 2 are fixed, displacements at these nodes are 0. But there will be reaction forces acting at these nodes, in order to maintain static equilibrium of the system. A force of 50 N acts at node 4 and at all other nodes forces are 0. Therefore, the global force, displacement and reaction matrices are as follows,

3.2(c): Partition of the system and Nodal displacements
The system can be expressed in the equation form as follows,

The above system of equations can be partitioned as,

where,

$$ \left[\textbf{K}_E \right] = \begin{bmatrix} k_1+k_3 &0 \\ 0 &k_2  \\ \end{bmatrix} \;\;\left[\textbf{K}_{EF} \right] = \begin{bmatrix} &-k_3 &-k_1 \\ &0 &-k_2 \\ \end{bmatrix}\; \; \left[\textbf{K}^{T}_{EF} \right] = \begin{bmatrix} k_3 &0  \\ -k_1 &-k_2\\ \end{bmatrix}\; \; \left[\textbf{K}_{F} \right] = \begin{bmatrix} &k_3+k_4 &-k_4 \\ &-k_4 &k_1+k_2+k_4\\ \end{bmatrix} $$

$$ \left[\textbf{d}_{E} \right] = \begin{bmatrix} 0\\ 0\\ \end{bmatrix}\; \; \left[\textbf{d}_{F} \right] = \begin{bmatrix} u_3\\ u_4\\ \end{bmatrix}\; \; \left\{\textbf{r}_{E} \right\} = \begin{Bmatrix} r_1\\ r_2\\ \end{Bmatrix}\; \; \left\{\textbf{f}_{F} \right\} = \begin{Bmatrix} 0\\ f_4\\ \end{Bmatrix}\;$$

The reduced system of equations is given by,

Therefore,

Substituting,  $$\;\;k_1 = 3k, k_2 = k, k_3 = k, k_4 = 2k \; and \; f_4 = 50\;\;$$, we get,

$$ \begin{bmatrix} &3k &-2k \\ &-2k &6k\\ \end{bmatrix}\;\begin{bmatrix} u_3\\ u_4\\ \end{bmatrix}\; \; = \begin{Bmatrix} 0\\ 50\\ \end{Bmatrix}\; $$

$$ \Rightarrow \;\; $$ $$3ku_3 - 2ku_4 = 0;$$ $$ -2ku_3 + 6ku_4 = 50 $$

Solving the two equations, we get,

3.2(d): Reaction Forces
The reaction forces are found using the equation,

Therefore, $$ \begin{bmatrix} &-k_3 &-k_1 \\ &0 &-k_2 \\ \end{bmatrix}\; \begin{bmatrix} u_3\\ u_4\\ \end{bmatrix}\; \; = \begin{Bmatrix} r_1\\ r_2\\ \end{Bmatrix}\; $$

Substituting,  $$\;\;k_1 = 3k, k_2 = k, k_3 = k, u_3 = \frac{50}{7k}\;\; u_4 = \frac{75}{7k} \;\;$$,  we get,

$$ \begin{bmatrix} &-k &-3k \\ &0 &-k\\ \end{bmatrix}\;\begin{bmatrix} 50/7k\\ 75/7k\\ \end{bmatrix}\; \; = \begin{Bmatrix} r_1\\ r_2\\ \end{Bmatrix}\; $$ $$ \Rightarrow \;\; $$ $$\;\; \frac{-50}{7}-\frac{225}{7}= r_1\;\; ;$$ $$ \frac{-75}{7} = r_2$$

Problem solved by vnarayanan

=Problem3.3 Find global stiffness matrix, nodal displacements, reactions and stresses of bar element=

Problem Statement
Consider the truss structure given below. Nodes A and B are fixed. A force equal to 10N acts in the positive x-direction at node C. Coordinates of joints are given in meters. Young's modulus is $$E = 10^{11}$$ Pa and the cross-sectional area for all bars are $$ A = 2*10^{-2} m^{2}$$.

a. Number the elements and nodes

b. Assemble the global stiffness and force matrix

c. Partition the system and solve for nodal displacements

d. Computer the stresses and reactions

Solution
a.

The elements and nodes numbers are displayed in the figure above, with red numbers represent nodes and blue numbers stand for elements.

b.

We have already subdivided the structure into elements, and numbered the nodes and elements. For constructing the global stiffness matrix and force matrix, we take following 2 steps. Step 1, deal with the formulation of each element.

Element 1:

Element 1 numbered with global nodes 1 and 3. It’s positioned at an angle $${{\phi }^{(1)}}={{90}^{\circ }}$$ with respect to the positive horizontal x-axis. Other relations are as follow:
 * {| style="width:100%" border="0"

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

\cos ({{90}^{\circ }})=0,\quad\sin ({{90}^{\circ }})=1,\quad{{l}^{(1)}}=1m,\quad{{k}^{(1)}}=\frac{AE}=2\cdot {{10}^{9}}_ – ^ – N/m

$$ Hence we can obtain the local stiffness matrix for element 1:
 * }
 * {| style="width:100%" border="0"

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

{{K}^{(1)}}=2\cdot {{10}^{9}}\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   1  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   -1  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   -1  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   1  \\ \end{matrix}  \\ \end{matrix} \right). $$
 * }.

Element 2:

Element 2 numbered with global nodes 2 and 4. It’s also positioned at an angle $${{\phi }^{(2)}}={{90}^{\circ }}$$ with respect to the positive horizontal x-axis. Other relations are as follow:
 * {| style="width:100%" border="0"

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

\cos ({{90}^{\circ }})=0,\quad\sin ({{90}^{\circ }})=1,\quad{{l}^{(2)}}=1m,\quad{{k}^{(2)}}=\frac{AE}=2\cdot {{10}^{9}}_ – ^ – N/m

$$ Hence we can obtain the local stiffness matrix for element 2:
 * }
 * {| style="width:100%" border="0"

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

{{K}^{(2)}}=2\cdot {{10}^{9}}\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   1  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   -1  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   -1  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   1  \\ \end{matrix}  \\ \end{matrix} \right)

$$
 * }.

Element 3:

Element 3 numbered with global nodes 2 and 3. It’s positioned at an angle $${{\phi }^{(3)}}={{135}^{\circ }}$$ with respect to the positive horizontal x-axis. Other relations are as follow:
 * {| style="width:100%" border="0"

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

\cos ({{135}^{\circ }})=-\frac{\sqrt{2}}{2},\quad\sin ({{90}^{\circ }})=\frac{\sqrt{2}}{2},\quad{{l}^{(3)}}=\sqrt{2}m,\quad{{k}^{(3)}}=\frac{AE}=\sqrt{2}\cdot {{10}^{9}}_ – ^ – N/m

$$ Hence we can obtain the local stiffness matrix for element 3:
 * }
 * {| style="width:100%" border="0"

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

{{K}^{(3)}}=\sqrt{2}\cdot {{10}^{9}}\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   \frac{1}{2}  \\   -\frac{1}{2}  \\ \end{matrix}  \\   -\frac{1}{2}  \\ \end{matrix}  \\   \frac{1}{2}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   -\frac{1}{2}  \\   \frac{1}{2}  \\ \end{matrix}  \\   \frac{1}{2}  \\ \end{matrix}  \\   -\frac{1}{2}  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   -\frac{1}{2}  \\   \frac{1}{2}  \\ \end{matrix}  \\   \frac{1}{2}  \\ \end{matrix}  \\   -\frac{1}{2}  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   \frac{1}{2}  \\   -\frac{1}{2}  \\ \end{matrix}  \\   -\frac{1}{2}  \\ \end{matrix}  \\   \frac{1}{2}  \\ \end{matrix}  \\ \end{matrix} \right)

$$
 * }.

Element 4:

Element 4 numbered with global nodes 3 and 4. It’s also positioned at an angle $${{\phi }^{(4)}}={{0}^{\circ }}$$ with respect to the positive horizontal x-axis. Other relations are as follow:
 * {| style="width:100%" border="0"

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

\cos ({{0}^{\circ }})=1,\quad\sin ({{0}^{\circ }})=0,\quad{{l}^{(4)}}=1m,\quad{{k}^{(4)}}=\frac{AE}=2\cdot {{10}^{9}}_ – ^ – N/m

$$ Hence we can obtain the local stiffness matrix for element 4:
 * }
 * {| style="width:100%" border="0"

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

{{K}^{(4)}}=2\cdot {{10}^{9}}\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   \begin{matrix}   1  \\   0  \\ \end{matrix}  \\   -1  \\ \end{matrix}  \\   0  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   -1  \\   0  \\ \end{matrix}  \\   1  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\ \end{matrix} \right)

$$
 * }.

Step 2, assemble the global matrix. Direct assemble:


 * {| style="width:95%" |


 * $$\displaystyle

K=2\cdot {{10}^{9}}\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   \begin{matrix}   0 & 0  \\ \end{matrix}  \\   \begin{matrix}   0 & 1  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   0  \\   0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   0 & 0  \\ \end{matrix}  \\   \begin{matrix}   0 & 0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   0  \\   -1  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   0 & 0  \\ \end{matrix}  \\   \begin{matrix}   0 & 0  \\ \end{matrix}  \\ \end{matrix}  \\   \begin{matrix}   \begin{matrix}   \begin{matrix}   0 & 0  \\ \end{matrix}  \\   \begin{matrix}   0 & 0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \frac{\sqrt{2}}{4}  \\   -\frac{\sqrt{2}}{4}  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix}   \begin{matrix}   -\frac{\sqrt{2}}{4} & -\frac{\sqrt{2}}{4}  \\ \end{matrix}  \\   \begin{matrix}   1+\frac{\sqrt{2}}{4} & \frac{\sqrt{2}}{4}  \\ \end{matrix} \\ \end{matrix} & \begin{matrix} \frac{\sqrt{2}}{4} \\ -\frac{\sqrt{2}}{4} \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 & 0 \\ \end{matrix}  \\ \begin{matrix} 0 & 1 \\ \end{matrix}  \\ \end{matrix} \\ \begin{matrix} \begin{matrix} \begin{matrix} 0 & 0 \\ \end{matrix}  \\ \begin{matrix} 0 & -1 \\ \end{matrix}  \\ \end{matrix} & \begin{matrix} -\frac{\sqrt{2}}{4} \\ \frac{\sqrt{2}}{4} \\ \end{matrix} \\ \end{matrix} & \begin{matrix} \begin{matrix} \frac{\sqrt{2}}{4} & 1+\frac{\sqrt{2}}{4} \\ \end{matrix} \\ \begin{matrix} -\frac{\sqrt{2}}{4} & -\frac{\sqrt{2}}{4} \\ \end{matrix} \\ \end{matrix} & \begin{matrix} -\frac{\sqrt{2}}{4} \\ 1+\frac{\sqrt{2}}{4} \\ \end{matrix} & \begin{matrix} \begin{matrix} -1 & 0 \\ \end{matrix}  \\ \begin{matrix} 0 & 0 \\ \end{matrix}  \\ \end{matrix} \\ \begin{matrix} \begin{matrix} \begin{matrix} 0 & 0 \\ \end{matrix}  \\ \begin{matrix} 0 & 0 \\ \end{matrix}  \\ \end{matrix} & \begin{matrix} 0 \\   0  \\ \end{matrix}  \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 & -1 \\ \end{matrix}  \\ \begin{matrix} -1 & 0 \\ \end{matrix}  \\ \end{matrix} & \begin{matrix} 0 \\   0  \\ \end{matrix} & \begin{matrix} \begin{matrix} 1 & 0 \\ \end{matrix}  \\ \begin{matrix} 0 & 1 \\ \end{matrix}  \\ \end{matrix} \\ \end{matrix} \right) $$


 * }.

c.

Hence we have the global system as:

The global system is partitioned after four rows and four columns:
 * {| style="width:100%" border="0"

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

{{d}_{E}}=\left( \begin{matrix}  0  \\   0  \\   0  \\   0  \\ \end{matrix} \right), {{d}_{F}}=\left( \begin{matrix}   {{u}_{3x}}  \\   {{u}_{3y}}  \\   {{u}_{4x}}  \\   {{u}_{4y}}  \\ \end{matrix} \right), {{f}_{F}}=\left( \begin{matrix}   0  \\   0  \\   10  \\   0  \\ \end{matrix} \right), {{K}_{F}}=\left( \begin{matrix}   1+\frac{\sqrt{2}}{4} & -\frac{\sqrt{2}}{4} & -1 & 0  \\   -\frac{\sqrt{2}}{4} & 1+\frac{\sqrt{2}}{4} & 0 & 0  \\   -1 & 0 & 1 & 0  \\   0 & 0 & 0 & 1  \\ \end{matrix} \right), {{K}_{EF}}=\left( \begin{matrix}   0 & 0 & 0 & 0  \\   0 & 1 & 0 & 0  \\   0 & 0 & \frac{\sqrt{2}}{4} & -\frac{\sqrt{2}}{4}  \\   0 & 0 & -\frac{\sqrt{2}}{4} & 1+\frac{\sqrt{2}}{4}  \\ \end{matrix} \right)

$$ The nodal displacements are solved from the global system as:
 * }.


 * {| style="width:95%" |


 * $$\displaystyle

d=\frac{1}{2\cdot {{10}^{9}}}\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\   0  \\ \end{matrix}  \\   20\sqrt{2}+10  \\ \end{matrix}  \\   10  \\   20\sqrt{2}+20  \\   0  \\ \end{matrix} \right) $$


 * }.

d.

The force matrix also can be solved from the global system as:


 * {| style="width:95%" |


 * $$\displaystyle

F=\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   \begin{matrix}   0  \\   -10  \\ \end{matrix}  \\   -10  \\   10  \\ \end{matrix}  \\   0  \\ \end{matrix}  \\   0  \\   10  \\   0  \\ \end{matrix} \right) $$


 * }.

Stresses in the elements:

From page 34 of the textbook, Jacob Fish and Ted Belytschko, A first course in Finite Elements, we can compute the stress in elements as follow: Element 1:


 * {| style="width:95%" |


 * $$\displaystyle

{{\sigma }^{(1)}}=\frac{E}\left( \begin{matrix}  0 & -1 & 0 & 1  \\ \end{matrix} \right)\frac{1}{2\cdot {{10}^{9}}}\left( \begin{matrix}   0  \\   0  \\   20\sqrt{2}+10  \\   10  \\ \end{matrix} \right)=500Pa $$

Element 2:
 * }.


 * {| style="width:95%" |


 * $$\displaystyle

{{\sigma }^{(2)}}=\frac{E}\left( \begin{matrix}  0 & -1 & 0 & 1  \\ \end{matrix} \right)\frac{1}{2\cdot {{10}^{9}}}\left( \begin{matrix}   0  \\   0  \\   20\sqrt{2}+20  \\   0  \\ \end{matrix} \right)=0Pa $$

Element 3:
 * }.


 * {| style="width:95%" |


 * $$\displaystyle

{{\sigma }^{(3)}}=\frac{E}\left( \begin{matrix}  \frac{\sqrt{2}}{2} & -\frac{\sqrt{2}}{2} & -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}  \\ \end{matrix} \right)\frac{1}{2\cdot {{10}^{9}}}\left( \begin{matrix}   0  \\   0  \\   20\sqrt{2}+10  \\   10  \\ \end{matrix} \right)=-500\sqrt{2}Pa $$

Element 4:
 * }.


 * {| style="width:95%" |


 * $$\displaystyle

{{\sigma }^{(4)}}=\frac{E}\left( \begin{matrix}  -1 & 0 & 1 & 0  \\ \end{matrix} \right)\frac{1}{2\cdot {{10}^{9}}}\left( \begin{matrix}   20\sqrt{2}+10  \\   10  \\   20\sqrt{2}+20  \\   0  \\ \end{matrix} \right)=500Pa $$


 * }.

Problem solved by Hailong Chen.

=Problem 3.4 Show that alpha and beta are not equal=

Problem Statement
$$ \mathbf{\alpha } = b_i\left [ a_2'b_j'+a_2b_j'' \right ]$$

$$\mathbf{\beta } = b_j\left [ a_2'b_i'+a_2b_i'' \right ] $$ Where, $$\displaystyle b_i(x) = \cos(ix)$$

$$b_j(x) = \cos(jx)\;; i\;\neq\;j$$ Show that $$ \;\; \mathbf{\alpha }\;\neq\; \mathbf{\beta }$$

Solution
$$\displaystyle b_i(x) = \cos(ix)$$

$$\displaystyle b_j(x) = \cos(jx)$$ $$b_i'(x) = -i\sin(ix)\; ; \;b_i''(x) = -i^2\cos(ix) $$ $$b_j'(x) = -j\sin(jx)\; ; \;b_j''(x) = -j^2\cos(jx) $$ Therefore, $$\mathbf{\alpha } = \cos(ix)\left [ a_2'(-j\sin(jx))+a_2(-j^2\cos(jx))\right] $$ $$\mathbf{\alpha } = \cos(ix)\left [ -ja_2'\sin(jx)-j^2a_2\cos(jx)\right] $$

Therefore, $$\mathbf{\beta } = \cos(jx)\left [ a_2'(-i\sin(ix))+a_2(-i^2\cos(ix))\right] $$ $$\mathbf{\beta } = \cos(jx)\left [ -ia_2'\sin(ix)-i^2a_2\cos(ix)\right] $$

From, Equations (3.4.2) and (3.4.4), it is evident that

If $$\;\; a_2 = a_2' \neq 0 \;\; then,$$

Let us consider a specific example, say $$i = 2 \; and \; j = 1$$ Therefore, $$\displaystyle b_2(x) = \cos(2x)$$

$$\displaystyle b_1(x) = \cos(x)$$ $$b_2'(x) = -2\sin(2x)\; ; \;b_2''(x) = -4\cos(2x) $$ $$b_j'(x) = -\sin(x)\; ; \;b_j''(x) = -\cos(x) $$ $$\mathbf{\alpha } = b_2\left [ a_2'b_1'+a_2b_1'' \right ]$$ Therefore, $$\mathbf{\alpha } = \cos(2x)\left [ a_2'(-\sin(x))+a_2(-\cos(x))\right] $$ $$\mathbf{\alpha } = \cos(2x)\left [ -a_2'\sin(x)-a_2\cos(x)\right] $$ $$\mathbf{\beta} = b_1\left [ a_2'b_2+a_2b_2'' \right ]$$ Therefore, $$\mathbf{\beta } = \cos(x)\left [ a_2'(-2\sin(x))+a_2(-4\cos(2x))\right] $$ $$\mathbf{\beta } = \cos(jx)\left [ -2a_2'\sin(2x)-4a_2\cos(2x)\right] $$

From Equations (3.4.6) and (3.4.7), it is clear that,

Problem solved by vnarayanan

= Problem3.5 Show the equivalent stiffness between a bar and a horizontal aligned spring =

Problem Statement
Show that the equivalent stiffness of a spring aligned in the x-direction for the bar of thickness $$t$$ with a centered square hole shown in above figure is:
 * {| style="width:100%" border="0"

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

k=\frac{5Etab}{(a+b)l}

$$ where $$E$$ is the Young’s Modulus and $$t$$ is the thickness of the bar.
 * }

Method 1

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

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

Spring_ – ^ – Equivalence_ – ^ – :=\frac{1}=\frac{1}+\frac{1}+\frac{1}

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

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

{{k}_{1}}=\frac{{{A}^{(1)}}E}=\frac{atE}{\frac{l}{10}}

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

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

{{k}_{2}}=\frac{{{A}^{(2)}}E}=\frac{4btE}{\frac{4l}{5}}

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

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

{{k}_{3}}=\frac{{{A}^{(3)}}E}=\frac{atE}{\frac{l}{10}}

$$ From $$\frac{1}=\frac{1}+\frac{1}+\frac{1}$$, using Wolfram Alpha, we can get :
 * }
 * {| style="width:100%" border="0"

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

{{k}_{eq}}=\frac{5Etab}{(a+b)l}

$$
 * }

Wolfram Alpha Code: (1/((a*t*Y)/(L/10)) + (1/((a*t*Y)/(L/10)) + (1/((4*b*t*Y)/(L-(L/5))))
 * Note: E=Y because Wolfram Alpha recognizes E as an exponent

Problem solved by perry, edited by hylon

Picture drawn by hylon

Reference
Wolfram Alpha

Method 2


We can equivalent the bar structure into above spring structure. In order to show the stiffness of the bar is $$k=\frac{5Etab}{(a+b)l}$$, we need to find the relationship between the external force $$f$$ and the displacement of node $$4$$, as depicted in the equivalent figure above.

First, we find out the global stiffness matrix of the system.

For element 1:

Element 1 is numbered with global node 1 and 2.
 * {| style="width:100%" border="0"

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

{{k}^{(1)}}=\frac{10atE}{l}

$$ Hence we have the local stiffness for element 1:
 * }
 * {| style="width:100%" border="0"

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

{{K}^{(1)}}={{k}^{(1)}}\left( \begin{matrix}  1 & -1  \\   -1 & 1  \\ \end{matrix} \right)

$$
 * }.

For element 2:

Element 2 is numbered with global node 2 and 3.
 * {| style="width:100%" border="0"

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

{{k}^{(2)}}=\frac{5btE}{l}

$$ Hence we have the local stiffness for element 2:
 * }
 * {| style="width:100%" border="0"

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

{{K}^{(2)}}={{k}^{(2)}}\left( \begin{matrix}  1 & -1  \\   -1 & 1  \\ \end{matrix} \right)

$$
 * }.

For element 3:

Element 3 is numbered with global node 3 and 4.
 * {| style="width:100%" border="0"

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

{{k}^{(3)}}=\frac{10atE}{l}

$$ Hence we have the local stiffness for element 3:
 * }
 * {| style="width:100%" border="0"

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

{{K}^{(3)}}={{k}^{(3)}}\left( \begin{matrix}  1 & -1  \\   -1 & 1  \\ \end{matrix} \right)

$$ Direct assemble, we can obtain the global stiffness matrix as:
 * }.
 * {| style="width:100%" border="0"

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

K=\left( \begin{matrix}  {{k}^{(1)}} & -{{k}^{(1)}} & 0 & 0  \\   -{{k}^{(1)}} & {{k}^{(1)}}+{{k}^{(2)}} & -{{k}^{(2)}} & 0  \\   0 & -{{k}^{(2)}} & {{k}^{(2)}}+{{k}^{(3)}} & -{{k}^{(3)}}  \\   0 & 0 & -{{k}^{(3)}} & {{k}^{(3)}}  \\ \end{matrix} \right)

$$ Thus, we have the global system as:
 * }.
 * {| style="width:100%" border="0"

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

\left( \begin{matrix}  {{k}^{(1)}} & -{{k}^{(1)}} & 0 & 0  \\   -{{k}^{(1)}} & {{k}^{(1)}}+{{k}^{(2)}} & -{{k}^{(2)}} & 0  \\   0 & -{{k}^{(2)}} & {{k}^{(2)}}+{{k}^{(3)}} & -{{k}^{(3)}}  \\   0 & 0 & -{{k}^{(3)}} & {{k}^{(3)}}  \\ \end{matrix} \right)\left( \begin{matrix}   0  \\   {{u}_{2}}  \\   {{u}_{3}}  \\   {{u}_{4}}  \\ \end{matrix} \right)=\left( \begin{matrix}   {{r}_{1}}  \\   0  \\   0  \\   f  \\ \end{matrix} \right)

$$ Partition the global stiffness matrix after the third row and the third column, we can solve above matrix by block. And finally, we can get:
 * }
 * {| style="width:100%" border="0"

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

\frac{5Etab}{(a+b)l}{{u}_{4}}=f

$$ Hence we have the stiffness of the bar:
 * }


 * {| style="width:95%" |


 * $$\displaystyle

k=\frac{5Etab}{(a+b)l} $$


 * }

Problem solved by Hailong Chen.

=Problem3.6 Find global stiffness matrix, nodal displacements, reactions and stresses of bar element=

Problem Statement
Given the three-bar structure subjected to the prescribed load at point C equal to $$ 10^{3} $$ N as shown in Figure 2.19 on page 38. The Young's modulus is $$ E = 10^{11} $$ Pa, the cross-sectional area of the bar BC is $$ 2*10^{-2} m^{2} $$ and that of BD and BF is $$ 10^{-2} m^{2} $$. Note that point D is free to move in the x-direction. Coordinates of the joints are given in meters.

a. Construct the global stiffness matrix and load matrix

b. Partition the matrices and solve for the unknown displacements at point B and displacement in the x-direction at point D

c. Find the stresses in the three bars

d. Find the reactions at the nodes C, D, and F

Solution
a.

We have already subdivided the structure into elements, and numbered the nodes and elements in the figure attached. For constructing the global stiffness matrix and load matrix, we take following 2 steps.

Step 1, deal with the formulation of each element.

Element 1:

Element 1 numbered with global nodes 1 and 4. It’s positioned at an angle $${{\phi }^{(1)}}={{135}^{\circ }}$$ with respect to the positive x-axis. Other relations are as follow:
 * {| style="width:100%" border="0"

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

\cos ({{135}^{\circ }})=-\frac{\sqrt{2}}{2},\quad\sin ({{135}^{\circ }})=\frac{\sqrt{2}}{2},\quad{{l}^{(1)}}=\sqrt{2}m,\quad{{k}^{(1)}}=\frac{{{A}^{(1)}}E}=\frac{\sqrt{2}}

$$ Hence we can obtain the local stiffness matrix for element 1:
 * }.
 * {| style="width:100%" border="0"

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

{{K}^{(1)}}=\frac{\sqrt{2}}\left( \begin{matrix}  \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2}  \\   -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}  \\   -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}  \\   \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2}  \\ \end{matrix} \right)

$$
 * }.

Element 2:

Element 2 numbered with global nodes 2 and 4. It’s positioned at an angle $${{\phi }^{(2)}}={{90}^{\circ }}$$ with respect to the positive x-axis. Other relations are as follow:
 * {| style="width:100%" border="0"

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

\cos ({{90}^{\circ }})=0,\quad\sin ({{90}^{\circ }})=1,\quad{{l}^{(2)}}=1m,\quad{{k}^{(2)}}=\frac{{{A}^{(2)}}E}=2\cdot {{10}^{9}}

$$ Hence we can obtain the local stiffness matrix for element 2:
 * }.
 * {| style="width:100%" border="0"

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

{{K}^{(2)}}=2\cdot {{10}^{9}}\left( \begin{matrix}  0 & 0 & 0 & 0  \\   0 & 1 & 0 & -1  \\   0 & 0 & 0 & 0  \\   0 & -1 & 0 & 1  \\ \end{matrix} \right)

$$
 * }.

Element 3:

Element 3 numbered with global nodes 3 and 4. It’s positioned at an angle $${{\phi }^{(3)}}={{45}^{\circ }}$$ with respect to the positive x-axis. Other relations are as follow:
 * {| style="width:100%" border="0"

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

\cos ({{45}^{\circ }})=\frac{\sqrt{2}}{2},\quad\sin ({{45}^{\circ }})=\frac{\sqrt{2}}{2},\quad{{l}^{(3)}}=\sqrt{2}m,\quad{{k}^{(3)}}=\frac{{{A}^{(3)}}E}=\frac{\sqrt{2}}

$$ Hence we can obtain the local stiffness matrix for element 3:
 * }.
 * {| style="width:100%" border="0"

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

{{K}^{(3)}}=\frac{\sqrt{2}}\left( \begin{matrix}  \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2}  \\   \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2}  \\   -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & \frac{1}{2}  \\   -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & \frac{1}{2}  \\ \end{matrix} \right)

$$
 * }.

Step 2, assemble the global matrix.

Direct assemble:


 * {| style="width:95%" |


 * $$\displaystyle

K={{10}^{9}}\left( \begin{matrix}  \begin{matrix}   \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0  \\   -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0 & 0  \\   0 & 0 & 0 & 0  \\   0 & 0 & 0 & 2  \\ \end{matrix} & \begin{matrix}   0 & 0 & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}}  \\   0 & 0 & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   0 & 0 & 0 & 0  \\   0 & 0 & 0 & -2  \\ \end{matrix}  \\   \begin{matrix}   0 & 0 & 0 & 0  \\   0 & 0 & 0 & 0  \\   -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0 & 0  \\   \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & -2  \\ \end{matrix} & \begin{matrix}   \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & \frac{1}{\sqrt{2}} & 0  \\   -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & \frac{1}{\sqrt{2}}+2  \\ \end{matrix} \\ \end{matrix} \right) $$

We also have the load matrix as:
 * }.


 * {| style="width:95%" |


 * $$\displaystyle

F=\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   {{r}_{1x}}  \\   {{r}_{1y}}  \\ \end{matrix}  \\   {{r}_{2x}}  \\   {{r}_{2y}}  \\   0  \\ \end{matrix}  \\   {{r}_{3y}}  \\   {{10}^{3}}  \\   0  \\ \end{matrix} \right) $$


 * }.

b.

Hence we have the global system as:
 * {| style="width:100%" border="0"

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

{{10}^{9}}\left( \begin{matrix}  \begin{matrix}   \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0  \\   -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0 & 0  \\   0 & 0 & 0 & 0  \\   0 & 0 & 0 & 2  \\ \end{matrix} & \begin{matrix}   0 & 0 & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}}  \\   0 & 0 & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   0 & 0 & 0 & 0  \\   0 & 0 & 0 & -2  \\ \end{matrix}  \\   \begin{matrix}   0 & 0 & 0 & 0  \\   0 & 0 & 0 & 0  \\   -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0 & 0  \\   \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & -2  \\ \end{matrix} & \begin{matrix}   \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & \frac{1}{\sqrt{2}} & 0  \\   -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & \frac{1}{\sqrt{2}}+2  \\ \end{matrix} \\ \end{matrix} \right)\left( \begin{matrix} \begin{matrix} \begin{matrix} 0 \\   0  \\ \end{matrix}  \\ 0 \\   0  \\   {{u}_{3x}}  \\ \end{matrix} \\ 0 \\   {{u}_{4x}}  \\ {{u}_{4y}} \\ \end{matrix} \right)=\left( \begin{matrix} \begin{matrix} \begin{matrix} {{r}_{1x}} \\ {{r}_{1y}} \\ \end{matrix} \\ {{r}_{2x}} \\ {{r}_{2y}} \\ 0 \\ \end{matrix}  \\ {{r}_{3y}} \\ {{10}^{3}} \\   0  \\ \end{matrix} \right)

$$ The global system is partitioned after four rows and four columns:
 * }.
 * {| style="width:100%" border="0"

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

{{d}_{E}}=\left( \begin{matrix}  0  \\   0  \\   0  \\   0  \\ \end{matrix} \right),{{d}_{F}}=\left( \begin{matrix}   {{u}_{3x}}  \\   0  \\   {{u}_{4x}}  \\   {{u}_{4y}}  \\ \end{matrix} \right),{{f}_{F}}=\left( \begin{matrix}   0  \\   {{r}_{3y}}  \\   {{10}^{3}}  \\   0  \\ \end{matrix} \right),{{K}_{F}}=\left( \begin{matrix}   \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}  \\   -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & \frac{1}{\sqrt{2}} & 0  \\   -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & \frac{1}{\sqrt{2}}+2  \\ \end{matrix} \right),{{K}_{EF}}=\left( \begin{matrix}   \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & 0  \\   -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0 & 0  \\   0 & 0 & 0 & 0  \\   0 & 0 & 0 & 2  \\ \end{matrix} \right)

$$ The nodal displacements are solved from the global system as:
 * }.
 * {| style="width:100%" border="0"

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

d=\frac{1}\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   0  \\   0  \\ \end{matrix}  \\   0  \\   0  \\   1+2\sqrt{2}  \\ \end{matrix}  \\   0  \\   \frac{1}{2}+2\sqrt{2}  \\   \frac{1}{2}  \\ \end{matrix} \right)m

$$ The force matrix also can be solved from the global system as:
 * }.
 * {| style="width:100%" border="0"

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

F=\left( \begin{matrix}  \begin{matrix}   \begin{matrix}   -{{10}^{3}}  \\   {{10}^{3}}  \\ \end{matrix}  \\   0  \\   -{{10}^{3}}  \\   0  \\ \end{matrix}  \\   0  \\   {{10}^{3}}  \\   0  \\ \end{matrix} \right)N

$$
 * }.

c. Stresses in the three elements:

Element 1:


 * {| style="width:95%" |


 * $$\displaystyle

{{\sigma }^{(1)}}=\frac{E}\left( \begin{matrix}  \frac{\sqrt{2}}{2} & -\frac{\sqrt{2}}{2} & -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}  \\ \end{matrix} \right)\frac{1}\left( \begin{matrix}   \frac{1}{2}+2\sqrt{2}  \\   \frac{1}{2}  \\   0  \\   0  \\

\end{matrix} \right)=\sqrt{2}\cdot {{10}^{5}}Pa $$


 * }.

Element 2:


 * {| style="width:95%" |


 * $$\displaystyle

{{\sigma }^{(2)}}=\frac{E}\left( \begin{matrix}  0 & -1 & 0 & 1  \\ \end{matrix} \right)\frac{1}\left( \begin{matrix}   \frac{1}{2}+2\sqrt{2}  \\   \frac{1}{2}  \\      0  \\   0  \\ \end{matrix} \right)=-5\cdot {{10}^{4}}Pa $$


 * }

Element 3:


 * {| style="width:95%" |


 * $$\displaystyle

{{\sigma }^{(3)}}=0_ – ^ – Pa $$


 * }

d. Hence, the displacement at point B is:


 * {| style="width:95%" |


 * $$\displaystyle

{{u}_{Bx}}=(\frac{1}{2}+2\sqrt{2}){{10}^{-6}}m,{{u}_{By}}=\frac{1}{2}\cdot {{10}^{-6}}m $$

The displacement in the x-direction at point D is:
 * }.


 * {| style="width:95%" |


 * $$\displaystyle

{{u}_{Dx}}=(1+2\sqrt{2})\cdot {{10}^{-6}}m $$

The reaction at nodes C, D and F are:
 * }.


 * {| style="width:95%" |


 * $$\displaystyle

{{r}_{Cx}}=0,{{r}_{Cy}}=-{{10}^{3}}N,{{r}_{Dx}}=0,{{r}_{Dy}}=0,{{r}_{Fx}}=-{{10}^{3}}N,{{r}_{Fy}}={{10}^{3}}N $$


 * }.

Problem solved by Hailong Chen.

=Problem3.7 Solve for trial solution of the PDE using Polynomial basis function with k=0 =

Problem Statement
Do Problem 2.9 with polynomial basis functions with $$k=0$$: $$\{{{(x+k)}^{j}},j=0,1,2,\cdots \}$$.

For more information about the problem statement, please refer to Homework2.9.

Solution
This is the special case of Problem 3.1 where k=0 hence the basis functions will be $$\displaystyle

\left\{ {{b}_{i}}\left( x \right)={{x}^{i}},i=0,1,2 \right\}

$$

1)

For $$n=2$$ :

for $$n=2$$ the basis function is $$b_{0}(x)=1,b_{1}(x)=x,b_{2}(x)=x^{2}$$

applying the boundary conditions as mentioned in Problem 3.1 the two equations satisfying boundary conditions will be

Essential boundary condition is

and

Natural boundary condition is

2) One more equation which is independent to the BC's and which is obtained by projecting the residue $$P({{u}^{h}})$$ on a basis function $${{b}_{k}}(x)$$ with k = 0, 1, 2

hence the third equation is

3) Matrix form

Hence Clearly $$\mathbf{K}$$ is non-Symmetric matrix

4) Using Matlab, we can easily get

5) Hence

exact solution as solved in Problem 3.1

The approximate solution happens to match completely with the exact solution

Matlab Code:

6) $$u^{h}(x)$$ obtained in part 5 is same as the one in Problem 3.1 hence $$u^{h}(x=0.5)$$ will be same as one found in Problem 3.1

and error is $${{e}_{n}}(0.5)=2.5625-2.5625=0$$

For $$n=4$$ :

1)

for $$n=4$$ the basis function is $$b_{0}(x)=1,b_{1}(x)=x,b_{2}(x)=x^{2},b_{3}(x)=x^{3},b_{4}(x)=x^{4}$$

Two equations that enforce boundary conditions are Essential boundary condition is

and

Natural boundary condition is

2) Additional equations to solve for $$d=\{{{d}_{j}}\}$$ by projecting the residue $$P({{u}^{h}})$$ on the basis function $${{b}_{k}}(x)$$ with k = 0, 1, 2 such that the additional equations are linearly independent from the above two equations in 1.)

3) Matrix form

Hence Clearly $$\mathbf{K}$$ is non-Symmetric matrix

4) Using Matlab, we can easily get

5) Hence

exact solution

The approximate solution happens to match completely with the exact solution

Matlab Code:

6)

and error is $${{e}_{n}}(0.5)=2.5625-2.5625=0$$

For $$n=6$$ :

1)

for $$n=6$$ the basis function is $$b_{0}(x)=1,b_{1}(x)=x,b_{2}(x)=x^{2},b_{3}(x)=x^{3},b_{4}(x)=x^{4},b_{5}(x)=x^{5},b_{6}(x)=x^{6}$$

Two equations that enforce boundary conditions are Essential boundary condition is

and

Natural boundary condition is

2) Additional equations to solve for $$d=\{{{d}_{j}}\}$$ by projecting the residue $$P({{u}^{h}})$$ on the basis function $${{b}_{k}}(x)$$ with k = 0, 1, 2,3,4 such that the additional equations are linearly independent from the above two equations in 1.)

3) Matrix form

Hence Clearly $$\mathbf{K}$$ is non-Symmetric matrix

4) Using Matlab, we can easily get

5) Hence

exact solution

The approximate solution happens to match completely with the exact solution

Plot of $$u^{h}(x)vs u(x) for n=2,n=4,n=6$$



the approximate $$u^{h}(x)$$and exact solution $$ u(x) $$ for n=2 ,n=4 and n=6 are found to be same, which is also evident from the plot where all the plots overlap on one another

Matlab Code:

6)

and error is $${{e}_{n}}(0.5)=2.5625-2.5625=0$$

Plot showing error between approximate and actual solution at x=0.5

since we got both approximate and actual solution to be the same hence the error is zero through out Problem solved by Avinash.

Comment
The exact solution of this problem is a second order polynomial. And we are using the polynomial basis to approximate it. From the view of interpolation, we can get the exact solution by interpolating it using polynomial basis. And since it's the second order polynomial, and when n=2 we can get the exact solution. So, for n=4 and n=6, the approximation result is the same as n=2. That's the same as what we get for this problem.

=Problem3.8 Derive weak form from strong form=

Problem Statement
This problem is extracted Problem 3.1 from the textbook of Fish J., and Belytschko T.,

Show that the weak form of

$$\frac{d}{dx}(AE\frac{du}{dx})+2x=0 ; \quad 1<x<3 $$

$$\sigma (1)=(E\frac{du}{dx})_{x=1} = 0.1,$$

$$\displaystyle u(3)=0.001,$$

is given by

$$\int\limits_{1}^{3}{\frac{dw}{dx}AE\frac{du}{dx}\text{dx}=-0.1{{(wA)}_{x=1}}}$$ +$$\int_{1}^{3}{2xw}\text{dx}$$

Solution
The strong form of the governing differential equation is:-
 * {| style="width:100%" border="0"

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

\frac{d}{dx}(AE\frac{du}{dx})+2x=0 ; 1)
 * 
 * }

The boundary conditions are given by the equations:-
 * {| style="width:100%" border="0"

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

\sigma (1)={{(E\frac{du}{dx})}_{x=1}}=0.1

$$ (Eq<2> )
 * 
 * }


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

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

u(3) = 0.001

$$ (Eq<3> )
 * 
 * }

To form the weak form, we first consider the strong form and multiply the governing equation Eq(2) by an arbitrary function $$\displaystyle w(x)$$ and integrate over the domain on which they hold.

The function $$\displaystyle w(x)$$ is called the weight function.

For the governing equation, the pertinent domain is the interval [1,3].

The resulting equation after multiplying the governing equation with $$\displaystyle w(x)$$ and integrating in the domain takes the form -


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

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

\int\limits_{1}^{3}{w[\frac{d}{dx}(AE\frac{du}{dx})+2x]}\text{dx} = 0 ; \forall w

$$ (Eq<4> )
 * 
 * }

In the above, $$\forall w$$ denotes that $$\displaystyle w(x)$$ is an arbitrary function which holds for all functions $$\displaystyle w(x)$$. Furthermore, it is convenient to consider that the weight function satisfies the condition


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

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

w(3)=0

$$ (Eq<5> )
 * 
 * }

For convenience, we write Eq(4) into the equivalent form :


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

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

\int\limits_{1}^{3}{w[\frac{d}{dx}(AE\frac{du}{dx})]\text{dx}+\int\limits_{1}^{3}{2wx\text{dx}=0}}

$$ (Eq<6> )
 * 
 * }

To the solve the first term of the above integral we use Integration by parts method:

The stress $$ \sigma\left (x\right)$$ can be expressed in terms of elastic modulus $$ E\left (x\right)$$ and strain$$(\frac{du}{dx})$$ as


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

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

\sigma (x)=E(x)\frac{du}{dx}

$$ (Eq<9> )
 * 
 * }

Using Eq(9), Eq(8) is written as


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

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

{{(wA\sigma )}_{x=3}}-{{(wA\sigma )}_{x=1}}- \int\limits_{1}^{3}{\frac{dw}{dx}(AE\frac{du}{dx})\text{dx}} +\int\limits_{1}^{3}{2wx\text{dx}=0}

$$ (Eq<10> )
 * 
 * }

Now   $${{(w)}_{x=3}}=0\quad and\quad{{(\sigma )}_{x=1}}=0.1$$ as defined earlier.

Applying these conditions to Eq(10), we obtain


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

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

-0.1{{(wA)}_{x=1}} - \int\limits_{1}^{3}{\frac{dw}{dx}(AE\frac{du}{dx})\text{dx}} +\int\limits_{1}^{3}{2wx\text{dx}=0} $$ (Eq<11> )
 * 
 * }

Rearranging the terms in Eq(11) we obtain the strong form which is given by:-

Problem solved by Ushnish.

2)
We know that To check the equilibrium equation in the strong form in Problem 3.1, we need to substitute Eq 10.14 into Eq 10.1a

3)
We know that the natural boundary condition is

To check the natural boundary condition; we need to substitute Eq 10.14 into Eq 10.1b

Problem solved by sahin and Ushnish

Reference
Matlab