text:
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

title: User:Eml5526.s11.team3/Homework 3


text:
User:Eml5526.s11.team3/Homework 4

 Homework 4

=Problem4.1 Solve for trial solution of the PDE using sine basis=

Given
$${{b}_{i}}(x)=[1,\sin (x),\sin (2x),\cdots ]$$.

Objective
$$1).$$ Find two equations that enforce boundary conditions for $${{u}^{h}}(x)=\sum\limits_{j=0}^{n}{{{d}_{j}}{{b}_{j}}(x)}$$.

Natural boundary condition: $$-\frac{d{{u}^{h}}(x=0)}{dx}=4$$;

essential boundary condition: $${{u}^{h}}(x=1)=0$$.

$$2).$$ Find one more equation to solve for $$\underline{d}=\left\{ {{d}_{j}} \right\} (j=0,1,2)$$by project the residue $$P({{u}^{h}})$$ on a basis function $${{b}_{k}}(x)$$ with $$k=0,1,2 $$, such that the additional equation is linear independent from the above two equations in $$(1)$$. $$3).$$ Display three equations in matrix form $$\underline{K}\underline{d}=\underline{F}$$, observe symmetric property of $$\underline{K}$$. $$4).$$ Solve for $$\underline{d}$$. $$5).$$ Construct $$u_{2}^{h}(x)$$ and plot $$u_{2}^{h}(x)$$ vs. $$u_{2}^ – (x)$$. $$6).$$Compute $$u_{2}^{h}(x=0.5)$$ and Error $${{e}_{2}}(x=0.5)=u(x=0.5)-{{u}^{h}}(x=0.5)$$. $$7).$$ Repeat $$1)-6)$$ for $$n=4,n=6$$. $$8).$$ Plot $${{e}_{2}}(x=0.5)$$ vs. $$n$$ Refer to lecture slide [[media:fe1.s11.mtg18.djvu|18-2]] for more information.

When n =2
1.

Equation enforcing natural boundary 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" |

{{d}_{1}}+2{{d}_{2}}=-4

$$
 * }

Equation enforcing essential boundary 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" |

{{d}_{0}}+0.8415{{d}_{1}}+0.9093{{d}_{2}}=0

$$
 * }

2.

Equation from projection:

 * {| 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.9194{{d}_{1}}+5.6646{{d}_{2}}=3

$$
 * }

3.

Display three equations in matrix 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" |

\left( \begin{matrix}  0 & 1 & 2  \\   1 & 0.8415 & 0.9093  \\   0 & 0.9194 & 5.6646  \\ \end{matrix} \right)\left( \begin{matrix}   {{d}_{0}}  \\   {{d}_{1}}  \\   {{d}_{2}}  \\ \end{matrix} \right)=\left( \begin{matrix}   -4  \\   0  \\   3  \\ \end{matrix} \right)

$$
 * }

$$K$$ is non-symmetric.

4.

 * {| 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}  {{d}_{0}}  \\   {{d}_{1}}  \\   {{d}_{2}}  \\ \end{matrix} \right)={{\left( \begin{matrix}   0 & 1 & 2  \\   1 & 0.8415 & 0.9093  \\   0 & 0.9194 & 5.6646  \\ \end{matrix} \right)}^{-1}}\left( \begin{matrix}   -4  \\   0  \\   3  \\ \end{matrix} \right)=\left( \begin{matrix}   4.7162  \\   -7.4908  \\   1.7454  \\ \end{matrix} \right)

$$
 * }

5.

The approximate solution:

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

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

\text{u}_{2}^{h}\text{(x)=4}\text{.7162-7}\text{.4908sin(x)+1}\text{.7454sin(2x)}

$$
 * }

6.

Error between exact solution and approximate solution at point x=0.5.

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

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

error(x=0.5)=-0.0311

$$
 * }

When n =4
1.

Equation enforcing natural boundary 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" |

{{d}_{1}}+2{{d}_{2}}+3{{d}_{3}}+4{{d}_{4}}=-4

$$
 * }

Equation enforcing essential boundary 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" |

{{d}_{0}}+0.8415{{d}_{1}}+0.9093{{d}_{2}}+0.1411{{d}_{3}}-0.7568{{d}_{4}}=0

$$
 * }

2.

Equations from projection:

 * {| 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.9194{{d}_{1}}+\text{5}.\text{6646}{{d}_{2}}+\text{11}.\text{9400}{{d}_{3}}+\text{13}.\text{2291}{{d}_{4}}=3

$$
 * }

 * {| 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.\text{5454}{{d}_{1}}+\text{3}.\text{1777}{{d}_{2}}+\text{5}.\text{7946}{{d}_{3}}+\text{3}.\text{8212}{{d}_{4}}=\text{1}.\text{379}0\text{9}

$$
 * }

 * {| 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.\text{7944}{{d}_{1}}+\text{4}.\text{7568}{{d}_{2}}+\text{9}.\text{2993}{{d}_{3}}+\text{8}.0\text{195}{{d}_{4}}=\text{2}.\text{1242}

$$
 * }

3.

Display five equations in matrix 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" |

\left( \begin{matrix}  0 & 1 & 2 & 3 & 4  \\   1 & 0.\text{8415} & 0.\text{9}0\text{93} & 0.\text{1411} & -0.\text{7568}0  \\   0 & 0.\text{9194} & \text{5}.\text{6646} & \text{11}.\text{9400} & \text{13}.\text{2291}  \\   0 & 0.\text{5454} & \text{3}.\text{1777} & \text{5}.\text{7946} & \text{3}.\text{8212}  \\   0 & 0.\text{7944} & \text{4}.\text{7568} & \text{9}.\text{2993} & \text{8}.0\text{195}  \\ \end{matrix} \right)\left( \begin{matrix}   {{d}_{0}}  \\   {{d}_{1}}  \\   {{d}_{2}}  \\   {{d}_{3}}  \\   {{d}_{4}}  \\ \end{matrix} \right)=\left( \begin{matrix}   -4  \\   0  \\   3  \\   \text{1}.\text{379}1  \\   \text{2}.\text{1242}  \\ \end{matrix} \right)

$$
 * }

$$K$$ is non-symmetric.

