User:Eml5526.s11.t1.ab/HW6

=Problem 2:Weak Form Using Quadratic Lagrangian elemental Basis Function With Uniform Discretization.=

Given
General 1-D Model 1 Data set 1. The Partial Differential equation for this data set is,

Find
1)For $$\displaystyle n_{el}=2$$ Compute $$\underline{\tilde{K}}$$ Using Equation 1 On page 30-5. Display $$\underline{\tilde{K}},e=1,2$$

2)Compute $$K^{e}$$, $$L^{e}$$ for $$e=1,2$$

3)Compute $$\underline{\tilde{K}}^{e}=\underline{\tilde{L}}^{e}\underline{K}^{e}\underline{L}^{e}, for e=1,2$$ & Compare to 1

4)Plot all QLEBF for $$n_{el}=3$$.

5)Plot $$\displaystyle \tilde{u_{n}}^{h}$$ vs $$\displaystyle u$$ and [$$\tilde{u_{n}}^{h}\left ( 0.5\right )-u\left (0.5\right )$$] vs $$\tilde{n}$$

Solution
Lets Obtain the Weak form first of all as follow

Performing Integration by parts technique we get

Solving this we get,

As we do not know $$ \frac{du}{dx}\left ( 1\right )$$ we choose $$w\left ( 1 \right )=0$$ and moreover $$u^{h}\left(1 \right )=4$$. To use quadratic Lagrangian elemental basis function we divide the dominium in elements, on each element we apply the weak form by using quadratic Lagrangian Polynomial. For that we need 3 nodes per element. Each element with global nodes i,i+1,i+2 will have 3 local nodes 1,2,3. Lets call $$x_{1},x_{2},x_{3}$$the co-ordinates at node 1,2,3 respectively.

To use natural co-ordinate $$\epsilon$$ instead of Co-ordinate $$x$$ we use the equation 4 on page 30-3 for that we know the natural co-ordinates of the nodes are -1,0,1 & local nodes 1,2,3 respectively. Using Lagrange Polynomial to map between two Co-ordinate system we have:

substituting the above values in the $$\displaystyle (6.2.7)$$ we get,

In this way we have found a shape function, and the later equation can be written as follow,

In the same way we can interpolate the function $$u$$ & $$w$$, therefore,

The Derivative of function $$u$$ with respect to $$\epsilon$$ is given by, In the same way, If we place the Node 2 at the midpoint of the element we can write the relation between the two co-ordinates as Taking the derivative on both sides we get, The derivative of $$u$$ & $$w$$ with respect to $$x$$ are calculated by using chain rule which is as below, Replacing all the terms in weak form as a function of $$\epsilon$$ we obtain This can be further simplified as

From the later Equation we can get the Stiffness Matrix of one element as follows,

and the force vector of the element will be,

From the later Equation we can get the Stiffness Matrix of one element as follows,

There will be additional following factor,

From the later Equation we can get the Stiffness Matrix of one element as follows,

Now only count for degree of freedom at x=1.

1) For $$\displaystyle n_{el}=2$$ Compute $$\underline{\tilde{K}}$$ Using Equation 1 On page 30-5. Display $$\underline{\tilde{K}},e=1,2$$.
Now according to the Lagrange Interpolation we have following equation,

Also,

Replacing in the weak form above we get, The last term is taken only for the last element. Further the stiffness matrix is

Performing the integration for the cone shown below for one element we obtain, Hence the global stiffness matrix will be,

2) Compute $$K^{e}$$, $$L^{e}$$ for $$e=1,2$$.
$$K^{e}$$ is already calculated above. For $$L^{e}$$ As the element has 3 degrees of freedom & the total number of degrees of freedom is 5 The L matrix will be of size $$3\times5$$. Therefore,

3) Compute $$\underline{\tilde{K}}^{e}=\underline{\tilde{L}}^{e}\underline{K}^{e}\underline{L}^{e}, for e=1,2$$ & Compare to 1.
Now, Solving the above matrix we get,

For element 2 we have,

Solving the above matrix we get,

Hence Combining equations $$\displaystyle (6.2.17)$$ & $$\displaystyle (6.2.18)$$ we get global stiffness matrix as follows,

Eventually Comparing $$\displaystyle (6.2.15)$$ with $$\displaystyle (6.2.19)$$ we found that the Stiffness matrices are identical.

4) Plot all QLEBF for $$n_{el}=3$$.
The Quadratic Lagrange Elemental basis Function is plotted using the following MATLAB code for $$\displaystyle n_{el}=3$$.

Appendix
 Matlab Code: 

The plot is as follow.



5) Plot $$\displaystyle \tilde{u_{n}}^{h}$$ vs $$\displaystyle u$$ and [$$\tilde{u_{n}}^{h}\left ( 0.5\right )-u\left (0.5\right )$$] vs $$\tilde{n}$$.
These Two plots are plotted as below Using the mentioned MATLAB code.

Appendix
 Matlab Code: 





===Comment: From the above plot we can observe that the solution converges at nel=16. The Convergence of $$10^{-6}$$ was difficult to achieve therefore the convergence criteria was modified to $$E=\left(u_{0.5}^{i}- u_{0.5}^{i-1}\right )$$.===

=Problem 10: Newtons Law of Colling=

Given
The partial Differential Equation is

1.Develop Weak form similar to (3)-(6) page 34-3.
The weak form will be as below,

Lets calculate the derivative of the following factor, Therefore, Replacing in the integral we obtain, Applying the Gauss Theorem on the first term we obtain, We can recognize that the Heat Flux is Given by, Eventually We have,

Organizing We get, Simplifying we get, For the Case when Conductivity if the Identity Matrix, we have following equations, Also We can Write as,

2) Develop Semidiscrete Equations (ODE's) similar to (3) page 23-3. Give Detailed Description for all the quantities.
Performing the Discretization on the Dominium we have,

Where N is the vector of the Interpolation (Lagrange Polynomial) and $$\displaystyle u_{e},w_{e}$$ are the Nodal Values of u & w Respectively. The Derivative With respect to $$\displaystyle x_{1}$$, $$\displaystyle x_{2}$$ are Replacing on the Semi Discrete Weak form we obtain, From the Last Equation we can Obtain the Conductivity Matrix, the Heat source vector and the Capacitance matrix. The Conductivity Matrix is Given by

Further the Capacitance Matrix is Given by,

And the Heat Source is

3.Solve G2DM2.0/D1 Using 2D LIBF (similar to HW 6.6) till $$10^{-6}$$ accuracy at center. Plot Solution in 3D w/Contour.
To solve G2DM2.0/D1 Following Routine was developed. The following figure shows the result for 16 points.

To solve the problem using LIBF we will devide the x-axis and Y-axis in 2 parts initially. In that way the Approximated solution by LIBF is given by:

We Have developed a MATLAB code to plot this routine. The MATLAB code is as follows

Appendix
 Matlab Code: 

The Plot for the Routine Developed along with the intended 3-D contour is as below.