User:Eml5526.s11.team6/hwk3

Given
General one dimensional Model 1.0/Data set 2 on page 1 of Mtg 12:

where $${a_2} = 2$$,  f = 3

Consider $$ {b_j}\left( x \right) = (x+k)^j $$ and k=1 to avoid $${b_j}'\left( 0 \right)=0 $$

To Find
A) If n=2, What is the D.O.F. of this problem?

B) Find two equations that enforce boundary conditions for $$ {u^h}\left( x \right) = \sum\limits_{j = 0}^n $$

C) Find one more equation to solve for $$ {\mathbf{d}} = {\left\{ \right\}_{3 \times 1}} $$ by projecting the residue $$ {\mathbf{P}}\left(  \right) $$ on a basis function $$ b_k\left( x \right)     $$ with  k = 0,1,2 such that the additional equation is linear independent from the equations found in the solution of the (B)

D) Display 3 equations in matrix form $$ {\mathbf{Kd}} = {\mathbf{F}} $$ and observe symmetric property of $$ {\mathbf{K}} $$

E) Solve for $$ {\mathbf{d}} $$

F) Construct $$ u_n^h\left( x \right) $$ and plot $$ u_n^h\left( x \right) \ vs \ u\left( x \right) $$ ....(where $$ u\left( x \right) $$ is the exact solution)

G) Repeat (A) to (F) for n=4 and n=6

H) Compare $$ u_n^h\left( {x = 0.5} \right)   for n = 2,4,6 $$

If error $$ {e_n}\left( {0.5} \right) = u\left( {0.5} \right) - {u^h}\left( {0.5} \right) $$ ,Plot $$ {e_n}\left( {0.5} \right) \ vs \ n $$

Find the degrees of freedom
Since for this problem,'n' starts from 0;

D.O.F.=n+1

D.O.F.=3

Find the equations that enforce the boundary conditions
From (Eq 3.1.2) we can get $$ \sum\limits_{j = 0}^2 {{d_j}{b_j}\left( 1 \right)} = 0 $$

Then rewrite the equation in matrix form.

where

From (3.1.3) we can get $$ \sum\limits_{j = 0}^2 {{d_j}b_j^{\left( 1 \right)}\left( 0 \right)} =  - 4 $$

Then rewrite the above equation in matrix form:

where

Equations 3.1.4 and 3.1.6 are the two equations that enforce the boundary conditions for $$ {u}^h\left( x \right) $$

Project the residue
Because the residue $$ P = \frac{d}\left( {{a_2}\frac{d}{b_j}\left( x \right)} \right) - f\left( x \right)$$, after projecting the residue we get:

$$ \sum\limits_{j = 0}^2 {\left\{ {\int_0^1 {{b_k}\left( x \right)\left[ {\frac{d}\left( {{a_2}\frac{d}{b_j}\left( x \right)} \right)} \right]dx} } \right\}} {d_j} = - \int_0^1 {{b_k}\left( x \right)f\left( x \right)dx} $$

Then rewrite the above equation in matrix form:

where

$$ {K_{kj}} = 2\int_0^1 {{b_k}\left( x \right)b_j^{\left( 2 \right)}\left( x \right)dx} $$

$$         = 2\int_0^1 (x+1)^k. j.(j-1).(x+1)^{j-2} dx $$

It is clear that matrix K is not symmetric.

Matrix d
In order to simplify the computation, the following fortran code was used based on equations (3.1.8),(3.1.9) and (3.1.10) for n=1-10

Using the above codes we can get solutions for all the possible inputs as follows:

Note: There is a slight difference in final answers when using fortran code and matlab code.

The exact solution of eq.3.1.1 by enforcing the boundary condition is

$$ u(x)= -\frac{3x^2}{2}-8x+\frac{19}{2} $$

Plot the exact and approximate solutions
The results are compared with the exact solution as following:



From the above figure it can be seen that the approximate solution almost equals to exact solution for n=4 and above.

Error function


The above figure shows the value of error function with respect to value of n.

From the figure it is clear that error function follows random path for lower values on n and for higher values error function decreases with increase in n.

Solution can be considered to be converging after n=8.

Compare to result using Matlab codes
Integration can be obtained using following matlab code:

Important notes:

From the Matlab results, the solution is exactly the same as the exact solution because we use second second order polynomial to approximate the original function which is also a second order polynomial. The error will be zero since n=2 and the vector d won't change when increasing n. This phenomenon will occur again in HW 3.7 later since both problems have very similar basis with a second order polynomial. The reason that we have some differences between our Fortran results and Matlab results is we do have some numerical errors when solving linear eqautions using Gauss's method by Fortran codes, since the error is so small and results are also very close to what we got from Matlab codes.

=Problem 3.2= For the spring system shown in the figure below in given data.

Find
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.

Given
Spring system with closed ends as shown in the figure



Solution
A: The elements and nodes are numbered as follows: Element numbers are indicated in parentheses and node numbers are indicated at the bottom of the figure.



B: The elements of the stiffness matrices are:

Where global numbers corresponding to the element nodes are indicated above each column and to the right of each row.

By direct assembly, we get the global stiffness matrix as follows

The displacement and force matrices are as follows

Combined force matrix can be written as

C: Now the global system of equations is

As first two displacements are prescribed, partition after two rows and two columns. We get,

Where,

From eq.(3.2.10), we can write

From equations (3.2.12), (3.2.13), (3.2.14), (3.2.15), we get,

and

Solving above two equations, we get

D:Calculation of the reaction forces

From eq.(3.2.10), we can also write

From equations (3.2.11), (3.2.13), (3.2.14), (3.2.19), we get,

and

Solving above two equations, we get

=Problem 3.3=

Given
Referred from - Pg 37, Reference (1)

Consider the truss shown in the following FIGURE 3.3.1



-Nodes A and B are fixed.

-A force F = 10 N acts in the positive x- direction at node C. Coordinates of the joints are given in meters.

-Young's modulus E = 1011 Pa

-Cross sectional area for all bars, A = '''2. 10-2 m2'''

Find
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 stresses & reactions

Number the Element and Nodes
Node and element numbering is done in the following fashion. Refer to the FIGURE 3.3.2 as given below-



Assemble the Global Stiffness and Force Matrix
We use the Following formula Referred from - Pg 30, Reference (1) for Elemental Stiffness Matrix of each element,

$$\displaystyle (Eq. 3.3.1) $$

Where, $${k^{e}} = \frac{l^{e}}$$ and the superscript e denotes the element number.

Now in this case as seen from the given data, A and E for all the 4 elements have same but distinct values. So we can eliminate the superscript e for A and E of each element.

Element 1

 * Element 1 goes from Node 1 to Node 3


 * Orientation : 90 degrees to the Positive X-direction, $$\phi^{(1)} =90^{o}$$ . Hence, $$cos(90^{o})=0,sin(90^{o})=1$$ and $${{l}^{(1)}}=1$$

Using Eq. 3.3.1 for Element 1,

Element 2

 * Element 2 goes from Node 2 to Node 4


 * Orientation : 90 degrees to the Positive X-direction, $$\phi^{(2)} =90^{o}$$ . Hence, $$cos(90^{o})=0,sin(90^{o})=1$$ and $${{l}^{(2)}}=1$$

Using Eq. 3.3.1 for Element 2,

Element 3

 * Element 3 goes from Node 2 to Node 3


 * Orientation : 135 degrees to the Positive X-direction, $$\phi^{(3)} =135^{o}$$ . Hence, $$cos(135^{o})=\frac{-1}{\sqrt{2}},sin(135^{o})=\frac{1}{\sqrt{2}}$$ and $$l^{(3)}=\sqrt{2}$$

Using Eq. 3.3.1 for Element 3,

Element 4

 * Element 4 goes from Node 3 to Node 4


 * Orientation : 0 degrees to the Positive X-direction,$$\phi^{(4)} =0^{o}$$, $$cos(0^{o})=1,sin(0^{o})=0$$, $$l^{(4)}=1$$,

Using Eq. 3.3.1 for Element 4,

GLOBAL STIFFNESS MATRIX
Using above Elemental Matrices, we construct the Global Stiffness Matrix for the given Truss System as follows,

FORCE MATRIX
The force matrix contains reactions and externally applied force components at various nodes.

Let Rx = Reaction at node in x -direction, Ry = Reaction at node in y -direction, subscripts denote the node number

Also node 3 is free to move in both x- and y- direction and is free of any external forces.