4.

 * {| 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}  {{d}_{0}}  \\   {{d}_{1}}  \\   {{d}_{2}}  \\   {{d}_{3}}  \\   {{d}_{4}}  \\ \end{matrix} \right)={{\left( \begin{matrix}   0 & 1 & 2 & 3 & 4  \\   1 & 0.\text{8415} & 0.\text{9}0\text{93} & 0.\text{1411} & -0.\text{7568}  \\   0 & 0.\text{9194} & \text{5}.\text{6646} & \text{11}.\text{9400} & \text{13}.\text{2291}  \\   0 & 0.\text{5454} & \text{3}.\text{1777} & \text{5}.\text{7946} & \text{3}.\text{8212}  \\   0 & 0.\text{7944} & \text{4}.\text{7568} & \text{9}.\text{2993} & \text{8}.0\text{1949}  \\ \end{matrix} \right)}^{-1}}\left( \begin{matrix}   -4  \\   0  \\   3  \\   \text{1}.\text{379}1  \\   \text{2}.\text{1242}  \\ \end{matrix} \right)=\left( \begin{matrix}   \text{4}.\text{7491}  \\   -\text{13}.\text{1144}  \\   \text{8}.\text{2}0\text{87}  \\   -\text{3}.\text{6116}  \\   0.\text{8830}  \\ \end{matrix} \right)

$$
 * }

5.

The approximate solution:

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

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

\text{u}_{4}^{h}\text{(x)=4}.\text{7491}-\text{ 13}.\text{1144sin}\left( \text{x} \right)\text{ }+\text{8}.\text{2}0\text{87sin}\left( \text{2x} \right)\text{ }-\text{ 3}.\text{6116sin}\left( \text{3x} \right)\text{ }+\text{ }0.\text{8830sin}\left( \text{4x} \right)

$$
 * }

6.

Error between exact solution and approximate solution at point x=0.5.

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

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

\text{ error(x=0}\text{.5)=0}\text{.0069}

$$
 * }

When n =6
1.

Equation enforcing natural boundary 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" |

1{{d}_{1}}+2{{d}_{2}}+3{{d}_{3}}+4{{d}_{4}}+5{{d}_{5}}+6{{d}_{6}}=-4

$$
 * }

Equation enforcing essential boundary 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" |

{{d}_{0}}+0.8415{{d}_{1}}+0.9093{{d}_{2}}+0.1411{{d}_{3}}-0.7568{{d}_{4}}-0.9589{{d}_{5}}-0.2794{{d}_{6}}=0

$$
 * }

2.

Equations from projection:

 * {| 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.9194{{d}_{1}}+5.6646{{d}_{2}}+11.9400{{d}_{3}}+13.2291{{d}_{4}}+7.1634{{d}_{5}}+0.4780{{d}_{6}}=3

$$
 * }

 * {| 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.5454{{d}_{1}}+3.1777{{d}_{2}}+5.7946{{d}_{3}}+3.8212{{d}_{4}}-3.5658{{d}_{5}}-10.2830{{d}_{6}}=1.3791

$$
 * }

 * {| 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.7944{{d}_{1}}+4.7568{{d}_{2}}+9.2993{{d}_{3}}+8.0195{{d}_{4}}-1.1704{{d}_{5}}-11.2633{{d}_{6}}=2.1242

$$
 * }

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

$$\displaystyle 0.6438{{d}_{1}}+4.1330{{d}_{2}}+9.4191{{d}_{3}}+11.9619{{d}_{4}}+8.2745{{d}_{5}}+0.0450{{d}_{6}}=1.9900
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$
 * }

 * {| 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.2388{{d}_{1}}+2.0049{{d}_{2}}+6.7285{{d}_{3}}+14.0213{{d}_{4}}+19.8920{{d}_{5}}+18.3258{{d}_{6}}=1.2402

$$
 * }

3.

Display seven equations in matrix form:

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

$$\displaystyle \left( \begin{matrix}  0 & 1 & 2 & 3 & 4 & 5 & 6  \\   1 & 0.8415 & 0.9093 & 0.1411 & -0.7568 & -0.9589 & -0.2794  \\   0 & 0.9194 & 5.6646 & 11.9400 & 13.2291 & 7.1634 & 0.4780  \\   0 & 0.5454 & 3.1777 & 5.7946 & 3.8212 & -3.5658 & -10.2830  \\   0 & 0.7944 & 4.7568 & 9.2993 & 8.0195 & -1.1704 & -11.2633  \\   0 & 0.6438 & 4.1330 & 9.4191 & 11.9619 & 8.2745 & 0.0450  \\   0 & 0.2388 & 2.0049 & 6.7285 & 14.0213 & 19.8920 & 18.3258  \\ \end{matrix} \right)\left( \begin{matrix}   {{d}_{0}}  \\   {{d}_{1}}  \\   {{d}_{2}}  \\   {{d}_{3}}  \\   {{d}_{4}}  \\   {{d}_{5}}  \\   {{d}_{6}}  \\ \end{matrix} \right)=\left( \begin{matrix}   -4  \\   0  \\   3  \\   1.3791  \\   2.1242  \\   1.9900  \\   1.2402  \\ \end{matrix} \right)
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$
 * }

$$K$$ is non-symmetric.

4.

 * {| 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}  {{d}_{0}}  \\   {{d}_{1}}  \\   {{d}_{2}}  \\   {{d}_{3}}  \\   {{d}_{4}}  \\   {{d}_{5}}  \\   {{d}_{6}}  \\ \end{matrix} \right)={{\left( \begin{matrix}   0 & 1 & 2 & 3 & 4 & 5 & 6  \\   1 & 0.8415 & 0.9093 & 0.1411 & -0.7568 & -0.9589 & -0.2794  \\   0 & 0.9194 & 5.6646 & 11.9400 & 13.2291 & 7.1634 & 0.4780  \\   0 & 0.5454 & 3.1777 & 5.7946 & 3.8212 & -3.5658 & -10.2830  \\   0 & 0.7944 & 4.7568 & 9.2993 & 8.0195 & -1.1704 & -11.2633  \\   0 & 0.6438 & 4.1330 & 9.4191 & 11.9619 & 8.2745 & 0.0450  \\   0 & 0.2388 & 2.0049 & 6.7285 & 14.0213 & 19.8920 & 18.3258  \\ \end{matrix} \right)}^{-1}}\left( \begin{matrix}   -4  \\   0  \\   3  \\   1.3791  \\   2.1242  \\   1.9900  \\   1.2402  \\ \end{matrix} \right)=\left( \begin{matrix}   4.7500  \\   -257.2133  \\   341.5391  \\   -261.6497  \\   128.7012  \\   -38.5828  \\   5.5312  \\ \end{matrix} \right)

$$
 * }

5.

The approximate solution:

 * {| 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_{6}^{h}(x)=4.75-257.2133\sin (x)+341.5391\sin (2x)-261.6497\sin (3x)+128.7012\sin (4x)-38.5828\sin (5x)+5.5312\sin (6x)

$$
 * }

6.

Error between exact solution and approximate solution at point x=0.5.

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

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

\text{ error(x=0}\text{.5)=0}\text{.0086}

$$
 * }

Plot of Error(x=0.5) Vs. n

=Problem 4.2 Deriving of the Weak Form for One Dimensional Heat Conduction =

Given: Equations of heat conduction in one dimension
A first Course in Finite Elements, Fish & Belychtschko, John Wiley & Sons Ltd, 2007, p.73, problem 3.5

Strong form for 1D heat conduction problems

Essential boundary conditions, Natural boundary conditions, The condition on the right is a convection condition

Find
Derive the weak form for one dimensional heat conduction.

Solution
We multiply the first equation in the strong form (Eq 2.1) by the weight function and integrate over the domain on $$\Omega$$, which yields

Using Integration by Parts of the first term in (Eq 2.4) gives

Since $$w$$ is arbitrary, we select $$w$$ such that $$w(0)=0$$ to drop out the unknown term at x=0. As we remember we already have essential boundary conditions at x=0 and substitute natural boundary condition (Eq 2.3) into (Eq 2.6). So we get

If we rearrange the equation, we obtain the weak form

=Problem 4.3 Derive weak form from given strong form, solve for trial solution and the differential equation respectively=

Statement
A first Course in Finite Elements, Fish & Belychtschko, John Wiley & Sons Ltd, 2007, p.73, problem 3.6

Given
The Strong form for the heat conduction problem in a circular plate is

Where R is the 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, and k,s and R are given

a) Construct the weak form for the given Strong form
Construct the weak form for the above strong form.

b) Solution to the weak form
Use the quadratic trial (candidate) solution of the form $$\alpha _{0}+\alpha _{1}r+\alpha _{2}r^{2}$$ and weight function of the same form to obtain the solution of the weak form.

c) Solution to the strong form
Solve the differential equation with the boundary conditions ans show that the temperature distribution along the radius is given by

a) Weak form of the Strong form for the heat conduction in a circular plate
Multiplying the given strong from with the weight function w and integrating it over the interval $$0< r\leq R$$ we get,

Integrating first term of Equation 4.3.3 by parts:

We know that the weight function should vanish at $$r=R$$ i.e it should satisfy the homogeneous essential boundary condition $$(w(R)=0)$$ and from the natural boundary condition $$\frac{\mathrm{d} T}{\mathrm{d} r}(r=0)=0$$

hence substituting the above boundary conditions and rearranging the terms we get,

hence the weak form for the strong form in Eq.4.3.1 is

b) Solution of the weak form using quadratic trial solution
Given the trial solution is of the form $$T=\alpha _{0}+\alpha _{1}r+\alpha _{2}r^{2}$$ and the weight function is also of the same form and can be written as $$w=\beta _{0}+\beta _{1}r+\beta _{2}r^{2}$$

From the given essential boundary condition we can write,

differentiating the trial solution of T we get

given natural boundary condition as $$\frac{\mathrm{d} T}{\mathrm{d} r}(r=0)=0$$ hence, $$\alpha_{1}=0$$ hence the trial solution is such that

from the homogeneous essential boundary condition we know $$w(R)=0$$ hence $$\beta _{0}+\beta _{1}R+\beta _{2}R^{2}=0$$ hence the weighting function is such that

substituting Eq 4..9 and 4.3.10 into the weak form derived in Eq 4.3.6 we get,

Integrating the left hand side of Eq 4.3.11 we get

Integrating the right hand side of Eq 4.3.11 we get

Hence combining Eq 4.3.13 and Eq 4.3.12 we can write,

Rearranging the terms in the above equation we get,

The above equation should be true for any value of $$\beta_{1}$$ and $$\beta_{2}$$ hence equating each of $$\beta_{1}$$ and $$\beta_{2}$$ to zero we get, both the equations yield the same result

Substituting $$\alpha_{2}$$ into Eq 4.3.7 we get, $$\alpha_{0}=\frac{sR^{2}}{4k}$$ Hence the trial solution to the weak form in Eq 4.3.6 is

c) Solution of the strong form
The Strong form for the heat conduction problem in a circular plate from Eq 4.3.1 is

Rearranging the terms in the strong form we get,

integrating the above Eq we get,

using the natural boundary condition we get, $$c_{1}=0$$ Agian integrating Eq 4.3.19 we get,

Applying essential boundary condition $$(T(r=R)=0)$$ we get, $$c_{2}=\frac{R^{2}s}{4k}$$

Hence the temperature distribution of the circular plate along the radius is given by

Note:- From Eq 4.3.17 and Eq 4.3.21 it is observed that the solution to the weak form and the strong form are the same in this case

=Problem 4.4 Construct weak form for the torsion bar and find integration constant using B.C.=

Statement
A first Course in Finite Elements, Fish & Belychtschko, John Wiley & Sons Ltd, 2007, p.73, problem 3.7

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

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

\frac{d}{dx}(JG\frac{d\phi }{dx})+m=0,_ – ^ – 0\le x\le l

$$
 * }

Natural boundary condition:

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

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

M(x=l)={{(JG\frac{d\phi }{dx})}_{l}}=\bar{M},

$$
 * }

Essential boundary condition:

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

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

\phi (x=0)=\bar{\phi }.

$$
 * }

Where $$m(x)$$ is a distributed moment per unit length, $$M$$ is the torsion moment, $$\phi $$ is the angle of rotation, $$G$$ is the shear modulus and $$J$$ is the polar moment of inertial given by $$J=\pi {{C}^{4}}/2$$, where $$C$$ is the radius of the circular shaft.

Construct the weak form
As the weight function must vanish on the essential boundaries, we consider all smooth weight functions $$w(x)$$ such that $$w(0)= 0$$.

Multiplying the given differential equation and the natural boundary condition over the domain specified, by an arbitrary weight function:

Replacing the coefficients and substituting the modified boundary conditions:

Now we integrate the Equation 4.4.3 by parts,

Integration by parts:

We have constructed the weight functions so that w(0) =0 therefore, the first term on the RHS of the above vanishes at x = 0. Substituting Equation 4.4.5 in Equation 4.4.1,

Substituting Equation 4.4.4 in Equation 4.4.6,

 * {| 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_{0}^{1}(3x+2)\frac{\mathrm{dw} }{\mathrm{d} x}\frac{\mathrm{d} \phi }{\mathrm{d} x}dx - \int_{0}^{1}5xw(x)dx\;\; - 12(w)_{x=1} = 0\;\; \forall\;\;w(x)\;\;with\;\;w(0)=0

$$
 * }

Determinate the integration constants
Assuming m(x) = 0, reduces the given strong form as follows:

Substituting for the coefficient,

Integrating with respect to x,

From Equation 4.4.4 ,

Hence,

Integrating Equation 4.4.9 with respect to x,

From Equation 4.4.2,

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

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

\phi = 4log(3x+2)+1.228

$$
 * }

=Problem 4.5=

Given
For each of the family of basis functions, and for $$\displaystyle{{\Gamma }_{g}}=\{\beta \}$$, find the corresponding basis functions satisfying

{ $$\displaystyle {{b}_{j}}(x), j=0,1,2...n $$ } such that $$\displaystyle {{b}_{0}}(\beta)\ne 0$$   and   $$\displaystyle {{b}_{j}}(\beta)= 0, j = 1,2,3...n $$

Polynomial Basis Function
$$\displaystyle i> {{F}_{p}}={{x}^{j}}, j=0,1,2...n $$

Cosine Basis Function
$$\displaystyle ii> {{F}_{c}}=\cos jx, j=0,1,2...n $$

Sine Basis Function
$$\displaystyle iii> {{F}_{s}}=1,\sin jx, j=1,2...n $$

Fourier Basis Function
$$\displaystyle iv> {{F}_{f}}=\cos jx,\sin kx, j=0,1,2...n , k=1,2...n$$

Exponential Basis Function
$$\displaystyle v> {{F}_{e}}={{e}^{jx}}, j=0,1,2...n $$

$$\Omega=\left[\alpha,\beta\right]=\left[-2,4\right]$$ is the domain where the family of basis functions satisfy the Constraint Breaking Solution.

Polynomial Basis Functions
$$\displaystyle i> $$ For the Constraint Breaking Function to be satisfied, we take the basis function for the Polynomial case so that it satisfies the Essential Boundary Condition $$\Gamma _{g}=\left \{ \beta\right \}$$ 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" |

\bar{{F}_{p}}={{(x-\beta)}^{j}}, j=0,1,2...n

$$ (Eq )<1>
 *
 * }

For $$\displaystyle j=0 $$ at  $$\displaystyle x= \beta, \bar{{F}_{p}} = {(x-\beta)^0} = 1$$

For $$\displaystyle j=1 $$ at  $$\displaystyle x= \beta, \bar{{F}_{p}} = {(\beta-\beta)^1} = 0$$

For $$\displaystyle j=2 $$ at  $$\displaystyle x= \beta, \bar{{F}_{p}} = {(\beta-\beta)^2} = 0$$

For $$\displaystyle j=n $$ at  $$\displaystyle x= \beta, \bar{{F}_{p}} = {(\beta-\beta)^n} = 0$$

Hence the Constraint Breaking Function is satisfied by this polynomial basis function where

$$\displaystyle {{b}_{0}}(\beta)\ne 0$$  and   $$\displaystyle {{b}_{j}}(\beta)= 0, j = 1,2,3...n $$

Now for the given problem, the value of $$\displaystyle \beta =4 $$

The plot of $$\displaystyle {{F}_{p}}$$ and $$\displaystyle\bar{{F}_{p}} $$ is generated from Matlab.

Matlab Code:

Cosine Basis Functions
$$\displaystyle ii>$$ While considering the Cosine basis function $$\displaystyle {{F}_{c}}=\cos jx $$, we are unable to find such a transformation for which $$\displaystyle {{b}_{0}}(\beta)= 0$$.

So we make the transformation of
 * {| style="width:100%" border="0"

$$\displaystyle {{b}_{0}}(\beta)= 1, \bar{{F}_{c}}= \cos (2j-1)\pi x/2\beta, j=1,2...n
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ (Eq )<2>
 *
 * }

For $$\displaystyle j=1 $$ at  $$\displaystyle x= \beta, \bar{{F}_{c}}= \cos (2-1)\pi \beta/2\beta= \cos \pi/2 = 0$$

For $$\displaystyle j=2 $$ at  $$\displaystyle x= \beta, \bar{{F}_{c}}= \cos (4-1)\pi \beta/2\beta= \cos 3\pi/2 = 0$$

For $$\displaystyle j=3 $$ at  $$\displaystyle x= \beta, \bar{{F}_{c}}= \cos (6-1)\pi \beta/2\beta= \cos 5\pi/2 = 0$$

For $$\displaystyle j=n $$ at  $$\displaystyle x= \beta, \bar{{F}_{c}}= \cos (2n-1)\pi \beta/2\beta= \cos (2n-1)\pi/2 = 0$$

Hence the Constraint Breaking Function is satisfied by this Cosine basis function where

$$\displaystyle {{b}_{0}}(\beta)\ne 0$$  and   $$\displaystyle {{b}_{j}}(\beta)= 0, j = 1,2,3...n $$

Now for the given problem, the value of $$\displaystyle \beta =4 $$

The plot of $$\displaystyle {{F}_{c}}$$ and $$\displaystyle\bar{{F}_{c}} $$ is generated from Matlab.

Matlab Code:

Sine Basis Functions
$$\displaystyle iii>$$ While considering the sine basis function $$\displaystyle {{F}_{s}}=\sin jx $$, we have taken the  transformation  $$\displaystyle {{b}_{0}}(\beta)= 1$$.

So we make the transformation for the sin function as of
 * {| style="width:100%" border="0"

$$\displaystyle {{b}_{0}}(\beta)= 1, \bar{{F}_{s}}= \sin (j\pi x/\beta), j=1,2...n
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$ (Eq )<3>
 *
 * }

For $$\displaystyle j=1 $$ at  $$\displaystyle x= \beta, \bar{{F}_{s}}= \sin 1\pi \beta/\beta= \sin \pi = 0$$