And node 4 is free to move in both x- and y- direction and has an external force of 10 N applied in the x- direction as per the given data in the problem statement.

The Force Matrix is assembled as follows -

NODAL DISPLACEMENT MATRIX
The noddal displacement matrix contains displacement components at various nodes.

Let ux = Displacement at node in x -direction, uy = Displacement at node in y -direction, subscripts denote the node number

Nodes 1 and 2 are fixed so they have all the components of displacement as zeros (0).

Also node 3 and 4 are free to move in both x- and y- direction.

The Nodal Displacement Matrix is assembled as follows -

Partition the system and solve for the Nodal Displacements.
Now, using the Global stiffness matrix (K), Force matrix (F) and the Nodal displacement matrix (d), we can write the Global System as follows,

$$\displaystyle (Eq.3.3.2) $$

As the displacements for first two nodes are prescribed, we partition the systems after four rows and columns.

System Partition
Also we use the following equation from - Pg 23, 'A first course in finite elements' - Fish & Belytschko -2007 and compare our global system with it to get the nodal reactions and elemental stresses using $$\displaystyle (Eq.3.3.2) $$ & $$\displaystyle (Eq.3.3.3)  $$ -

$$\displaystyle (Eq.3.3.3) $$

On comparison then we get the following partitioned matrices -

Nodal Displacements
We solve the following system of equations to get the nodal displacement matrix,

But $$\bar{d}_{E}$$ is a Zero Matrix, which will yield, $$d_{F}=K_{F}^{-1}.f_{F}$$ Solving using the above equation, the unknown nodal displacements are

$$\displaystyle (Eq.3.3.5) $$ -

Nodal reactions
Using the above equation the unknown reaction values are

$$\displaystyle (Eq.3.3.7) $$

Elemental Stresses
The stresses in the elements are given by the following expression from page 34 of the reference (1)

For element 1:


 * $$ cos90^{\circ} = 0,\quad sin90^{\circ}=1, $$



\mathbf{d}^{(3)}=l^{(3)}/AE\begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x}\\ u_{2y} \end{bmatrix}, \qquad \sigma ^{(3)}=1/A\begin{bmatrix} 0 & -1 &0 &1 \end{bmatrix}\begin{bmatrix} 0\\ 0\\ 38.2843\\ 10 \end{bmatrix} = 500N/m^{2}

$$

For element 2:


 * $$ cos135^{\circ} = -0.707,\quad sin135^{\circ}=0.707, $$



\mathbf{d}^{(2)}=l^{(2)}/AE\begin{bmatrix} u_{2x}\\ u_{2y}\\ u_{3x}\\ u_{3y} \end{bmatrix}, \qquad \sigma ^{(2)}=1/A\begin{bmatrix} 0 & -1 &-0 &1 \end{bmatrix}\begin{bmatrix} 0\\ 0\\ 48.2843\\ 0 \end{bmatrix} = 0 N/m^{2}

$$

For element 3: 


 * $$ cos90^{\circ} = 0,\quad sin90^{\circ}=1, $$



\mathbf{d}^{(1)}=l^{(1)}/AE\begin{bmatrix} u_{2x}\\ u_{2y}\\ u_{4x}\\ u_{4y} \end{bmatrix}, \qquad \sigma ^{(1)}=1/A\begin{bmatrix} 0.5 & -0.5 &-0.5 &0.5 \end{bmatrix}\begin{bmatrix} 0\\ 0\\ 38.2843\\ 10 \end{bmatrix} = -1000N/m^{2}

$$

For element 4:


 * $$ cos0^{\circ} = 1,\quad sin0^{\circ}=0, $$



\mathbf{d}^{(4)}=l^{(4)}/AE\begin{bmatrix} u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix}, \qquad \sigma ^{(4)}=1/A\begin{bmatrix} -1& 0 &1 &0 \end{bmatrix}\begin{bmatrix} 38.2843\\ 1\\ 48.2843\\ 0 \end{bmatrix} = 10/(2*10^{-2})= 500N/m^{2}

$$

Plot exact and approximate solutions


As you can see the approximate solutions is same as exact solution

Note:Since problems 3.1 and 3.7 have similar basis, the graph of solution follows similar path. Approximate solution equals exact solution for n=2 and then result remains same for higher values of n.