For $$\displaystyle j=2 $$ at  $$\displaystyle x= \beta, \bar{{F}_{s}}= \sin 2\pi \beta/\beta= \sin 2\pi = 0$$

For $$\displaystyle j=3 $$ at  $$\displaystyle x= \beta, \bar{{F}_{s}}= \sin 3\pi \beta/\beta= \sin 3\pi = 0$$

For $$\displaystyle j=n $$ at  $$\displaystyle x= \beta, \bar{{F}_{s}}= \sin n\pi \beta/\beta= \sin n\pi = 0$$

Hence the Constraint Breaking Function is satisfied by this Sine basis function where

$$\displaystyle {{b}_{0}}(\beta)\ne 0$$  and   $$\displaystyle {{b}_{j}}(\beta)= 0, j = 1,2,3...n $$

Now for the given problem, the value of $$\displaystyle \beta =4 $$

The plot of $$\displaystyle {{F}_{s}}$$ and $$\displaystyle \bar{{F}_{s}} $$ is generated from Matlab.

Matlab Code:

Fourier Basis Functions
$$\displaystyle iv>$$ While considering the Fourier basis function $$\displaystyle {{F}_{f}}=\cos jx,\sin kx $$, we are unable to find such a transformation for which $$\displaystyle {{b}_{0}}(\beta)= 0$$.

So we make the transformation of

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

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

{{b}_{0}}(\beta)= 1, \bar{{F}_{f}}= \cos (2j-1)\pi x/2\beta, j=1,2...n , \sin (k\pi x/\beta), k=1,2...n

$$ (Eq )<4>
 *
 * }

For $$\displaystyle j=1 $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \cos (2-1)\pi \beta/2\beta= \cos \pi/2 = 0$$

For $$\displaystyle k=1 $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \sin 1\pi \beta/\beta= \sin \pi = 0$$

For $$\displaystyle j=2 $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \cos (4-1)\pi \beta/2\beta= \cos 3\pi/2 = 0$$

For $$\displaystyle k=2 $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \sin 2\pi \beta/\beta= \sin 2\pi = 0$$

For $$\displaystyle j=3 $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \cos (6-1)\pi \beta/2\beta= \cos 5\pi/2 = 0$$

For $$\displaystyle k=3 $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \sin 3\pi \beta/\beta= \sin 3\pi = 0$$

For $$\displaystyle j=n $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \cos (2n-1)\pi \beta/2\beta= \cos (2n-1)\pi/2 = 0$$

For $$\displaystyle k=n $$ at  $$\displaystyle x= \beta, \bar{{F}_{f}}= \sin n\pi \beta/\beta= \sin n\pi = 0$$

Hence the Constraint Breaking Function is satisfied by this Fourier basis function where

$$\displaystyle {{b}_{0}}(\beta)\ne 0$$  and   $$\displaystyle {{b}_{j}}(\beta)= 0, j = 1,2,3...n $$

Now for the given problem, the value of $$\displaystyle \beta =4 $$

The plot of $$\displaystyle {{F}_{f}}$$ and $$\displaystyle\bar{{F}_{f}} $$ is generated from Matlab.

Matlab Code:

Exponential Basis Functions
$$\displaystyle v>$$ While considering the Exponential basis function $$\displaystyle {{F}_{e}}={{e}^{jx}}$$, we are unable to find such a transformation for which $$\displaystyle {{b}_{0}}(\beta)= 0$$.

So we make the transformation of

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

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

{{b}_{0}}(\beta)= 1, \bar{{F}_{e}}= {{e}^{j(x-\beta)}}-1, j=1,2...n

$$ (Eq )<5>
 *
 * }

For $$\displaystyle j=1 $$ at  $$\displaystyle x= \beta, \bar{{F}_{e}}= {{e}^{1(\beta-\beta)}}-1= {e}^0 - 1 = 1-1 = 0$$

For $$\displaystyle j=2 $$ at  $$\displaystyle x= \beta, \bar{{F}_{e}}= {{e}^{2(\beta-\beta)}}-1= {e}^0 - 1 = 1-1 = 0$$

For $$\displaystyle j=3 $$ at  $$\displaystyle x= \beta,\bar{{F}_{e}}= {{e}^{3(\beta-\beta)}}-1= {e}^0 - 1 = 1-1 = 0$$

For $$\displaystyle j=n $$ at  $$\displaystyle x= \beta,\bar{{F}_{e}}= {{e}^{n(\beta-\beta)}}-1= {e}^0 - 1 = 1-1 = 0$$

Hence the Constraint Breaking Function is satisfied by this polynomial basis function where

$$\displaystyle {{b}_{0}}(\beta)\ne 0$$  and   $$\displaystyle {{b}_{j}}(\beta)= 0, j = 1,2,3...n $$

Now for the given problem, the value of $$\displaystyle \beta =4 $$

The plot of $$\displaystyle {{F}_{e}}$$ and $$\displaystyle\bar{{F}_{e}} $$ is generated from Matlab.

Matlab Code:

2).
$$\displaystyle {{F}_{e}}={{e}^{jx}}, j=0,1,2 $$

F = {$$\displaystyle 1,e^x,e^{2x}$$} = {$$\displaystyle b_{1}(x),b_{2}(x),b_{3}(x)$$} $$\Gamma = \begin{bmatrix} \Gamma_{11}&\Gamma_{12}&\Gamma_{13} \\ \Gamma_{21}&\Gamma_{22}&\Gamma_{23} \\ \Gamma_{31}&\Gamma_{32}&\Gamma_{33} \\ \end{bmatrix}$$

$$\;\Gamma_{ij}=\int_{\Omega}b_{i}(x)b_{j}(x)dx$$ $$\Gamma_{11} =  = \int_{-2}^{4}1 dx = 6$$ $$\Gamma_{12} = \Gamma_{21}= = \int_{-2}^{4}e^x dx =54.4628 $$ $$\Gamma_{13} = \Gamma_{31}= = \int_{-2}^{4}e^{2x} dx =1490.47 $$ $$\Gamma_{22} =  = \int_{-2}^{4} e^{2x} dx = 1490.47$$ $$\Gamma_{23} = \Gamma_{32}= = \int_{-2}^{4}e^{3x} dx =54251.6 $$ $$\Gamma_{33} =  = \int_{-2}^{4}e^{4x} dx =2.2215*10^6 $$

$$\Gamma = \begin{bmatrix} 6 &54.4628&1490.47 \\ 54.4628&1490.47&54251.6  \\ 1490.47&54251.6&2.2215*10^6  \\ \end{bmatrix}$$

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

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

\Gamma = \begin{bmatrix} 6 &54.4628&1490.47 \\ 54.4628&1490.47&54251.6  \\ 1490.47&54251.6&2.2215*10^6  \\ \end{bmatrix}

$$ (Eq )<6>
 *
 * }

The determinant of the matrix is given by:-

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

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

Det \begin{bmatrix} 6 &54.4628&1490.47 \\ 54.4628&1490.47&54251.6  \\ 1490.47&54251.6&2.2215*10^6  \\ \end{bmatrix} =1.13603*10^9

$$ (Eq )<7>
 *
 * }

Hence we can say that the function $$\displaystyle e^{jx} $$ is independent on $$\displaystyle \Omega $$

Refereneces
The determinant is solved using Wolframalpha.com

=Problem 4.6 Derive the strong form and weak form of a elastic bar and spring system=

Statement
A first Course in Finite Elements, Fish & Belychtschko, John Wiley & Sons Ltd, 2007, p.74, problem 3.9

Given
Consider an elastic bar with a variable distributed spring $$p(x) $$along its length as shown in the figure below. The distributed spring imposes an axial force on the bar in proportion to the displacement. Consider a bar of length $$l$$, cross sectional area $$A(x)$$, Young’s Modulus $$E(x) $$ with body force $$b(x)$$ and boundary conditions as shown in the figure below.

A.) Construct the strong form
From the free body diagram, we can obtain the equilibrium equation for the system.
 * {| style="width:100%" border="0"

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

n(x)A(x)\sigma (x)+n(x+dx)A(x+dx)\sigma (x+dx)+[b(x)-p\left( x \right)]dx=0

$$ Where $$n(x)$$ and $$n(x+dx)$$ are the normal of the two cross sections, $$\sigma (x)$$ is the normal stress on that cross section.
 * }

Using Taylor series to expand above equilibrium equation and neglect the high order term, we have
 * {| style="width:100%" border="0"

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

\frac{d}{dx}[A(x)\sigma (x)]+b(x)-p(x)=0

$$ Since $$\sigma (x)=E(x)\frac{du}{dx}$$, plugging this into above equation, we have the strong form as
 * }

 * {| style="width:95%"

 * $$\displaystyle

\frac{d}{dx}[A(x)E(x)\frac{du}{dx}]+b(x)-p(x)=0 $$

The two boundary conditions are:
 * }.

Natural Boundary Condition: $${{\left. E(x)\frac{du}{dx} \right|}_{x=l}}=\bar{t}$$ ;

Essential Boundary Condition: $$u(x=0)=0$$.

B.) Construct the weak form
We now multiply the strong form by the weight function and integrate over the domains over which they hold. This gives us:
 * {| style="width:100%" border="0"

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

\int\limits_{0}^{l}{w(x)\{\frac{d}{dx}[A(x)E(x)\frac{du}{dx}]+b(x)-p(x)\}dx}=0

$$ Simplifying above equation using integrate by part yields the general weak form:
 * }.

 * {| style="width:95%"

 * $$\displaystyle

\int\limits_{0}^{l}{\frac{dw}{dx}A(x)E(x)\frac{du}{dx}dx}=\int\limits_{0}^{l}{w(x)[b(x)-p(x)]}+w(l)A(l)\bar{t} $$

with $$\forall w(0)=0$$
 * },

=Problem 4.7 Installing Calculix and reproducing basic examples to show its usage.=

Statement
1. Installing CGX ( Calculix Graphics Module)

2. Read the manual; signup with usergroup to ask questions if any. Also access archive.

3. Reproduce the basic examples. : disc, cylinder, sphere, sphere- volume, airfoil.

4. Write report for "dummies".Explain to novices how to run cgx ( all cgx commands in basic examples, screenshots,..)

Steps for Installation
a. For installing Calculix, first we need to download the zip file of Calculix from www.bconverged.com

b.Unzip the downloaded Calculix file in a folder of your choice.

c. Run the executable(.exe) file to install the program in your computer.

d. Choose the directory of the folder from where the program would be executed.

e. Finish the installation process.

f. Then we can run the Calculix command by going to Start > Programs > bConverged > Calculix > Calculix Command.

Documentation
a. For reading the manual, go to the Resources page of the website.

b.Click on the Getting Started Guide in the Documentation section.

c. The Getting Started Guide page is opened where some basic commands are explained.

d. To sign up for the Yahoo User's Group, click on the Yahoo User's Group on the Relevant Links on the Resources page of the website.

e. The Yahoo User's Group webpage opens where you can click on "Join this Group" to become a member of the group and ask questions if you have any doubts.

Sample Examples
The basic examples given in the CalculiX CGX were reproduced and the following meshes were obtained. 1. Disc: 2. Cylinder:

3. Sphere:

4. Sphere Volume:

5. Airfoil:

Report
Let us take the disc example and proceed step by step to create a meshed disc.

1. Go to Start - Programs - SciTE, to open the SciTE editor, where the geometry of the example can be modeled using various commands as described below:

2. Once the programming is done, go to File- SaveAs and then save the file with ".fbd" extension.

3. Go to Tools - Pre-processing in the SciTE editor tool bar, to execute the program.

4. Once the program is executed, the CalculiX Graphics window opens up with the model. Now left click on the bottom of the graphics window and proceed as per the following steps.

a. Generating a disc

The CGX commands used for the generation of a disc are- PNT, LINE, GSUR, ELTY, MESH and plus.

The explanation of the commands are given as follows-

i> PNT- This command plots a point on the graphics window. The syntax to print a point with co-ordinates is

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

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

'pnt' |'!' [ <y> <z>]|

$$
 * }

where x, y, z are the co-ordinates of the point "name".

An example 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" |

PNT \quad P001 \quad   0.70711 \quad      -0.00000 \quad      -0.70711

$$
 * }

Here P001 is the name of the point and 0.70711, 0 and -0.70711 are its co-ordinates.

ii> LINE- This command plots a line on the graphics window by using already defined points. The syntax to print a line is

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

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

'line' <name(char<9)>|'!' <p1> <p2> <cp|seq>