Plot the error function
As error is 0 throughout, graph will be same as x-axis.

=Problem 3.8=

Given
Show that the weak form of

With Essential and Natural Boundary Conditions given by is given by

Solution
Choose $$ w\left ( x \right ) $$ such that $$ w \left ( 3 \right )=0 $$ Multiply the governing equation and the natural boundary condition by an arbitrary weight function:

Use integration by parts to simplify (Eq.3.8.5)

where therefore, (Eq.3.8.5) simplifies to Given that $$ w\left ( 3 \right )=0 $$, (Eq.3.8.9) reduces to Finally, apply the natural boundary condition (Eq.3.8.2) to reduce (Eq. 3.8.10) to Rearrange the terms in (Eq.3.8.11) to yield If (Eq.3.8.12) is equivalent to (Eq.3.8.4), then (Eq.3.8.1-3) from the Given section are equivalent to (Eq.3.8.4).

=Problem 3.9=

Given
Consider a trial (candidate) solution of the form $$u\left ( x \right )=\alpha _{0}+\alpha _{1}\left ( x-3 \right )$$ and a weight function of the same form. Obtain a solution to the weak form in Problem 3.8 (Eq.3.8.4). Check the equilibrium equation in the strong form in Problem 3.8; is it satisfied? Check the natural boundary condition; is it satisfied?

Solution
Given the candidate solution,

Assume,

The weight function must disappear at x=3, therefore

For the trial solution to be admissable, the function must satisfy the essential Boundary Conditions (Eq.3.8.3),

Substitute (Eq.3.9.4) and (Eq.3.9.5) into (Eq.3.9.1) and (Eq.3.9.2) to get

Also, we will need the derivatives of (Eq.3.9.6) and (Eq.3.9.7) in order to obtain the solution to the weak equation given.

Substitute (Eq.3.9.3,6-9) into the weak form (Eq.3.8.4)

Reduce and isolate $$ \beta _{1} $$

Solution must hold for all $$ \beta _{1} $$, therefore, $$ 2\ast 10^{5}\alpha _{1} + 6.47=0 $$

Substitute (Eq.3.9.12) into (Eq.3.9.6) to obtain the linearly approximated solution, i.e.,

Given

Next, substitute (Eq.3.9.13-14) into the Equilibrium Equation to determine whether or not the Equilibrium Equation is satisfied.

$$2x\neq 0 \; \therefore $$ The Equilibrium Equation is not satisfied.

Finally, we check to see if the Natural Boundary Conditions are satisfied.

$$ \sigma \left ( 1 \right )=\left ( E\frac{du}{dx} \right )_{x=1}=0.1\neq -3.23 \; \therefore $$ The Natural Boundary Conditions are not satisfied.

To summarize,

=Problem 3.10= Refer to lecture slide [[media:Fe1.s11.mtg17.djvu|17-2]] for problem assignment.

Given
Strong form:

Weak form:

Find
Consider a trial (candidate) solution and weight functions of the form

Find a solution to the weak form, then check if the equilibrium equation and natural boundary condition in the strong form are satisfied.

Solution
From $$\displaystyle (Eq. 3.10.2)$$ $$w(3) = 0$$ to satisfy the homogeneous essential boundary condition.

Likewise, to satisfy the essential boundary condition of the strong form, $$u(3) = 0.001$$

Next we will evaluate the derivatives of the trial solution and the weight function.

Plugging $$\displaystyle (Eq. 3.10.7)$$ and $$\displaystyle (Eq. 3.10.8)$$ into $$\displaystyle (Eq. 3.10.2)$$ and assuming A and E are uniform in the x-direction (i.e. not functions of x):

For legibility, each term in $$\displaystyle (Eq. 3.10.9)$$ will be evaluated separately.

Now, substituting equations $$\displaystyle (Eq. 3.10.10)$$,$$\displaystyle (Eq. 3.10.11)$$, and $$\displaystyle (Eq. 3.10.12)$$ back into $$\displaystyle (Eq. 3.10.9)$$:

Rearranging and isolating $$\beta_1$$ and $$\beta_2$$:

Equation $$\displaystyle (Eq. 3.10.14)$$ must hold for any arbitrary weight function $$w(x)$$, so long as it satisfies the homogeneous essential boundary condition. It then also follows that it must also hold for any arbitrary $$\beta_1$$ and $$\beta_2$$, leading to the following algebraic equation:

Pre-multiplying both sides $$\displaystyle (Eq. 3.10.15)$$ by the inverse of the coefficient matrix on the left side yields:

Substituting $$\displaystyle (Eq. 3.10.17)$$ back into $$\displaystyle (Eq. 3.10.3)$$ yields the quadratically approximated solution to the weak form $$\displaystyle (Eq. 3.10.2)$$.

Substituting $$\displaystyle (Eq. 3.10.19)$$ back into the equilibrium equation from the strong form $$\displaystyle (Eq. 3.10.1)$$

Therefore the equilibrium equation is not satisfied over the whole interval ]1,3[.

Now using $$\displaystyle (Eq. 3.10.19)$$ to check the natural boundary condition

Therefore, the solution $$\displaystyle u^{\text{quad}}(x)$$ satisfies the natural boundary condition only in the case where $$ A = \frac{40}{3}$$.

=Problem 3.11=

Given
General one dimensional Model 1.0/Data set 2 on page 1 of Mtg 12:

where $${a_2} = 2$$,  f = 3

Consider $$ {b_j}\left( x \right) = 1+{a_j}cos(jx)+{b_j}sin(jx) $$

To find
A) If n=2, What is the D.O.F. of this problem?

B) Find two equations that enforce boundary conditions for $$ {u^h}\left( x \right) = \sum\limits_{j = 0}^n $$

C) Find one more equation to solve for $$ {\mathbf{d}} = {\left\{ \right\}_{3 \times 1}} $$ by projecting the residue $$ {\mathbf{P}}\left(  \right) $$ on a basis function $$ b_k\left( x \right)     $$ with  k = 0,1,2 such that the additional equation is linear independent from the equations found in the solution of the (B)

D) Display 3 equations in matrix form $$ {\mathbf{Kd}} = {\mathbf{F}} $$ and observe symmetric property of $$ {\mathbf{K}} $$

E) Solve for $$ {\mathbf{d}} $$

F) Construct $$ u_n^h\left( x \right) $$ and plot $$ u_n^h\left( x \right) \ vs \ u\left( x \right) $$ ....(where $$ u\left( x \right) $$ is the exact solution)

G) Repeat (A) to (F) for n=4 and n=6

H) Compare $$ u_n^h\left( {x = 0.5} \right)   for n = 2,4,6 $$

If error $$ {e_n}\left( {0.5} \right) = u\left( {0.5} \right) - {u^h}\left( {0.5} \right) $$ ,Plot $$ {e_n}\left( {0.5} \right) \ vs \ n $$

Find the degrees of freedom
Since for this problem,'n' starts from 0;

D.O.F.=n+1

D.O.F.=3

Find the equations that enforce the boundary condition
From (Eq 2.9.2) we can get $$ \sum\limits_{j = 0}^2 {{d_j}{b_j}\left( 1 \right)} = 0 $$

Then rewrite the above equation in matrix form:

where

From (3.1.3) we can get $$ \sum\limits_{j = 0}^2 {{d_j}b_j^{\left( 1 \right)}\left( 0 \right)} =  - 4 $$

Then rewrite the above equation in matrix form:

where

Equations 3.1.4 and 3.1.6 are the two equations that enforce the boundary conditions for $$ {u}^h\left( x \right) $$

Matrix d
Following matlab code is used to compute integration

Plot exact and approximate solutions


The above figure compares the exact solution and the approximated solution

It can be seen that at n=5, approximate solution almost equals exact solution

Plot error function


The above figure compares error between the exact and approximate solution at x=0.5

=Contributing Team Members= Eml5526.s11.team6.vork 03:54, 16 February 2011 (UTC)

Eml5526.s11.team6.deshpande 05:36, 16 February 2011 (UTC)

Eml5526.s11.team6.joglekar 06:17, 16 February 2011 (UTC)

Eml5526.s11.team6.tupsakhare 08:49, 16 February 2011 (UTC)

Eml5526.s11.team6.gravois 16:48, 16 February 2011 (UTC)

Eml5526.s11.team6.kurth 18:21, 16 February 2011 (UTC)

eml5526.s11.team6.lee 20:31, 16 February 2011 (UTC)