$$
 * }

where p1, p2 the points. " " is defined as the name of line. "cp" is the center-point and "seq" is sequential-set.

An example 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" |

LINE \quad L002 \quad P001 \quad P003 \quad p0 \quad 4

$$
 * }

Here P001 and P003 are the names of the points and L002 is the name of the line.

iii> GSUR- This command plots a surface on the graphics window by using already defined lines. The syntax to print a surface is

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

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

'gsur' <name(char<9)>|'!' '+|-' 'BLEND| ' '+|-' <line|lcmb> '+|-' -> <line|lcmb> .. (3-5 times)

$$
 * }

where are the lines defining the surface. " " is defined as the name of surface. The keyword BLEND shows what type of curves (Coons or a NURBS) used for the interior of the surface. '+' or '-' in front of the lines or lcmbs shows the orientation.

An example 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" |

GSUR \quad A001\quad +\quad BLEND \quad - L003 \quad - L002 \quad- L001 \quad- L004

$$
 * }

Here L001, L002 etc. are the lines defining the surface.

iv> ELTY- This command helps to define a specific element type. The syntax to assign element type is

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

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

'elty' [ ] ['be2'|'be3'|'tr3'|'tr3u'|'tr6'|'qu4'|'qu8'|'he8'|'he8f'||'he8i'|'he8r'|'he20'|'he20r']

$$
 * }

where is the set of entities we assigned element type. Two letters in element name 'be2' shows the shape of element such as be: beam, tr: triangle, qu: quadrangle, he: hexahedra, numbers in element name shows the number of nodes. And some other properties can be defined u: unstructured mesh, r: reduced integration, i: incompatible modes, f: fluid element for ccx.

An example 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" |

ELTY \quad all \quad QU4

$$
 * }

v> MESH- This command helps to create mesh of the model. The syntax to make meshing is

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

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

'mesh' ['fast'] ['block'|'lonly'|'nolength'|'noangle'|'length'|'angle']

$$
 * }

where is the set of entities we create mesh. Other features can be used for other parameters.

An example 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" |

MESH \quad all

$$
 * }

vi> plus- This command is used to display the entities of an additional set on the Calculix Graphics page. The syntax of plus is

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

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

'plus' ['n'|'e'|'f'|'p'|'l'|'s'|'b'|'S'|'L']['a'|'d'|'p'|'q'|'v'] -> 'w'|'k'|'r'|'g'|'b'|'y'|'m'|'i'

$$
 * }

, where Nodes =n, Elements =e, Faces =f, Points =p, Lines =l, Surfaces =s, Bodies =b, Nurb Surfaces =S and Nurb Lines =L and White =w, Black =k, Red =r, Green =g, Blue =b, Yellow =y, Magenta =m

An example 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" |

plus \quad ea\quad all

$$
 * }

This command will display the whole image on the Calculix graphics window.

b. Generating a cylinder

The CGX commands used for the generation of a disc are- PNT, LINE, GSUR, ELTY, MESH and plus.

The commands are the same used to generate a disc and are explained above.-

c. Generating a sphere

The CGX commands used for the generation of a sphere are- PNT, LINE, GSUR, ELTY, MESH.

The commands are the same used to generate a disc and are explained above.-

d. Generating a sphere-volume

The CGX commands used for the generation of a sphere are- PNT, LINE, GSUR, ELTY, MESH and GBOD.

The commands are the same used to generate a disc and are explained above except GBOD.-

i> GBOD- This command is used to define or redefine a body in the most basic way. Each body must have five to seven surfaces to be mesh-able. The syntax to print a point with co-ordinates is

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

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

'gbod' <name(char<9)>|'!' 'NORM' '+|-' '+|-' -> .. ( 5-7 surfaces )

$$
 * }

An example 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" |

GBOD \quad B001 \quad NORM \quad - \quad S001 \quad + \quad S002 \quad - \quad S005 \quad - \quad S004 \quad - \quad S003\quad -\quad  S006

$$
 * }

This command will create a body with name B001. NORM is a necessary placeholder which may be used for functionality in the future. S001, S002 etc. are the surfaces which are meshed. The surfaces are followed with a + or - sign in front which in an indication of the orientation of each surface.

e. Generating a Airfoil

The CGX commands used for the generation of a airfoil are- PNT, LINE, GSUR, SEQA, SETA, LCMB, ELTY, MESH and PLUS.

The commands are the same used to generate a disc and are explained above except SEQA, SETA, LCMB.-

i> SEQA- This command is used to make or edit a predefined set marked as a consecutive. The syntax to create a set is

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

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

'seqa' ['pnt' .. <=>]|['afte'|'befo' .. <=>]|['end' .. <=>] $$
 * }

where is the set of entities. ['pnt' .. <=>] is the point we applied to create a consecutive one. The parameter AFTE will insert additional points after the first specified point. The parameter BEFO will insert additional points before the first specified point and the parameter END will add additional points to a sequential set.

An example 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" |

SEQA\quad Q003\quad PNT\quad P004\quad P005\quad P006\quad P00M\quad P00N

$$
 * }

ii> SETA- This command is used to create or edit a set. The syntax to create a set is

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

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

'seta' '!'|'n'|'e'|'p'|'l'|'c'|'s'|'b'|'S'|'L'|'se' <name ..> |['n'|'e' '-' ]

$$
 * }

where is the set of entities. Name of some entities are Nodes n, Elements e, Faces f, Points p, Lines l, Surfaces s, Bodies b, Nurb Surfaces S, Nurb Lines L and names of other sets se

An example 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" |

SETA\quad POINTS\quad P\quad P1\quad P2 $$
 * }

Here it adds the points p1 and p2 to the set 'points'.

iii> LCMB- This command is used to create or edit a combined line. The syntax to create a combined line is

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

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

'lcmb' <name(char<9)>|'!' ['+|-' '+|-' '+|-' ..(up to 14 lines)]| ['ADD' '+|-' '+|-' '+|-' ..(up to 14 lines)]

$$
 * }

where is the name of the combined line we create or edit.

An example 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" |

LCMB\quad U260\quad ADD\quad - U218\quad - U217

$$
 * }

Here it extends U260 by two additional lines U218 and U217.

=Problem 4.8 Find mass matrix and force matrix=

Given
The strong form is

Given $$A=1$$, $$E=2$$,$$\tilde{m}=3$$. The cubic trial solution is given as $$u^{h}(x)=\alpha _{0}+\alpha _{1}(x-3)+\alpha _{2}(x-3)^{2}+\alpha _{3}(x-3)^{3}$$

Where the basis function is given by $$\left\{b_{i}\right\}=\left\{1,(x-3),(x-3)^{2},(x-3)^{3};i=0,1,2,3\right\}$$

a) Mass Matrix
Mass matrix $$\mathbf{\tilde{M}}$$

b) Force Matrix
Force matrix $$\mathbf{F}$$ for dynamics with $$u^{h}(\beta ,t)=g(t)=sin(2t)$$

a) Mass Matrix
The mas matrix $$\mathbf{\tilde{M}}$$ is given by

where $$\mathbf{M_{EE}}=[M_{00}]$$, $$M_{00}=\tilde{M}_{00}$$ and $$\tilde{M}_{00}=\tilde{m}(b_{0},b_{0})$$

$$\tilde{m}(b_{0},b_{0})=\int_{1}^{3}b_{0}\bar{m}b_{0}dx$$ as given in Lecture 19-1

Substituting for $$b_{0}\quad and\quad \bar{m} $$ we get $$\tilde{m}(b_{0},b_{0})=\int_{1}^{3}3dx$$ hence

Similarly $$\mathbf{M_{EF}}=[M_{0j};j=1,2,3]$$, $$M_{0j}=\tilde{M}_{0j}$$ and $$\tilde{M}_{0j}=\tilde{m}(b_{0},b_{j})$$

Clearly we can see that $$\mathbf{M_{EF}}$$ is row matrix

hence we can write $$\tilde{m}(b_{0},b_{1})=\int_{1}^{3}3(x-3)dx$$ which on integration gives $$\tilde{m}(b_{0},b_{1})=-6$$

similarly $$\tilde{m}(b_{0},b_{2})=\int_{1}^{3}3(x-3)^{2}dx=8$$

and $$\tilde{m}(b_{0},b_{3})=\int_{1}^{3}3(x-3)^{3}dx=-12$$

hence

it is observed that $$\mathbf{M_{EF}}=\mathbf{M_{FE}^{T}}$$ hence $$\mathbf{M_{FE}}$$ is column matrix given by

and $$\mathbf{M_{FF}}=[M_{ij};i,j=1,2,3]$$, $$M_{ij}=\tilde{M}_{ij}$$ and $$\tilde{M}_{ij}=\tilde{m}(b_{i},b_{j})$$ where $$\tilde{m}(b_{1},b_{1})=\int_{1}^{3}(x-3)3(x-3)dx=8$$ $$\tilde{m}(b_{1},b_{2})=\int_{1}^{3}(x-3)3(x-3)^{2}dx=-12$$ $$\tilde{m}(b_{1},b_{3})=\int_{1}^{3}(x-3)3(x-3)^{3}dx=\frac{96}{5}$$ $$\tilde{m}(b_{2},b_{1})=\int_{1}^{3}(x-3)^{2}3(x-3)dx=-12$$ $$\tilde{m}(b_{2},b_{2})=\int_{1}^{3}(x-3)^{2}3(x-3)^{2}dx=\frac{96}{5}$$ $$\tilde{m}(b_{2},b_{3})=\int_{1}^{3}(x-3)^{2}3(x-3)^{3}dx=-32$$ $$\tilde{m}(b_{3},b_{1})=\int_{1}^{3}(x-3)^{3}3(x-3)dx=\frac{96}{5}$$ $$\tilde{m}(b_{3},b_{2})=\int_{1}^{3}(x-3)^{3}3(x-3)^{2}dx=-32$$ $$\tilde{m}(b_{3},b_{3})=\int_{1}^{3}(x-3)^{3}3(x-3)^{3}dx=\frac{384}{7}$$

hence

Substituting Eq 4.8.3,Eq 4.8.4,Eq 4.8.5 and Eq 4.8.6 into Eq 4.8.2 we get

b) Force Matrix
The force matrix is given by

where $$\mathbf{\tilde{F}}=\left \{ \frac{\mathbf{F_{E}}}{\mathbf{F_{F}}} \right \}$$ and $$\mathbf{F_{F}}=\left \{ F_{i},i=1,2,3 \right \}$$, $$ F_{i}=\tilde{F_{i}}$$ and $$\tilde{F_{i}}=\tilde{f}(b_{i})$$  as given in Lecture 19-1

$$\tilde{f}(b_{i})=b_{i}(\alpha)h+\int_{1}^{3}b_{i}fdx$$

at $$\alpha$$ natural boundary condition is specified, in this case is $$\alpha=1$$ and $$h$$ is the value of the of the natural boundary condition given by $$h=0.1$$ and $$f$$ in this problem is $$2x$$

hence $$\tilde{f}(b_{1})=(1-3)0.1+\int_{1}^{3}(x-3)2xdx=\frac{-103}{15}$$ $$\tilde{f}(b_{2})=(1-3)^{2}0.1+\int_{1}^{3}(x-3)^{2}2xdx=\frac{42}{5}$$ $$\tilde{f}(b_{3})=(1-3)^{3}0.1+\int_{1}^{3}(x-3)^{3}2xdx=-12$$

stiffness matrix matrix is given by $$\mathbf{K_{EF}}=[K_{0j};j=1,2,3]$$, $$K_{0j}=\tilde{K}_{0j}$$ and $$\tilde{K}_{0j}=\tilde{k}(b_{0},b_{j})$$ and $$\tilde{k}(b_{0},b_{j})=\int_{1}^{3}{b_{0}}'AE{b_{j}}'dx$$

but $$b_{0}=1$$ hence $${b_{0}}'=0$$ therefore $$\tilde{k}(b_{0},b_{j})=0$$ $$\implies\quad\mathbf{K_{EF}}=\begin{bmatrix} 0 & 0 & 0 \end{bmatrix}$$

since $$\mathbf{K_{EF}}=\mathbf{K_{FE}^{T}}$$

from Eq 4.8.1 $$g=sin(2t)$$ $$\implies \quad {g}''=-4sin(2t)$$

hence substituting Eq.4.8.5, Eq4.8.9 and Eq 4.8.10 into Eq 4.8.8 we get

=Contributing Members & Referenced Lecture=

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

title: User:Eml5526.s11.team3/Homework 4