text:
User:Egm6321.f10.team2.oztekin/Newwwwwwwww

=Problem 7.1: Two Dimensional Heat problem with fixed boundary considering static and transient case =

Given: The Following 2D Data Set
Given the data set provided below.

Domain: Bi-Unit Square
1) The following domain information $$ \Omega = \bar \omega  = \square $$ which is for a bi-unit square

Equation: PDE, K, f, and Initial Condition
2) The partial differential equation of the form $$ \displaystyle \operatorname{div} \left ( K  \cdot \nabla{u} \right ) + f = \rho c\frac $$ 2.a) $$  {K}= {I}=I_n = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix} $$

Find: 2-D Lagrange Interpolation Basis Functions (LIBF), Linear Lagrange Interpolation Basis Functions(LLEBF) and Solve the PDE for $$ m = n = 2,4,6,8 $$
Use 2-dimensional Lagrange Interpolation Basis Functions (LIBF) and Linear Lagrange Interpolation Basis Functions(LLEBF) with $$ \displaystyle m = n = 2,4,6,8 $$ nodes per element to solve the data set above. Ensure the accuracy of the approximate solution $$ \displaystyle u^h(x = 0, y = 0) = 10^{-6} $$ at the center of the element $$ \displaystyle (x,y) = (0,0) $$

1-) Static case (steady state) $$f\left(\boldsymbol{x} \right)=1$$ in $$\Omega =\square $$

1(a) using 2D LIBF

'''1(b). 2D LLEBF'''

1(b).1 uniform mesh

1(b).2 non-uniform mesh

2-) Dynamic Case(Transient)

$$\rho c=3$$

initial condition: $$u\left( x,t=0 \right)=xy$$, $$\forall x\in \Omega $$

2(a). using 2D LIBF.

$$f\left(\boldsymbol{x},t \right)=0$$, $$\forall \boldsymbol{x}\in \Omega $$, $$\forall t>0$$

$$g\left( \boldsymbol{x},t \right)=2$$ on $${{\Gamma }_{g}}=\partial \Omega $$

2(b). With new conditions:

$$f\left( \boldsymbol{x},t \right)=1$$, $$\forall \boldsymbol{x}\in \Omega $$, $$\forall t>0$$

$$g\left( \boldsymbol{x},t \right)=2$$ on $${{\Gamma }_{g}}=\partial \Omega $$

2(b).1 2D LIBF

2(b).2 2D LLEBF

2(b).2.1 uniform mesh

2(b).2.1 non uniform mesh

Solution
The 2-D LIBF can be defined as in (6.5.1)

Where $$ \displaystyle L_{i,m}(x) $$ and  $$ \displaystyle L_{j,n}(y) $$ represent Lagrange Interpolation Basis Functions

For LLEBF we will have only four basis function for each element, which can be constructed from the coordinates of that particular element in the same way as LIBF are constructed.

Initial condition is given for transient case as

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

$$  \displaystyle \begin{matrix} u(x,t=0)=u_{o}(x)=xy & \forall x\in \Omega \end{matrix} $$     (7.1.3) Now in order to find values at nodes at time zero we have to use projection method to get matrices.
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle = $$     (7.1.4) Our goal is
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

Find $$\displaystyle \begin{matrix} d_{j}(0) & st \forall x_{j} \in \Omega \end{matrix}$$

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

$$  \displaystyle \sum_{i} c_{i}\left [ \sum_{j}- \right ]=0 $$     (7.1.5)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle \begin{matrix} \forall c_{i} & \forall x_{i}\in \Omega =c_{F} \end{matrix} $$     (7.1.6)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle G_{i}= $$     (7.1.7)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle c_{F}[\Gamma _{FF}d_{F}-G]=0 $$     (7.1.8)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle d(0)=d_{F}(0)=\Gamma _{FF}^{-1}G $$     (7.1.9)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \Gamma _{FF}=[,i,j\in \eta _{F}] $$     (7.1.10)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

Where $$\displaystyle\eta _{F}$$ is a set of global node numbers not on essential boundary.

DWF form was given in Meeting 23

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

$$  \displaystyle c_{E}\left [ (M_{EE}g^{(s)}+K_{EE}g)+M_{EF}d_{F}^{(s)}+K_{EF}d_{F}-F_{E} \right ]=0 $$     (7.1.11)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle c_{F}\left [ (M_{FE}g^{(s)}+K_{FE}g)+M_{FF}d_{F}^{(s)}+K_{FF}d_{F}-F_{F} \right ]=0 $$     (7.1.12)
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

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

$$  \displaystyle M_{FF}d_{F}^{(s)}+K_{FF}d_{F}=F_{F}-(M_{FE}g^{(s)}+K_{FE}g) $$     (7.1.13)
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

$$  \displaystyle Md^{(s)}+Kd=F $$     (7.1.14) For static case
 * style="width:95%" |
 * style="width:95%" |
 *

 * }

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

$$  \displaystyle \begin{matrix} Kd=F\\ where\\ F=F_{F}-K_{FE}g \end{matrix} $$     (7.1.15) For transient heat transfer case
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 *

 * }

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

$$  \displaystyle \displaystyle \dot{d}=M^{-1}[F-Kd] $$     (7.1.16)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 *

 * }

1(a) using 2D LIBF

Temperature (u) at the center when m=n=5. This value is taken as accurate and error is calculated with reference to this value

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

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

{T_{center}} = {\text{2}}{\text{.294921821297876}} $$     (N)
 *

 * }

1(b) using 2D LLEBF

comment
The plot of $$u-u_h$$ is zig-zag in this case because in this plot both odd number of global node(3*3,5*5,7*7) and even number of global nodes (4*4,6*6) are considered together. If we consider only odd numbers and even numbers separately, curve will smoothly approach to accurate solution (Tcenter = 2.294921821297876).

2(a) using 2D LIBF

first five plot represent the state of temperature with respect to time. Last plot represent system approaching steady state with time. In this case as heat source is zero system should approach to (u=2) at center with time

comment
First five plot represent the state of temperature with respect to time. Last plot represent system approaching steady state with time. In this case as heat source is zero system should approach to (u=2) at center with time.

2(b)1 using 2D LIBF

comment
First three plot represent the state of temperature with respect to time. Last plot represent system approaching steady state with time. In this case as heat source is one system should approach to (u=2.294921821297876) at center with time.

2(b)2 using 2D LLEBF

comment
First three plot represent the state of temperature with respect to time. Last plot represent system approaching steady state with time. In this case as heat source is one system should approach to (u=2.294921821297876) at center with time


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

title: File:Eml5526.s11.team6.HW7.1.1.jpg


title: File:Eml5526.s11.team6.HW7.1.2.jpg


title: File:Eml5526.s11.team6.HW7.1.3.jpg


title: File:Eml5526.s11.team6.HW7.1.5.jpg


title: File:Eml5526.s11.team6.HW7.1.6.jpg


title: File:Eml5526.s11.team6.HW7.1.7.jpg


text: File:FEA 7 2 1.png

title: File:FEA 7 2 1.png


text:
Communication and Identities in Institutional Arenas - Part I/Liinanna final paper

In my final paper on the graduate course KOIIA I, I will deal with some aspects of portraying and negotiating identities in multilingual contexts (through literacy). The theories and ideas presented in the paper that will be published here by May 9, 2011, will form a basis for the "kappa" part of my dissertation.

20 April, 2011

Annaliina Gynne


text: File:X-shaped.png

title: File:X-shaped.png


text:
Talk:Rogers 19 point theory

Development
I hope the formatting I have done makes it a bit more readable. Pehaps we could sort the words and phrases out a bit more. But what would be the best way to do this? above comment added by User:Leutha, 09:16, 21 April 2011. --Abd 19:27, 15 June 2011 (UTC)]]

title: Talk:Rogers 19 point theory


text:
User:Egm5526.s11.team4.zou/NM HW7

=Problem 7.1.3: FE and BE Algorithm Applied to the Integrations of the Logistic Growth Equation=

From the lecture slide Mtg 40

 This problem was solved with referring to S10 homework: S10 Team4's Problem 7.5

Given
The Verhulst Model as stated in the Problem 7.1.2,

Two cases of initial condition,

 * Case 1: initial condition $$ \displaystyle x_0 = 3 < \frac{1}{2}x_{max} $$.
 * Case 2: initial condition $$ \displaystyle x_0 = 9 < \frac{1}{2}x_{max} $$.

Find
Integrate the logistic equation with increasing step size $$\displaystyle h=\hat{h}\cdot {{2}^{k}}$$ by

(1) Foward Euler algorithm

(2) Backward Euler algorithm

Integrate using Backward Euler Algorithm
Do a scaling such that $$\displaystyle \bar{x}(t)=\frac{x(t)}$$, then the Verhulst Model turns into,

Then by Forward Euler Algorithm,

Case 1: $$\displaystyle {{x}_{0}}=3$$
Result using Forward Euler Algorithm when $$\displaystyle k=1$$

Result using Forward Euler Algorithm when $$\displaystyle k=4$$

Result using Forward Euler Algorithm when $$\displaystyle k=7$$

Case 2: $$\displaystyle {{x}_{0}}=9$$
Result using Forward Euler Algorithm when $$\displaystyle k=1$$

Result using Forward Euler Algorithm when $$\displaystyle k=4$$

Result using Forward Euler Algorithm when $$\displaystyle k=7$$

Integrate using Backward Euler Algorithm
From Eq 1.3.2 and Backward Euler Algorithm we have,

Case 1: $$\displaystyle {{x}_{0}}=3$$
Result using Forward Euler Algorithm when $$\displaystyle k=1$$

Result using Forward Euler Algorithm when $$\displaystyle k=7$$

Case 2: $$\displaystyle {{x}_{0}}=9$$
Result using Forward Euler Algorithm when $$\displaystyle k=4$$

Result using Forward Euler Algorithm when $$\displaystyle k=7$$

Comment
As we can see in the Foward Euler Algorithm the result become indented when $$\displaystyle k=7$$ reaches 7. This means the Foward Euler Algorithm's solution may become unstable with step size increasing. While the result by Backward Euler Algorithm still remain stable with $$\displaystyle k=7$$.

title: User:Egm5526.s11.team4.zou/NM HW7


text: File:Feuler.png

title: File:Feuler.png


title: File:Feuler9.png


text: File:Be 3 k=1.png

title: File:Be 3 k=1.png


text: File:Be 3 k=7.png

title: File:Be 3 k=7.png


text: File:Be 9 k=1.png

title: File:Be 9 k=1.png


text: File:Be 9 k=7.png

title: File:Be 9 k=7.png


text: File:Fe 3 k=1.png

title: File:Fe 3 k=1.png


text: File:Fe 3 k=4.png

title: File:Fe 3 k=4.png


text: File:Fe 3 k=7.png

title: File:Fe 3 k=7.png


text: File:Fe 9 k=1.png

title: File:Fe 9 k=1.png


text: File:Fe 9 k=4.png

title: File:Fe 9 k=4.png


text: File:Fe 9 k=7.png

title: File:Fe 9 k=7.png


text: File:Beuler.png

title: File:Beuler.png


title: File:Beuler3-2.png


title: File:Beuler3-3.png


title: File:Beuler9.png


title: File:Beuler9-2.png


title: File:Beuler9-3.png


title: File:Feuler3-2.png


title: File:Feuler3-3.png


title: File:Feuler9-2.png


title: File:Feuler9-3.png


title: File:Egm6341.s11.7.2.2.jpg


title: File:Egm6341.s11.7.2.3.jpg


title: File:Egm6341.s11.7.2.1.png


title: File:Egm6341.s11.7.2.4.jpg


title: File:Egm6341.s11.7.2.5.jpg


text:
User:Lil wayne~enwikiversity

wats gud

title: User:Lil wayne~enwikiversity


text: File:Nm1 s11 team1 HW7 p0.png

title: File:Nm1 s11 team1 HW7 p0.png


title: File:Nm1.s11.team4.HW7.3.fig1.jpg


title: File:Nm1.s11.team4.HW7.3.fig2.jpg


title: File:Nm1.s11.team4.HW7.3.fig3.jpg


title: File:Nm1.s11.team4.HW7.3.fig4.jpg


title: File:Egm6341.s11.7.2.6.jpg


title: File:Egm6341.s11.7.2.7.jpg


title: File:Egm6341.s11.7.2.8.jpg


title: File:Egm6341.s11.7.2.9.jpg


title: File:Egm6341.s11.7.2.10.jpg


text:
User:Eml5526.s11.team3.sahin/Homework 7

=Problem 7.2 Static and Dynamic Finite Element Modelling and Analysis for Vibrating Membrane=

Given: Problem defined on [[media:fe1.s11.mtg40.djvu|Mtg 40 (c)]]
Strong Form

Essential boundary conditions,

Find
1. Find the static solution $$ \quad U_{s}^h $$ to $$ \quad 10^{-6} $$ accuracy at x = 0.

2. Dynamic (Transient) Solution

2a. Free Vibration Analysis

2b. Plot $$ u^h(0,t) \ vs \ t \quad  for \ t \ \epsilon \ [0,2T_1]$$  and produce a movie of the vibrating membrane for symmetric

2c. Plot $$ u^h(0,t) \ vs \ t \quad  for \ t \ \epsilon \ [0,2T_1]$$  and produce a movie of the vibrating membrane for non-symmetric

i. Weak Form
From the PDE the weak form can be written as:

Now from the theory of Integration by parts, we can write

Hence Eq. 2.1 can be now written as

Now we can write $$\displaystyle d\Omega $$ as $$\displaystyle d\Omega ={{\Gamma }_{g}}\cup {{\Gamma }_{h}}$$ , where

$$\displaystyle {{\Gamma }_{g}}$$ = essential boundary condition. $$\displaystyle {{\Gamma }_{h}}$$ = natural boundary condition.

So by applying the Gauss theorem on the first term we obtain,

We select $$\displaystyle w $$ such that $$\displaystyle w=0 $$ on $$\displaystyle {{\Gamma }_{g}}$$

This Equation is of the form:-

where,

ii. Finding approximate solution
Approximated solution $$ u^{h} \left ( x \right ) $$ and $$ w^{h} \left ( x \right ) $$ :

where $$ \left ( c_{i} \right ) $$ and $$ \left ( d_{j} \right ) $$ are constants and $$ \left ( b_{i} \right ) $$ and $$ \left (b_{j}\right ) $$ are the LIBF. We will use these LIBF as a shape function satisfying the CBS.

2DLIBF are given by

$$L_{i,m}$$ is the Lagrange basis function along x and $$L_{i,m}$$ is the Lagrange basis function along y given by

The LIBF be defined by a function

where,

Capacitance matrix

Conductivity matrix

Force matrix

1. Static solution
We have strong form of vibrating membrane equation (Eq. 2.1) where T = 4, $$f\left( x \right){\text{  =  1}}$$,  $$\quad $$   $$\quad \rho = 3 $$,  $$\quad $$ $$\qquad\frac{\partial^2u}{\partial t^2}=0$$

we can also write for 2D

$$\frac{{{\partial }^{2}}u}{\partial x_{i}^{2}}= \frac{\partial^2u}{\partial x_{1}^2} + \frac{\partial^2u}{\partial x_{2}^2}$$

If we rearrange Eq 2.1 for static case conditions, we get

title: User:Eml5526.s11.team3.sahin/Homework 7


text:
User:Eml5526.s11.team3.akj/Homework 7

=Problem 7.1 Static and Dynamic Finite Element Modelling and Analysis for Heat Problem=

Given: Problem defined on [[media:fe1.s11.mtg40.djvu|Mtg 39 (a)]]
Strong Form

for the static case

Initial condition for dynamic case

Boundary conditions for Dynamic case

Find
1. Find the static solution $$ \quad U_{s}^h $$

2. Dynamic (Transient) Solution

2b. Plot $$ u^h(0,t) \ vs \ t $$ at center till convergence to steady state

ii. Finding approximate solution
2DLIBF are given by

$$L_{i,m}$$ is the Lagrange basis function along x and $$L_{i,m}$$ is the Lagrange basis function along y given by

The LIBF be defined by a function

In order to solve for the free degree of freedom $${{d}_{F}},_ – ^ – $$, we construct the Conductance matrix, Heat Source  matrix and Capacitance matrix as follow:
 * {| style="width:100%" border="0"

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

\tilde{\kappa}(\omega,u)=\int_{\Omega}\bigtriangledown \omega.\boldsymbol{\kappa }.\bigtriangledown ud\Omega

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

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

\tilde{f}(\omega)=-\int_{\Gamma _{h}}\omega h d\Gamma _{h}+\int_{\Omega}\omega f d\Omega $$
 * }

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

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

\tilde{m}(\omega,u)=\int_{\Omega}\rho c\omega \frac{\partial u}{\partial t} d\Omega $$
 * }

To solve for the static case from Eq. 1.2 since $$\frac{\partial u}{\partial t}\left( {x,t} \right) = 0, f\left( {x,t} \right) = 1 $$ hence
 * {| style="width:100%" border="0"

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

\tilde{m}(\omega,u)=0 $$
 * }

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

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

\tilde{f}(\omega)= \int_{\Omega}\omega d\Omega $$
 * }

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

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

\tilde{\kappa}(\omega,u)=\int_{\Omega}\bigtriangledown \omega.\bigtriangledown ud\Omega

$$
 * }

$$\displaystyle

\tilde{m}(\omega),\tilde{f}(\omega,u),\tilde{\kappa}(\omega,u)$$ calculate all the degrees of freedom of the system.

Approximated solution $$ u^{h} \left ( x \right ) $$ and $$ w^{h} \left ( x \right ) $$ :

where $$ \left ( c_{i} \right ) $$ and $$ \left ( d_{j} \right ) $$ are constants and $$ \left ( b_{i} \right ) $$ and $$ \left (b_{j}\right ) $$ are the LIBF. We will use these LIBF as a shape function satisfying the CBS.

where,

Capacitance matrix

Conductivity matrix

Force matrix

But we are interested in finding out the degree of freedom at the free ends.

So have to solve for the free degrees of freedom. For this, the MATLAB code is generated such that free ends are collected at one end.

The mas matrix $$\mathbf{\tilde{M}}$$ is given by

The $$\displaystyle \tilde{\kappa} $$ matrix is given by

The force $$\displaystyle  \tilde{F} $$ matrix is given by

Since we only solve for the the free degree of freedom, we only use the equation:-

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

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

\mathbf{M_{FF}}\mathbf\dot+ \mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-(\mathbf{M_{FF}}\mathbf\dot{g}+\mathbf{K_{FE}}\mathbf{g})

$$ (Eq 1.17 )
 *
 * }

In the existing problem, $$\displaystyle \mathbf\dot=0 $$ and $$\displaystyle \mathbf{g}$$ = constant. Thus $$\displaystyle \mathbf\dot{g} = 0.$$

Taking this considerations, Eq 1.17 reduces to :-
 * {| style="width:100%" border="0"

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

\mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-\mathbf{K_{FE}}\mathbf{g}

$$ (Eq 1.18 )
 *
 * }

We can reaarange to find $$\displaystyle \mathbf{d_{F}} $$ as:-

Taking this considerations, Eq 1.17 reduces to :-
 * {| style="width:100%" border="0"

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

\mathbf{d_{F}}=\mathbf{K_{FF}^{-1}}(\mathbf{F_{F}}-\mathbf{K_{FE}}\mathbf{g})

$$ (Eq 1.19 )
 *
 * }

The Matlab code to generate the plots and LIBF is:

Solution for dynamic case
In the dynamic case, $$\displaystyle \mathbf\dot \ne 0 $$ and $$\displaystyle \mathbf\dot{g} = 0.$$

So, Eq 1.17 can be re-written as:-

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

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

\mathbf{M_{FF}}\mathbf\dot+ \mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-(\mathbf{K_{FE}}\mathbf{g})$$

$$\displaystyle \mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-(\mathbf{K_{FE}}\mathbf{g}+\mathbf{M_{FF}}\mathbf\dot)$$

$$\displaystyle \mathbf{d_{F}}=\mathbf{K_{FF}^{-1}}(\mathbf{F_{F}}-(\mathbf{K_{FE}}\mathbf{g}+\mathbf{M_{FF}}\mathbf\dot)) $$ (Eq 1.20 )
 *
 * }

the initial conditions obtained is:

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

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

\mathbf{d(o)}=\mathbf{M_{FF}^{-1}}\mathbf{G}

$$ where
 * }

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

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

G={{\langle {{b}_{i}},{{u}_{0}}\rangle }_{M}}=\int\limits_{\Omega }\rho c{{u}_{o}}d\Omega $$

$$\displaystyle \mathbf{u}=\mathbf{xy}

$$ (Eq 1.21)
 *
 * }

Using MATLAB code "ode45" and the Equation 1.12 using the boundary conditions in 1.13, we solve the problem.

The Matlab code to generate the plots and LIBF is:

title: User:Eml5526.s11.team3.akj/Homework 7


title: File:Hw7q4.png


title: File:Hw7q4dmatrix.PNG


text:
User:EML5526.S11.Team3.vnarayanan/Homework7

=Problem 7.3 Find static solution to $$10^{-6}\;\;$$ accuracy using Quadrilateral and Triangular elements=

Problem Statement
Let $$\;\;\Omega\;\;$$ be a circular domain. Find static solution such that T = 4 and $$\;\;f(x) = 1\;\;\;\;$$in $$\;\;\;\;\Omega\;\; and \;\;g = 0 \;\;on \;\;\Gamma_{g}= \partial\Omega\;\;$$. Use both quadrilateral and triangular elements. Plot the deformed shape in 3-D.

Solution
The nodes and the elements information were obtained from Calculix and are shown below, for both the cases, quadrilateral elements and triangular elements.

We ran the program for the case of quadrilateral elements. We were unable to get the code to run till the end as it took a lot of time. Hence we terminated the code halfway while it was running and we got the values of 'd'(Displacement). The values are listed below.


text:
Survey research and design in psychology/Tutorials/ANOVA and qualitative analysis/Instructor notes

Timing
Approximate timings:
 * 1) ANOVA
 * 2) General steps (including assumptions) - 20 mins
 * 3) Visualising ANOVA - 10 mins
 * 4) Descriptive statistics - 10 mins
 * 5) One-way repeated measures - 20 mins
 * 6) Mixed design ANOVA - 20 mins
 * 7) ANCOVA - 15 mins


title: File:Hs3.png


title: File:Hs9.png


text:
Talk:Linux/linux-0.01.tar

I could not understand about the kernel. Iam a learner in programming. please explain me about boot.s file.

````veleducationo2````

title: Talk:Linux/linux-0.01.tar


text:
User:Charles Jeffrey Danoff/Mr. Danoff's Teaching Laboratory

This page lists Open Educational Resources (OER) that I, Mr. Danoff of Mr. Danoff's Teaching Laboratory, have developed or am developing. Wikiversity specific resources I have developed or am developing are on my main user page. This page is designed to let others see my involvement with the OER community, for me to keep track myself, as well as a disclaimer for this community so people know my affiliation(s). I am the Founder & CEO of the laboratory. It is a for-profit educational business that I am managing. Along for-fee services, OER are also developed in the lab for free. Unless otherwise noted, all OER listed below are copyleft licensed (e.g. GNU-FDL, CCASA or CeCILL).

How to Build an NBA Dynasty

 * Notes

Paragogy

 * Paragogy Bibliography

Roto 5x5 Snake Drafting 101

 * Month long course starting in P2PU Cycle 6


text: File:Deformed Shape in 3D.png

title: File:Deformed Shape in 3D.png


text: File:Displacement contour.png

title: File:Displacement contour.png


title: File:Feahw7q1aLLEBFplot7.png


text: File:Feahw7q1aLLEBFu h.png

title: File:Feahw7q1aLLEBFu h.png


title: File:Feahw7q1aLIBF.png


text: File:Assignment 7 4 plot.jpg

title: File:Assignment 7 4 plot.jpg


text: File:Feahwq1LIBFu h.png

title: File:Feahwq1LIBFu h.png


text:
User:Eml5526.s11.team6.deshpande/hwk7

=Problem 7.3 Static Solution for Unit Circle with quadratic and triangular elements=

Given
Let $$ \displaystyle \Omega = $$ circle with unit radius.

$$ \displaystyle T = 4 $$

$$ \displaystyle f (x) = 1 \quad in \quad \Omega $$

$$ \displaystyle g = 0 \quad on \quad \Gamma_g = \partial \Omega \quad to \quad 10^{-6} $$ accuracy at center.

Find
For both quad and triangle elements.

sym $$ \displaystyle \Rightarrow $$ use $$\displaystyle \frac{1}{4} $$ of meshes.

Plot deformed shape in 3-D.

Solution
By using MATLAB, quarter circle is created with 9 elements and the datapoints are also located and calculated using MATLAB. Following figure shows the mesing of the unit circle.

The co-ordinates are obtained using matlab code and the nodes are numbered from 1 to 16 for 9 elements.

The boundary conditions are defined for harizontal and vertical axes above as

Hence, we can say that nodes on the arc of the circle are fixed and don't have any displacement and remaining nodes can be moved. The MATLAB code is as follows:

The result of the MATLAB code in the form of Deformed 3D shape is shown below.

As we can see, the boundary of the circle is not deflected but other nodes show deflection.

The displacement contour is as shown below

title: User:Eml5526.s11.team6.deshpande/hwk7


text:
User:Eml5526.s11.team6/hwk7

=Problem 7.1 Solve G2DM1.0/D1 (Static and Transient) using 2D LIBF and LLEBF=

Given
The G2Dm1.0/D1 PDE

where: $$\left[ \right]$$ is an identity matrix The essential boundary condition is $$g(x,t) = 2\ \ \ \ on\ \ \ \ \partial \Omega $$.

Find
1a) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for steady state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

1b1) Solve G2DM1.0/D1 using 2D LLEBF with uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for steady state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

1b2) Solve G2DM1.0/D1 using 2D LLEBF with non-uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for steady state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

2a) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 0\ \ \ \ in\ \ \ \ \Omega = \square $$.

2b1) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

2b2a) Solve G2DM1.0/D1 using 2D LLEBF with uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

2b2a) Solve G2DM1.0/D1 using 2D LLEBF with non-uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

Solution
For Eq.6.5.1, we can write the weak form

1a) 2D LIBF
For 2D LIBF, we choose Lagrange Interpolation Basis Function for each element. For example, we choose n=m=2

 * $$N_{I}(x,y)=L_{i,m}(x)\cdot L_{j,n}(y),\;\; where \;\;I=i+(j-1)m$$

After substituting Eq.7.1.3 into Eq.7.1.2 in which Lagrange Interpolation Basis Function are both trial function and weight function, we integrate over the whole domain to obtain the capacitance, conductance, and flux matrices. These matrices are used to achieve final results for the temperatures at each node using the following Matlab codes.

1b) 2D LLEBF
For the 2D LLEBF case, we choose Linear Lagrange elemental basis functions for each element, i.e.,

After substituting Eq.7.1.4 into Eq.7.1.2 in which Linear Lagrange elemental basis function are both trial function and weight function, we integrate over the whole domain to obtain the capacitance, conductance, and flux matrices. These matrices are used to achieve final results for the temperatures at each node using the following Matlab codes.

1b1) Uniform Mesh
The first solution is calculated using a uniform mesh (a mesh where the nodes are equidistant from each other).

One sign that the solution was correctly coded is that the essential boundary conditions are met in these figures. Notice that along the edges of the temperature distribution graph the temperature is equal to 2. Additionally, because there is a positive heat source within the boundaries of our model we would expect the temperature to rise from the boundaries of the model, which is also apparent. These are good indicators that the solution provided is correct, or at least we are heading in the correct direction.

1b2) Non-Uniform Mesh
The second solution is calculated using a non-uniform mesh (an irregularly shaped mesh that displays the abilities of the robust code adopted from the textbook by Fish and Belytschko).

Similarly, the same boundary conditions are evident, and the rise in temperature in the center of the object is obvious.

2b2b) Non-Uniform Mesh
=Problem 7.2: Using 2D LIBF to Solve For Transverse Displacement of an Elastic Membrane=

Find
Solve the given problem using 2D LIBF until the relative error when evaluated at the point $$(x,y)=(0,0)$$ converges to O(10-6).

For the dynamic cases prepare a movie to show $$u^h(x,t)$$

Solution
2D LIBF take the following form:

where $$L_{i,m}(x)$$ and $$L_{j,n}(y)$$ are 1D LIBF with m and n nodes in the x and y directions, respectively.

Static Case
For the static case the strong form reduces to

Introducing a trial solution and weight function as follows

the discrete weak form becomes

Comments
The code still had unresolved errors at the time of the hw submission deadline, as such the desired convergence (or any for that matter) was not achieved. However the essential boundary condition was satisfied. This code is very computationally expensive as it depends on symbolic integration and therefore higher order basis functions were not able to be calculated due to the code timing out.

=Problem 7.3 Static Solution for Unit Circle with quadratic and triangular elements=

Given
Let $$ \displaystyle \Omega = $$ circle with unit radius.

$$ \displaystyle T = 4 $$

$$ \displaystyle f (x) = 1 \quad in \quad \Omega $$

$$ \displaystyle g = 0 \quad on \quad \Gamma_g = \partial \Omega \quad to \quad 10^{-6} $$ accuracy at center.

Find
For both quad and triangle elements.

sym $$ \displaystyle \Rightarrow $$ use $$\displaystyle \frac{1}{4} $$ of meshes.

Plot deformed shape in 3-D.

Solution
By using MATLAB, quarter circle is created with 9 elements and the datapoints are also located and calculated using MATLAB. Following figure shows the mesing of the unit circle.

The co-ordinates are obtained using matlab code and the nodes are numbered from 1 to 16 for 9 elements.

The boundary conditions are defined for harizontal and vertical axes above as

Hence, we can say that nodes on the arc of the circle are fixed and don't have any displacement and remaining nodes can be moved. The MATLAB code is as follows:

The result of the MATLAB code in the form of Deformed 3D shape is shown below.

As we can see, the boundary of the circle is not deflected but other nodes show deflection.

The displacement contour is as shown below

=Problem 7.4-Solution to a general 2-D problem (Heat Transfer)with Lagrange Interpolation Basis Functions(LIBF) =

Given
Problem statement initiated under problem 6.10 (Re-Done here) in mtgs - [[media:Fe1.s11.mtg37.djvu|37-2]] and [[media:Fe1.s11.mtg37.djvu|38-2]]

2D space with boundary convection
Consider a heat conduction problem in 2-dimensional space with boundary convection as shown in the below Figure 6.10.1.

Boundary Conditions
Above 2D space has temperature $$ T(x,y)$$ as space unknown to solve for. We can state it as a boundary value problem with following boundary conditions specified at 3 different kind of boundaries as : $${\Gamma _g}$$, $${\Gamma _h}$$ and $${\Gamma _H}$$

Essential BC on $${\Gamma _g}$$
$$\displaystyle (Eq. 7.4.1) $$

Note : $$T (x,y)$$ can also can be defined as a varying quantity along $${\Gamma _g}$$

Natural BC on $${\Gamma _h}$$
Refer to the Figure 7.4.2 below. We specify the heat flux on boundary $${\Gamma _h}$$ as $${q_n} = {q_0}$$.

This is as we have : $$ {q_n} = \overline q .\overline n $$ with $$\overline n $$ as the unit normal vector on $${\Gamma _h}$$ Hence we can also write : $${q_n} = ( - k.\overline \nabla T).\overline n  =  - k(\overline \nabla  T.\overline n )$$ OR we then have :

$$\displaystyle (Eq. 6.10.2) $$ Note- $$ \nabla T$$ is positive as explained in the - Assumption- subsection below

We are to consider only the normal component of $$\overline q $$, i.e. $$ {q_n}$$, for heat transactions as that is the only component causing heat flow in or out of the element. The tangential component of $$\overline q $$, i.e. $$ {q_t}$$, is just circulating over the boundary.

Assumption
$$\overline n  $$ is unit vector pointing in the outward direction. This implies that $$ \nabla T$$ is positive and temperature increases as we go out away from the surface towards inside of the body. Also we note that $$\overline n  $$ is a vector quantity normal to the $$ T(x,y)$$ contours in 2D space. As $$ T(x,y)$$ increases from inside to outside, heat flux $${q_n}$$ goes into the volume from outside of the boundary $${\Gamma _h}$$ This implies that when $${q_0}$$ is positive, heat comes into the system

Convection BC on $${\Gamma _H}$$
We have as the convection boundary where heat transfer occurs via heat convection Hence we can say that the normal component $${q_n}$$ of the heat flux $$\overline q $$ equals the heat transfer due to heat convection at $${\Gamma _H}$$. Hence we write-

$${q_n} = h(T - Ta)$$ where, we have - $$T$$ : Temperature at the surface of $${\Gamma _H}$$ $${T_a} = {T_\infty } =$$ Say, Ambient temperature $$h$$ :Coefficient of heat convection at $${\Gamma _H}$$ Thus :

$$\displaystyle (Eq. 7.4.3) $$

Given domain with 2-dimensional model data set 2
The following data set needs to be solved using 2D LIBFs -

Solution
Referenced from Section 6.3 (pg 142 )and Section 8.1 (pg 189 ) of text book - ' A first course in finite element method' by authors Fish & Belytschco

Deriving the Equilibrium Equation for 2D- heat transfer
We first derive the equilibrium equation using elemental plane surface as shown in Figure- 7.4.3 below. The elemental surface has length as 'dx' and height as 'dy'. The inflow and outflow of heat is shown symbolically with arrows at various sections.

Note- Heat flux ' $$ \underline q $$ ' is a vector quantity in real 2D space defined per unit surface area and also $$f = $$ Heat generation per unit volume for the given elemental surface We assume unit dimension along the direction perpendicular to the screen or paper

Now the equilibrium of heat can be written in following manner -

Now we use | Taylor Series Expansion technique and can then write -

$$ {({q_x})_{x + dx}} = {({q_x})_x} + {\left( {\frac} \right)_x}dx + higher\_order\_terms$$ OR we can say with the higher order terms neglected :

Similarly we can write for y- direction -

And now $$ Eq.7.4.4 $$ can be modified as below - Canceling out $$dxdy $$ throughout as we have used it for elemental surface and it has to be true over the entire 2D space,

$$ Eq. 7.4.8 $$ can be written with different notation in the form of Divergence of $$\overline q $$ as - Note that $$ Eq. 7.4.9 $$ is valid for 2D as well as 3D heat transfer with $$\overline q $$ as 2D or 3D vector respectively in real 2D or 3D space. Here we are solving for temperature $$T(x,y)$$ so we use 2D space with $$T(x,y)$$

Constitutive Equation for 2D- heat transfer
Here we use the Fourier Law for 2D space which is stated as : $$\overline q  =  - \left[ k \right]\overline \nabla  T $$ It can be written as - Note- Usually, $$=0 $$ and $$=0 $$ for most of the materials as we do not see flow of heat perpendicular to the direction of the temperature gradient line. Also, we have : $$= $$ for isotropic materials and $$ \ne  $$ for anisotropic materials We consider isotropic materials here hence we can assume : $$= $$. Hence $$ Eq. 7.4.10 $$ now becomes -

Hence we write then :

Governing Equation for 2D- heat transfer
Using $$ Eq.7.4.9 $$ and $$ Eq. 7.4.10 $$ we can say that :

Thus the strong form for 2D - Heat Transfer can be written as follows-

$$\displaystyle (Eq. 7.4.11) $$

Deriving the Weak Form for 2D- heat transfer with boundary convection
Let us select the weighting function $$ w(x,y) $$ such that $$ w(x,y)=0$$ at boundary $$ {{\Gamma }_{g}} $$ where the essential boundary condition is applied. This is actually obvious to say as the essential boundary condition is specified at $$ {{\Gamma }_{g}} $$, the value of $$ T (x,y)$$ will be constant there and will no longer require any weighting function multiplication to approximate it in the given domain. Also this will eliminate the unknown flux at the boundary $${{\Gamma }_{g}} $$

Integrating over volume (for general 3D - case) after multiplying with the weighting function w(x), in the form of weighted residual, we get -

Note- For 2D- case : Put $$ dz $$ = $$1$$ which will give us the Area Integral Now we use Green's Theorem to perform the integration by parts. Our aim to get rid of the derivative on $$ \overline{\nabla }.({{K}_{ij}}\overline{\nabla }T)$$ and a derivative should appear on $$ w(x) $$ We state the general Green's Theorem for scalar fields $$ \varphi $$, T - is as follows -

Here, $${{\nabla }^{2}} $$ serves as Laplacian or Divergence of the gradient. Also, $$\overline{n} $$ is the unit normal at boundary and is variable at each point on the boundary. Thus ,we have a volume integral converted to a volume plus a surface integral Hence, $$ Eq. 6.10.12 $$ can now be modified as follows -

This is general equation for 3D as well as 2D case and can be converted to specific case depending upon whether the gradient is 2D or 3D Now again we note that in the weak form the natural boundary conditions get pulled in automatically, for example- the surface integral in above equation will be split into integral over the 3 boundaries - $${{\Gamma }_{g}} $$, $${{\Gamma }_{H}} $$ , and$${{\Gamma }_{h}} $$ , as follows :

Note- In the above Equation, the first term on the RHS will vanish as the weighting function $$ w (x) $$ vanishes on $$ {{\Gamma }_{g}}$$ as explained. Also here we flipped the signs for convection boundary condition on boundary on $${{\Gamma }_{g}} $$ Now, we put $$ Eq. 7.4.15 $$ into $$ eq. 6.10.14$$ and keep the terms required to get the stiffness matrix $$ \left[ K \right]$$ along with the terms involving temperature on LHS and all the terms with known values on the RHS.

$$\displaystyle (Eq. 7.4.16) $$ $$ Eq. 7.4.16 $$ is the weak form for the given 2D Heat Transfer domain with convection boundary specified as $$ {{\Gamma }_{H}}$$ This equation is of the following form -

It is the continuous form of the weak equation and we need to convert it into a discrete form using the method of Desensitization in FEM. Also note that - 1. In $$ Eq. 7.4.16 $$, temperature is a continuous function of field parameters (x,y). 2. This equation is good for 2D as well as 3D domains. 3. We can handle anisotropic materials by converting conductivity 'k' into matrix form as explained.

Developing the semi-discrete equations
Now the temperature distribution in 2D space can be written as follows with subscript 'e' denoting to an element i.e. temperature at a point (x,y) in the 2D space is a sum over that of the elemental nodes. Also here $$ N_i $$ denotes shape function and $$ T_i $$ denotes nodal temperature values.

This can be written in matrix form as follows where the nodal values matrix is shown by $$\left\{ {{X}_{e}} \right\} $$

Also applying the same for the weighting function, we can write : This can be written in matrix form as follows -

We note that we have got gradient terms in the weak form as seen in $$Eq. 7.4.16 $$. The $$\overline{\nabla }T $$ can be written as :

Here $$ \left[ B \right]$$ denotes the matrix of derivatives of the shape functions as shown in below $$ Eq. 7.4.22 $$.

Similar equation can be written for the weighting function as follows -

Now the LHS of the weak form in $$ Eq. 7.4.16 $$ can be written as for a single element -

Hence, LHS becomes for $$ Eq. 7.4.16$$ as :

Depending upon various boundary conditions the RHS of $$ Eq. 7.4.16 $$ can be written as $$\left[ {{F}_{e}} \right]{{\left\{ \partial {{X}_{e}} \right\}}^{T}} $$ Hence the discrete form Equation for first term of the stiffness matrix becomes -

where, $$ \left[ {{K}_{e1}} \right]=\int\limits_{V}{{{\left[ B \right]}^{T}}\left[ K \right]}\left[ B \right]dv$$

The second term of the stiffness can be written in discrete form as follows-

Now the mass matrix can be written in discretized form from the first term of the continuous weak form of $$ Eq. 7.4.16 $$ as follows -

Now the first term on RHS can be discretized as -

Now the second term on RHS can be discretized as -

And the last term on RHS of the weak form can be discretized as -

Capacitance matrix

Conductivity matrix

Heat source matrix

Semi- discrete equations form the given data set G2DM2.0/D1
Referring to mtg - [[media:Fe1.s11.mtg37.djvu|41-5]], the discrete form equations for the given data set as explained in the 'Given' section above, can be written as follows -

Note- In all the matrices below, we have overlap at Node 13 of the domain as shown in Figure 7.4.5 below -

Now the conduction contribution from $${{\Gamma }_{H}} $$ can be written as follows -

Solve the weak form
The following Matlab codes are used to solve the weak form given above using LLEBF.

compare results for different numbers of elements
We use uniform elements and the results of different numbers of elements are as following

From the above results, we can find these results are identical to each other, which means that our codes are correct to achieve the same results even with different numbers of elements.

We can find "negative temperature" at the intersection corner of natural boundary condition and mixed boundary condition. Although it's quite unphysical to have "negative temperature", we try to adjust the boundary condition to find out the reason why the "negative temperature" exists.

compare results for different essential boundary conditions
First we compare the results of different essential boundary conditions with g=2,4,10,300

From the above results, we can find that we won't have "negative temperature" if g>4. From the 2D temperature distribution, we can find that whatever the essential boundary is, the distribution is similar with each other but the magnitude of maximum and minimum are different. However, the distance between maximum and minimum temperature is identical to each other, which means that the essential boundary can "shift" the results. Therefore maybe the "negative temperature" exists because the essential boundary is a little low for this case.

compare results for different natural boundary conditions
Then we compare the results of different natural boundary conditions with h=6,3,2,1

From the above results, we can find that the lower outward heat flux is, the minimum temperature is, because with lower outward heat flux, the temperature will dissipate slowly. When h = 1, the "negative temperature" disappears because less and less temperature will flow out of the area.

Why do we have "negative temperature"?
It's true that the "negative temperature" exist without any physical meaning, because in the universe the minimum temperature is zero K where "kelvin" is a unit of measurement for temperature. But if we suppose the unit of the temperature is Celsius scale, which is another unit scale for thermal temperature, we can have negative temperature just like what we have in winter. In other word, the actual governing equation we solved in this problem is a Laplace equation which can not guarantee that all the results on all the nodes are not negative, unless we can give it the "true" boundary condition which are physical in real world.

And also from the above procedures that we compare results from different boundary conditions, we can find that the boundary conditions do dominate the temperature distribution and "negative temperature" disappears with certain boundary conditions.

Therefore, We believe that we can have physical temperature distribution with physical boundary conditions, even sometimes "negative temperature" can exist with using different aspects of measure unit.

=Contributing Members=

Eml5526.s11.team6.vork 17:52, 21 April 2011 (UTC)

Eml5526.s11.team6.deshpande 19:27, 21 April 2011 (UTC)

Eml5526.s11.team6.kurth 20:31, 21 April 2011 (UTC)

Eml5526.s11.team6.gravois 20:49, 21 April 2011 (UTC)

Eml5526.s11.team6.tupsakhare 20:54, 21 April 2011 (UTC)

eml5526.s11.team6.lee 20:58, 21 April 2011 (UTC)

Eml5526.s11.team6.joglekar 23:01, 21 April 2011 (UTC)

title: User:Eml5526.s11.team6/hwk7


text: File:Mesh representation.png

title: File:Mesh representation.png


text: File:Hw7 1 b 1.jpg

title: File:Hw7 1 b 1.jpg


title: File:Fe1.s11.team1.HW7.fig3.3.jpg


text: File:Hw7 1 b 2.jpg

title: File:Hw7 1 b 2.jpg


title: File:Fe1.s11.team6.HW7.3.fig1.JPG


text: File:ThreeD plot.JPG

title: File:ThreeD plot.JPG


text: File:Egm6341.s11.team3.HW7.3 x0=3.png

title: File:Egm6341.s11.team3.HW7.3 x0=3.png


text: File:Egm6341.s11.team3.HW7.3 x0=9.png

title: File:Egm6341.s11.team3.HW7.3 x0=9.png


text: File:Disp contour.JPG

title: File:Disp contour.JPG


text: File:7q1b1LIBFplot 3.png

title: File:7q1b1LIBFplot 3.png


text: File:7q1b1LIBFplot 2.png

title: File:7q1b1LIBFplot 2.png


text: File:7q1b1LIBFplot 1.png

title: File:7q1b1LIBFplot 1.png


text: File:7q1b1LIBFplot u h.png

title: File:7q1b1LIBFplot u h.png


text: File:Fig quad4.png

title: File:Fig quad4.png


text: File:2bLLEBFplot 3.png

title: File:2bLLEBFplot 3.png


text: File:2bLLEBFplot 2.png

title: File:2bLLEBFplot 2.png


text: File:Fig quad1.png

title: File:Fig quad1.png


text: File:2bLLEBFplot 1.png

title: File:2bLLEBFplot 1.png


text: File:Fig quad2.png

title: File:Fig quad2.png


text: File:2bLLEBFplot u h.png

title: File:2bLLEBFplot u h.png


text: File:Fig tri4.png

title: File:Fig tri4.png


text: File:Fig tri1.png

title: File:Fig tri1.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.n2.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.n3.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.n4.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.n5.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.n6.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.n7.png


text: File:7q1bLIBFplot 5.png

title: File:7q1bLIBFplot 5.png


text: File:7q1bLIBFplot 4.png

title: File:7q1bLIBFplot 4.png


text: File:Fig tri2.png

title: File:Fig tri2.png


text: File:7q1bLIBFplot 3.png

title: File:7q1bLIBFplot 3.png


text: File:7q1bLIBFplot 2.png

title: File:7q1bLIBFplot 2.png


text: File:7q1bLIBFplot 1.png

title: File:7q1bLIBFplot 1.png


text: File:7q1bLIBFplot u h.png

title: File:7q1bLIBFplot u h.png


text: File:Eml5526.s11.team4.ccook.hw7.problem4.error log.png


text: File:Eml 5526 s11 hw7.jpg

title: File:Eml 5526 s11 hw7.jpg


text: File:Eml 5526 s11 hw7angle.jpg

title: File:Eml 5526 s11 hw7angle.jpg


text: File:Eml 5526 s11 hw7displacement.jpg

title: File:Eml 5526 s11 hw7displacement.jpg


title: File:Eml.5526.s11.tm7.7.4.2.Slide1.png


title: File:Eml.5526.s11.tm7.hw7.1.Slide1.png


text: File:7q1bLIBFplot u h new.png

title: File:7q1bLIBFplot u h new.png


text:
User:Egm5526.s11.team4.zou/HW7

Alternative Solution by Fangfang Zhu:
First derive the equilibrium equation for an arbitrary choose domain showns as:

As showns in the above figure,the equilibrium of heat can be written as:

Using the Taylor series expansion we can write:

So the equilibrium equation can be rewritten as:

Dividing dxdy both side yields:

Written in compact form:

Consider the constitutive equation for 2D in heat transfer. The 2D Fourier Law for 2D space can be written as: $$\displaystyle \vec{q}=-[k]\nabla T$$

Written in cartisan matrix form:

For isotropic material we have $$\displaystyle {{k}_{xy}}={{k}_{yx}}=0$$ and $$\displaystyle {{k}_{xx}}={{k}_{yy}}=k$$ So we can write:

or

Thus the strong form for this problem in index notation can be written as:

Following the routine of deriving weak form we can easily get the weak form for this problem:

Let

The weak form can be written as:

Descritization process: Set the approximated solution as a linear combination of shape functions $$\displaystyle {{N}_{i}}(x)$$ as:

And the approximated weighting function is alinear combination of shape functions $$\displaystyle {{N}_{i}}(x)$$ as:

Where $$\displaystyle \tilde{K}={{\int\limits_{D}{\left[ B \right]}}^{T}}\cdot \left[ K \right]\cdot \left[ B \right]$$

Similoarily:

Where $$\displaystyle \tilde{M}=\int\limits_{D}{{{N}^{T}}\cdot \rho c\cdot N}dD$$

where $$\displaystyle {{F}_{i}}=\left[ \int\limits_{-hTN_{i}^{T}ds}+\int\limits_{{{q}_{0}}N_{i}^{T}ds}+\int\limits_{h{{T}_{a}}{{N}_{i}}ds+\int\limits_{D} – f{{N}_{i}}}dD \right]$$

Thus we can write the discrete weak form:

The temprature distribution

title: User:Egm5526.s11.team4.zou/HW7


text: File:4 A 1.png

title: File:4 A 1.png


text: File:4 A 2.png

title: File:4 A 2.png


text: File:4 A 3.png

title: File:4 A 3.png


text: File:HW7 3 3nel.png

title: File:HW7 3 3nel.png


text: File:HW7 2 2nel.png

title: File:HW7 2 2nel.png


text: File:HW7 5 5nel.png

title: File:HW7 5 5nel.png


text:
User:Eml5526.s11.team7.jb/as

Unfortunately performing the analysis with more than 5x5 nodes was too computational costly since the routine relied on symbolic integration. A routine using adaptive Simpson quadrature was adapted from 7.1 to compute the solution with more than 5x5 nodes. As seen from the values of displacement at the center for each case, no convergence was found.

Results with 6x6 Nodes

 * $$\displaystyle u(0,0)=1.946157094594595$$

title: User:Eml5526.s11.team7.jb/as


text: File:Fe1 s11 team3 HW7 1 fig1.jpg

title: File:Fe1 s11 team3 HW7 1 fig1.jpg


text: File:Fe1 s11 team3 HW7 1 fig2.jpg

title: File:Fe1 s11 team3 HW7 1 fig2.jpg


text:
User:Egm6341.s11.team3/sub7

= Cauchy Distribution =

Given
Cauchy Random Noise
 * {| style="width:100%" border="0"

$$\underline{F} = \underline{I} + \triangle \underline{A} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 7.0.1)
 * }

$$\underline{I} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix},\quad \underline{A} = \begin{pmatrix} -0.2& 1 \\ -1 & -0.2 \end{pmatrix}, \quad \triangle = 0.02 \quad \underline{x_0}= \dbinom{{{x}_{0}}^{1}}{{{x}_{0}}^{2}} $$

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

$$LSSM: \quad \underline{x}_{k+1} = \underline{F} \quad \underline{x}_{k}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 7.0.2)
 * }

Find
Plot $$ ( \underline{{x}_{j}}, \quad j=0,1,2,...)$$ in single-slit diffraction

Solution

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

$$ LSSMRN: \quad \underline{{x}_{k+1}} = \underline{F} \quad \underline{x_k} + \underline{G} \quad \underline{w}_{k+1}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 7.0.1)
 * }

Result is

Given
An ideal mass-spring-damper system with mass m, spring constant k and viscous damper of damping coefficient c is subject to a force F=u.
 * {| style="width:100%" border="0"

$$ \left\{ {\begin{array}{*{20}{c}} \dot x_1 \\ \dot x_2 \\ \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} 0 & 1 \\  -\frac{k}{m} & -\frac{c}{m} \\ \end{array}} \right] \left\{ {\begin{array}{*{20}{c}} x_1 \\ x_2 \\ \end{array}} \right\} + \left[ {\begin{array}{*{20}{c}} 0 & 0 \\  1 & 0 \\ \end{array}} \right] \left\{ {\begin{array}{*{20}{c}} u \\ 0 \\ \end{array}} \right\} $$   $$ (Eq.7.0.2) \ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

Finding
Let $$k=1,/ m=1/2,/ x_0=[0.8,-0.4]^T$$ For $$u=0.5$$ Cauchy noise, and $$c=\frac{3}{2}c_{cr}$$, plot $$x_k$$

Solution
Egm6341.s11.team3.sungsikkim 03:29, 21 April 2011 (UTC)

= =

Given: Verhulst Equation with initial conditions
Verhulst equation:

$$ \frac{d}{dt}x(t) = rx(1 - \frac{x}{x_{max}}) \ $$

$$ x_{max} = 15 \ $$ and $$ r = 1.4 \ $$

with initial condition of $$ x_0 \ $$

for time $$ t \in [0, 10] \ $$

Case 1:

$$ x_0 = 3 < \frac{1}{2} x_{max} \ $$

Case 2:

$$ x_0 = 9 > \frac{1}{2} x_{max} \ $$

Find: By using Hermite - Simpson rule, integrate given Verhulst equation
From given Prof. Vu Quoc's notes for

$$ z_{i+ 1} = z_i + \frac{h/2}{3} \left ( f_i + 4f_{i+\frac{1}{2}}[g(z_i,z_{i+1})],t_{i+\frac{1}{2}} + f_{i+1} \right ) \ $$

integrate this function

Solution
Main reference pages were meeting 38 and 39 notes

From given Prof. Vu Quoc's notes for

$$ z_{i+ 1} = z_i + \frac{h/2}{3} \left ( f_i + 4f_{i+\frac{1}{2}}[g(z_i,z_{i+1})],t_{i+\frac{1}{2}} + f_{i+1} \right ) \ $$

,where

$$ f_i = f(z_i, f_i) \ $$

$$ f_{i+1} = f(z_{i+1}, t_{i+1}) \ $$

$$ f_{i+\frac{1}{2}} = f[g(z_i,z_{i+1},t_{i+\frac{1}{2} })] = f(z_{i+\frac{1}{2}}, t_{i+\frac{1}{2}}) \ $$

,where $$ t_{i+\frac{1}{2}} = t_i + \frac{h}{2} \ $$ and as it was proven in HW 6.6 as well as given by Prof. Vu Quoc's

$$ z_{i+ \frac{1}{2}} = \frac{1}{2}(z_i + z_{i+1}) + \frac{h}{8}(f_i - f_{i+1}) \ $$

,note that $$ h \ $$ is the time step

Also specifying the equation to be of a form $$ F(z_{i+ 1})= 0 \ $$

$$ F(z_{i+ 1})= -z_{i+ 1} + z_i + \frac{h/2}{3} \left ( f_i + 4f([g(z_i,z_{i+1})],t_{i+\frac{1}{2}}) + f_{i+1} \right ) = 0\ $$

and, finally, as mentioned in

$$ z_{i+1}^{(k+1)} = z_{i+1}^{(k)} - \left [ \frac{d}{dz}F(z_{i+1}^{(k)}) \right ]^{-1} F(z_{i+1}^{(k)}) \ $$

in our case, for a given function, we have first order derivative. Therefore

$$ z_{i+1}^{(1)} = z_{i+1} - \left [ \frac{d}{dz}F(z_{i+1}) \right ]^{-1} F(z_{i+1}) \ $$

Then check the condition for convergence as it is shown in

$$ z_{i+1}^{(1)} - z_{i+1} ?< Tolerance \ $$

If condition is not met, choose $$ z_{i+1}^{(1)} = z_{i+2} \ $$ and iterate again to get $$ z_{i+2}^{(1)} \ $$ until the tolerance is met

Egm6341.s11.team3.rakesh 02:24, 21 April 2011 (UTC)

Egm6341.s11.team3.JA 17:00, 21 April 2011 (UTC)

Solution Part 2:
All of the necessary information is given in the problem statement. Therefore, all that is left to do is perform the calculations using the proposed methods. The same MatLAB code from Solution Part 1 may be used, with the exception of the file main.m. This file will be replaced by mainB.m as shown below, where the only difference between these files is that mainB.m uses the angle of attack governed by Equation (7.2.2.B), and mainB does not calculate J or ode45.

Below are the new resulting state variable plots, followed by the new M-File that produced them:

Egm6341.s11.team3.russo 03:50, 20 April 2011 (UTC)

title: User:Egm6341.s11.team3/sub7


title: File:Eml.5526.s11.tm7.7.1.HW71nel.png


text:
User:Egm6341.s11.team3/sub7a

= Cauchy Distribution =

Given
Cauchy Random Noise
 * {| style="width:100%" border="0"

$$\underline{F} = \underline{I} + \triangle \underline{A} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 7.0.1)
 * }

$$\underline{I} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix},\quad \underline{A} = \begin{pmatrix} -0.2& 1 \\ -1 & -0.2 \end{pmatrix}, \quad \triangle = 0.02 \quad \underline{x_0}= \dbinom{{{x}_{0}}^{1}}{{{x}_{0}}^{2}} $$

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

$$LSSM: \quad \underline{x}_{k+1} = \underline{F} \quad \underline{x}_{k}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 7.0.2)
 * }

Find
Plot $$ ( \underline{{x}_{j}}, \quad j=0,1,2,...)$$ in single-slit diffraction

Solution

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

$$ LSSMRN: \quad \underline{{x}_{k+1}} = \underline{F} \quad \underline{x_k} + \underline{G} \quad \underline{w}_{k+1}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 7.0.1)
 * }

Result is

Given
An ideal mass-spring-damper system with mass m, spring constant k and viscous damper of damping coefficient c is subject to a force F=u.
 * {| style="width:100%" border="0"

$$ \left\{ {\begin{array}{*{20}{c}} \dot x_1 \\ \dot x_2 \\ \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} 0 & 1 \\  -\frac{k}{m} & -\frac{c}{m} \\ \end{array}} \right] \left\{ {\begin{array}{*{20}{c}} x_1 \\ x_2 \\ \end{array}} \right\} + \left[ {\begin{array}{*{20}{c}} 0 & 0 \\  1 & 0 \\ \end{array}} \right] \left\{ {\begin{array}{*{20}{c}} u \\ 0 \\ \end{array}} \right\} $$   $$ (Eq.7.0.2) \ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 * }

Finding
Let $$k=1,/ m=1/2,/ x_0=[0.8,-0.4]^T$$ For $$u=0.5$$ Cauchy noise, and $$c=\frac{3}{2}c_{cr}$$, plot $$x_k$$

Solution
Egm6341.s11.team3.sungsikkim 03:29, 21 April 2011 (UTC)

= =

Given: Verhulst Equation with initial conditions
Verhulst equation:

$$ \frac{d}{dt}x(t) = rx(1 - \frac{x}{x_{max}}) \ $$

$$ x_{max} = 15 \ $$ and $$ r = 1.4 \ $$

with initial condition of $$ x_0 \ $$

for time $$ t \in [0, 10] \ $$

Case 1:

$$ x_0 = 3 < \frac{1}{2} x_{max} \ $$

Case 2:

$$ x_0 = 9 > \frac{1}{2} x_{max} \ $$

Find: By using Hermite - Simpson rule, integrate given Verhulst equation
From given Prof. Vu Quoc's notes for

$$ z_{i+ 1} = z_i + \frac{h/2}{3} \left ( f_i + 4f_{i+\frac{1}{2}}[g(z_i,z_{i+1})],t_{i+\frac{1}{2}} + f_{i+1} \right ) \ $$

integrate this function

Solution
Main reference pages were meeting 38 and 39 notes

From given Prof. Vu Quoc's notes for

$$ z_{i+ 1} = z_i + \frac{h/2}{3} \left ( f_i + 4f_{i+\frac{1}{2}}[g(z_i,z_{i+1})],t_{i+\frac{1}{2}} + f_{i+1} \right ) \ $$

,where

$$ f_i = f(z_i, f_i) \ $$

$$ f_{i+1} = f(z_{i+1}, t_{i+1}) \ $$

$$ f_{i+\frac{1}{2}} = f[g(z_i,z_{i+1},t_{i+\frac{1}{2} })] = f(z_{i+\frac{1}{2}}, t_{i+\frac{1}{2}}) \ $$

,where $$ t_{i+\frac{1}{2}} = t_i + \frac{h}{2} \ $$ and as it was proven in HW 6.6 as well as given by Prof. Vu Quoc's

$$ z_{i+ \frac{1}{2}} = \frac{1}{2}(z_i + z_{i+1}) + \frac{h}{8}(f_i - f_{i+1}) \ $$

,note that $$ h \ $$ is the time step

Also specifying the equation to be of a form $$ F(z_{i+ 1})= 0 \ $$

$$ F(z_{i+ 1})= -z_{i+ 1} + z_i + \frac{h/2}{3} \left ( f_i + 4f([g(z_i,z_{i+1})],t_{i+\frac{1}{2}}) + f_{i+1} \right ) = 0\ $$

and, finally, as mentioned in

$$ z_{i+1}^{(k+1)} = z_{i+1}^{(k)} - \left [ \frac{d}{dz}F(z_{i+1}^{(k)}) \right ]^{-1} F(z_{i+1}^{(k)}) \ $$

in our case, for a given function, we have first order derivative. Therefore

$$ z_{i+1}^{(1)} = z_{i+1} - \left [ \frac{d}{dz}F(z_{i+1}) \right ]^{-1} F(z_{i+1}) \ $$

Then check the condition for convergence as it is shown in

$$ z_{i+1}^{(1)} - z_{i+1} ?< Tolerance \ $$

If condition is not met, choose $$ z_{i+1}^{(1)} = z_{i+2} \ $$ and iterate again to get $$ z_{i+2}^{(1)} \ $$ until the tolerance is met

Egm6341.s11.team3.rakesh 02:24, 21 April 2011 (UTC)

Egm6341.s11.team3.JA 17:00, 21 April 2011 (UTC)

Solution Part 2:
All of the necessary information is given in the problem statement. Therefore, all that is left to do is perform the calculations using the proposed methods. The same MatLAB code from Solution Part 1 may be used, with the exception of the file main.m. This file will be replaced by mainB.m as shown below, where the only difference between these files is that mainB.m uses the angle of attack governed by Equation (7.2.2.B), and mainB does not calculate J or ode45.

Below are the new resulting state variable plots, followed by the new M-File that produced them:

Egm6341.s11.team3.russo 03:50, 20 April 2011 (UTC)

title: User:Egm6341.s11.team3/sub7a


title: File:Fe1.s11.tm7.HW7.fig7.2.3a.png


text: File:Abaqus 7 3.png

title: File:Abaqus 7 3.png


text: File:FEA 7 2 2.png

title: File:FEA 7 2 2.png


text: File:Hw7 1 a.jpg

title: File:Hw7 1 a.jpg


text: File:Hcpy 1.png

title: File:Hcpy 1.png


title: File:Fe1.s11.team6.kurth.HW7.2.jpg


title: File:Fe1.s11.team6.kurth.HW7.2.2.jpg


text: File:Hcpy 2.png

title: File:Hcpy 2.png


text: File:Eml5526.s11.team7.hw7.2a n5sim.jpg

title: File:Eml5526.s11.team7.hw7.2a n5sim.jpg


text: File:Triangle n.png

title: File:Triangle n.png


text: File:Trielements.png

title: File:Trielements.png


text: File:Eml5526.s11.team7.hw7.2a n3sim.jpg

title: File:Eml5526.s11.team7.hw7.2a n3sim.jpg


text:
User:Eml5526.s11.team7/HW7

Problem Description
Following the solution to HW 6.5 found hereUse 2-D Lagrange Interpolation Basis Functions until accuracy of 10-6 at the center is found (x,y)=(0,0).

Given
The boundary is given by a biunit square $$\Omega =\bar{\omega }=\square$$

The governing PDE for the General 2-Dim Model is

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

$$\displaystyle div(\mathbf K\cdot grad(u))+f=\rho c\frac{\partial u}{\partial t}$$ (7.1.1)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

Where

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

$$\displaystyle \mathbf K=\mathbf I,
 * style="width:10%; |
 * style="width:10%; |

f=0,

\frac{\partial u}{\partial t}=0$$ (7.1.2)
 *
 * }

The Essential Boundary Condition is given by

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

$$\displaystyle g=2|_{\partial\Omega}$$ (7.1.3)
 * style="width:10%; |
 * style="width:10%; |
 *
 * }

7.1.a) Find Solution using Lagrange Interpolation Basis Functions

7.1.b) Find Solution using Linear Lagrange Element Basis Functions

7.1.a:LIBF Solution
For the x and y one dimensional cases:

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

$$ \displaystyle
 * style="width:10%; |
 * style="width:10%; |

{L_{i,m}}(x) = \prod\limits_{\begin{array}{ccccccccccccccc} {k = 1} \\ {k \ne i} \end{array}}^m {\frac}

$$     (Eq 7.1.4)
 *
 * }

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

$$ \displaystyle
 * style="width:10%; |
 * style="width:10%; |

{L_{j,n}}(y) = \prod\limits_{\begin{array}{ccccccccccccccc} {l = 1} \\ {l \ne j} \end{array}}^n {\frac}

$$     (Eq 7.1.5)
 *
 * }

Combining these for the two dimension case:

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

$$ \displaystyle
 * style="width:10%; |
 * style="width:10%; |

{N_I}(x,y) = {L_{i,m}}(x) \cdot {L_{j,n}}(y),

$$     (Eq 7.1.6)
 *
 * }

Where I = i + (j - 1)m.

Moving to the Matrix Equations:

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

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

m(w,u) + k(w,u) = f(w,u)

$$     (Eq 7.1.7)
 *
 * }

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

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

m(w,u)= \int\limits_\Omega {w\rho c\frac} d\Omega

$$     (Eq 7.1.8)
 *
 * }

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

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

k(w,u) = \int\limits_\Omega {\nabla wK\nabla u} d\Omega

$$     (Eq 7.1.9)
 *
 * }

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

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

f(w,u) = - \int_{\Gamma h} {whd{\Gamma _h}}  + \int\limits_\Omega  {wf} d\Omega

$$     (Eq 7.1.10)
 *
 * }

For this problem $$\displaystyle \frac{\partial u}{\partial t}=0$$ and $$\displaystyle f=1$$. Moving to the matrix form:

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

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

Kd = F

$$     (Eq 7.1.11)
 *
 * }

Now from Eq 7.1.9 we have

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

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

{K_{ij}} = \int_ {(\nabla N_I)I(\nabla N_J)d{w}}

$$     (Eq 7.1.12) The value of u at x=0 and y=0 for n=m=3 is 2.3125 The value of u at x=0 and y=0 for n=m=4 is 2.2469 The value of u at x=0 and y=0 for n=m=5 is 2.2949
 *
 * }

7.1b:LLEBF
For the x and y one dimensional cases:

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

$$ \displaystyle
 * style="width:10%; |
 * style="width:10%; |

{L_{i,m}^e}(x) = \prod\limits_{\begin{array}{ccccccccccccccc} {k = 1} \\ {k \ne i} \end{array}}^m {\frac}

$$     (Eq 7.1.13)
 *
 * }

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

$$ \displaystyle
 * style="width:10%; |
 * style="width:10%; |

{L_{j,n}^e}(y) = \prod\limits_{\begin{array}{ccccccccccccccc} {k = 1} \\ {k \ne j} \end{array}}^n {\frac}

$$     (Eq 7.1.14)
 *
 * }

Combining these for the two dimension element case:

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

$$ \displaystyle
 * style="width:10%; |
 * style="width:10%; |

{N_I^e}(x,y) = {L_{i,m}^e}(x) \cdot {L_{j,n}^e}(y),

$$     (Eq 7.1.15)
 *
 * }

Where I = i + (j - 1)m.

Moving to the Matrix Equations:

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

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

m^e(w,u) + k^e(w,u) = f^e(w,u)

$$     (Eq 7.1.16)
 *
 * }

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

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

m^e(w,u)= \int\limits_{\Omega^e} {w\rho c\frac} d\Omega^e

$$     (Eq 7.1.17)
 *

 * }

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

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

k^e(w,u) = \int\limits_{\Omega^e} {\nabla wK^e\nabla u} d\Omega^e

$$     (Eq 7.1.18)
 *

 * }

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

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

f^e(w,u) = - \int_{\Gamma h} {whd{\Gamma _h}}  + \int\limits_{\Omega^e}  {wf} d\Omega^e

$$     (Eq 7.1.19)
 *

 * }

For this problem $$\displaystyle \frac{\partial u}{\partial t}=0$$ and $$\displaystyle f=1$$. Moving to the matrix form:

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

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

K^ed^e = F^e

$$     (Eq 7.1.20)
 *

 * }

Here's the Matlab code for the solution to this problem:

Giving the following Plot:

Comments
This figure is obviously incorrect, however, we were unable to run down the errors. Analysis of our scattering matrices showed correct alignment with the global node designation (global nodes enumerated initially on the essential boundary, and then the remaining nodes). Additionally the values do not make any sense either. Displaying the d matrix while running the code showed that the coefficients are being designated correctly for the essential boundary. It is likely the error is either in the multiplication of the coefficients to the basis functions, or the substituting in of values for making the plot. More time would be needed to track down the errors in the code

Part 2a: Dynamic(Transient) case for f=0 by using 2D LIBF
For Transient analysis (heat problem):

$${{\mathbf{M}}_{\mathbf{FF}}}{{{\mathbf{\dot{d}}}}_{\mathbf{F}}}+{{\mathbf{K}}_{\mathbf{FF}}}{{\mathbf{d}}_{\mathbf{F}}}={{\mathbf{F}}_{\mathbf{F}}}-\left( {{\mathbf{M}}_{\mathbf{FE}}}\mathbf{\dot{g}}+{{\mathbf{K}}_{\mathbf{FE}}}\mathbf{g} \right)$$

In this problem, $$\mathbf{g}$$ is a constant, so $$\mathbf{\dot{g}}=0$$.

Thus, $${{{\mathbf{\dot{d}}}}_{\mathbf{F}}}=\mathbf{M}_{\mathbf{FF}}^{\mathbf{-1}}\left( {{\mathbf{F}}_{\mathbf{F}}}-{{\mathbf{K}}_{\mathbf{FE}}}\mathbf{g}-{{\mathbf{K}}_{\mathbf{FF}}}{{\mathbf{d}}_{\mathbf{F}}} \right)$$

The initial condition: $$u(\mathbf{x},t=0)=xy$$

$$\rho c=3$$

Heat source: $$f=0$$

The plot for uh vs.t at center:

Part 2b1: Dynamic(Transient) case for f=1 by using 2D LIBF
Similar to part 2a.

Heat source: $$f=1$$

The plot for uh vs.t at center:

Problem Description
Solve the vibrating membrane problem using 2-D Lagrangian Interpolation Basis Functions (LIBF) with $$\displaystyle n=m=2,4,6,8...$$ by until an accuracy of $$\displaystyle 10^{-6}$$ is reached at the center $$\displaystyle (x,y)=(0,0)$$.

For the dynamic case, Solve the free vibration problem:
 * $$ \mathbf{K \phi}=\lambda \mathbf{ M \phi}$$

For the first three eigenpairs:
 * $$ (\lambda_i,\mathbf{\phi_i}), i=1,2,3.$$

Given
The governing partial differential equation (PDE) is

 * $$\displaystyle T\frac{\partial}{\partial x_{i}} \left[ \frac{\partial u}{\partial x_{j}} \right ] + f = \rho \frac{\partial^2 u}{\partial t^2} $$

 * where
 * The domain $$\displaystyle \Omega = \overline{\omega}=\Box$$ is a biunit square

 * $$\displaystyle f(\mathbf{x}) = 1$$ in $$\displaystyle \Omega$$

 * 1) Static Case:
 * $$\displaystyle \frac{\partial^2 u}{\partial t^2 } = 0 $$
 * 2) Dynamic/Transient Case:
 * $$\displaystyle \rho = 3 $$

The essential boundary condition is defined over the entire boundary:

 * $$\displaystyle g=0$$ on $$\displaystyle \partial u$$

Part 1: Static case
The weak form for the vibrating membrane can be written as follows:

\displaystyle

m(w,u) + k(w,u) = f(w,u)

$$

Examining each term and simplifying with the given conditions yields:

\displaystyle

m(w,u) = \int\limits_\Omega {w\rho \frac} d\Omega =0

$$

\displaystyle

k(w,u) = \int\limits_\Omega {T\nabla w\nabla u} d\Omega

$$

\displaystyle

f(w,u) = - \int_{\Gamma h} {whd{\Gamma _h}}  + \int\limits_\Omega  {wf} d\Omega= \int\limits_\Omega  {wf} d\Omega

$$

\displaystyle

Kd = F

$$ Now the basis functions are developed to plug into the preceding equations. For the 1D case the LIBF Basis function can be written:

\displaystyle

{L_{i,m}}(x) = \prod\limits_{\begin{array}{ccccccccccccccc} {k = 1} \\ {k \ne i} \end{array}}^m {\frac}

$$

It is extended to the the 2D case as:

\displaystyle

{N_I}(x,y) = {L_{i,m}}(x) \cdot {L_{j,n}}(y)

$$

where

\displaystyle

I = i + (j - 1)m

$$ Two examples of basis functions, $$\displaystyle N(x,y)$$ are shown below:

(A)n=m=4, I=1

(B)n=4, I=7

One can see from the plots of the basis functions that they both satisfy the constraint breaking solution (CBS). For example the first plot shows the value to be one at the first node (I=1, (x,y)=(-1,-1)) and zero at all other nodes.

Plugging in the basis fucntions yields:

\displaystyle

{K_{ij}} = \int_ T{(\nabla N_I^e)(\nabla N_J^e)d{w^e}}

$$ Since $$\displaystyle T$$ is a constant it can be pulled out of the integral.

\displaystyle

{K_{ij}} = T\int_ {(\nabla N_I^e)(\nabla N_J^e)d{w^e}}

$$ Considering the essential boundary condition is the outline of the biunit, one special treatment is needed. We reconstruct $$\displaystyle {K_f}$$ and $$\displaystyle {F_f}$$, so that we can solve for matrix $$\displaystyle d$$. Firstly, considering

\displaystyle

{L_{i,m}}({x_j}) = {\delta _{i,m}} = \left\{ {\begin{array}{ccccccccccccccc} {1,forj = i} \\ {0,forj \ne i} \end{array}} \right.

$$ For all $$\displaystyle {d_i}$$ belong to the nodes on the essential boundary condition, we can get the value directly:

\displaystyle

{d_i} = 0,where(i \in ess.b.c)

$$ For the unknown or free coefficients$$\displaystyle {d_F}$$, which are located in the middle area of the biunit and number is $${(n - 2)^2}$$, we need to construct $$\displaystyle {K_F}$$ and $$\displaystyle {F_F}$$ to solve for the free coefficients.

\displaystyle

{K_F}{d_F} = {F_F}

$$ Note: The dimension of $$\displaystyle {K_F}$$ is $$\displaystyle {n^2} \times {(n - 2)^2}$$, the dimension of $$\displaystyle {d_F}$$ is $$\displaystyle {(n - 2)^2} \times 1$$ , the dimension of $$\displaystyle {F_F}$$ is $$\displaystyle {n^2} \times 1$$

The following matlab routine was adapted from Team 7 #6.5 to perform the integration and solve for the unknown displacements.

Results with 3x3 Nodes

 * $$\displaystyle u(0,0)=1.946157094594595$$

Results with 4x4 Nodes

 * $$\displaystyle u(0,0)=2.492872807017544$$

Results with 5x5 Nodes

 * $$\displaystyle u(0,0)=1.976986310175030$$

Comments on Results
One can see that each solution satisfies the essential boundary condition. That is, the displacement is equal to zero around the entire boundary. The maximum displacement occurs at the center which is also an expected result. The point furthest from the essential boundary (the center) in this case should have the largest displacement because the tension and transverse forces are constant across the entire domain.

Unfortunately performing the analysis with more than 5x5 nodes was too computational costly since the routine relied on symbolic integration. Also, as seen from the values of displacement at the center for each case, no convergence was found. An alternative routine using adaptive Simpson quadrature was adapted from 7.1 to compute the solution.

Results with 5x5 Nodes using adaptive Simpson quadrature

 * $$\displaystyle u(0,0)=0.073730557049503$$

Part 2: Dynamic case
The initial and boundary conditions for the dynamic case are defined by

 * $$\displaystyle u(\mathbf{x},t=0) = u_s^h (x) \quad \forall \mathbf{x} \in \Omega$$
 * $$\displaystyle \dot{u}(\mathbf{x},t=0) = 0 \quad \forall \mathbf{x} \in \Omega$$
 * $$\displaystyle T=4$$ (tension per unit length)
 * $$\displaystyle f(\mathbf{x}, t)=0$$ (no applied transverse force)

 * where $$\displaystyle u_s^h (x)$$ refers to the calculated transverse displacement of the static case (above).

Before the transient analysis can be performed, the first three eigenpairs ($$\displaystyle \lambda_i, \ \mathbf{\phi}_i$$) need to be determined from the generalized eigenvalue problem. From this the natural frequencies of vibration can be determined. Matlab may be used in this case by the "eig" operator. However the generalized eigenvalue problem must be of the form $$\displaystyle \mathbf{A} \phi = \lambda \phi$$. In our case $$\displaystyle \mathbf{A} = \mathbf{M}^{-1} \mathbf{K}$$.

Part 3: Dynamic case with asymmetric u0
This analysis is identical to the previous section (Problem 7.2-2) with the exception that the initial displacement is asymmetric and is defined and plotted below.

 * $$\displaystyle u(\mathbf{x},t=0) = (x+1) \left(y+\frac{1}{2} \right ) cos \left(\frac{\pi}{2}x \right ) cos \left(\frac{\pi}{2}y \right )\quad \forall \mathbf{x} \in \Omega$$..

7.3:Solve for a circle with unit radius through quadratic and triangular element meshes
Let $$\Omega$$ = circle with unit radius. Find static solution such that T=4, f(x)=1 in $$\Omega$$, g=0 on boundaries to 10^-6 accuracy at center. Use Both quad and triangular elements. Use symmetry 1/4 of the total meshes. plot deformed shape in 3-D

Solution
Nodal distribution of quad elements Element distribution of quad elements Elements	nodes numbers 1 ->	1	2	3	4 2 ->	4	3	5	6 3 ->	6	5	7	8 4 -> 8	7	9	10 5->	2	11	12	3 6->	3	12	13	5 7->	5	13	14	7 8->	7	14	15	9 9->	11	16	17	12 10->	12	17	18	13 11->	13	18	19	14 12->	14	19	20	15 13->	16	21	22	17 14->	17	22	23	18 15->	18	23	24	19 16->	19	24	25	20 Nodal coordinates X              Y 1	0.00E+00	0.00E+00 2	2.50E-01	0.00E+00 3	2.50E-01	2.50E-01 4	0.00E+00	2.50E-01 5	2.40E-01	4.98E-01 6	0.00E+00	5.00E-01 7	2.22E-01	7.42E-01 8	0.00E+00	7.50E-01 9	1.95E-01	9.81E-01 10	0.00E+00	1.00E+00 11	5.00E-01	0.00E+00 12	4.98E-01	2.40E-01 13	4.77E-01	4.77E-01 14	4.38E-01	7.06E-01 15	3.83E-01	9.24E-01 16	7.50E-01	0.00E+00 17	7.42E-01	2.22E-01 18	7.06E-01	4.38E-01 19	6.43E-01	6.43E-01 20	5.56E-01	8.31E-01 21	1.00E+00	0.00E+00 22	9.81E-01	1.95E-01 23	9.24E-01	3.83E-01 24	8.31E-01	5.56E-01 25	7.07E-01	7.07E-01

The displacements obtained... nodes		displacements 1	>	-0.0034 2	>	-0.0009 3	>	-0.0001 4	>	-0.0009 5	>	-0.0001 6	>	-0.0007 7	>	-0.0001 8	>	-0.0007 9	>	-0.0001 10	>	-0.0003 11	>	-0.0007 12	>	-0.0001 13	>	0 14	>	0 15	>	0 16	>	-0.0006 17	>	-0.0001 18	>	0 19	>	0 20	>	0.0005 21	>	-0.0003 22	>	-0.0001 23	>	0.0001 24	>	0.0007 25	>	0.0019

Elements distribution 1	>	1	2	3	 2	>	3	2	4	 3	>	3	4	5	 4	>	5	4	6	 5	>	5	6	7	 6	>	7	6	8	 7	>	2	9	4	 8	>	4	9	10	 9	>	4	10	6	 10	>	6	10	11	 11	>	6	11	8	 12	>	8	11	12	 13	>	9	13	10	 14	>	10	13	14	 15	>	10	14	11	 16	>	11	14	15	 17	>	11	15	12	 18	>	12	15	16

Nodal coordinates of triangular distribution

X	         Y		X	Y 1	>	0.00E+00	0.00E+00 2	>	3.33E-01	0.00E+00 3	>	0.00E+00	3.33E-01 4	>	3.30E-01	3.30E-01 5	>	0.00E+00	6.67E-01 6	>	3.04E-01	6.53E-01 7	>	0.00E+00	1.00E+00 8	>	2.59E-01	9.66E-01 9	>	6.67E-01	0.00E+00 10	>	6.53E-01	3.04E-01 11	>	5.96E-01	5.96E-01 12	>	5.00E-01	8.66E-01 13	>	1.00E+00	0.00E+00 14	>	9.66E-01	2.59E-01 15	>	8.66E-01	5.00E-01 16	>	7.07E-01	7.07E-01 matlab code used cgx code use pnt p0      0.00000        0.00000        0.00000 pnt p00A    0.00000        0.00000        1.00000 pnt p00G    0.70711        0.00000        0.70711 pnt p00I    1.00000        0.00000        0.00000 line l00A p0 p00I 3 line l00C p00A p0 3 line l00G p00A p00G p0 3 line l00I p00G p00I p0 3 gsur a004 + BlEND + l00A - l00I - l00G + l00C ElTY all tr3 MESH all plot ea all

Result
The above program did not execute correctly. The plot and displacements obtained are not correct!! The nodes and elements were obtained by calculiX CGX

Problem Description
An energy balance on a two dimensional unit square is performed. This square has a prescribed temperature of 2 for all y at x=1. The square also has a temperature of 2 for all x with y=-1. This is given by the following equations and is seen in the figure below with $$\displaystyle \Gamma_g$$.

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

$$\displaystyle T(-1,y)=T(x,1)=2$$ B.C.1 and 2
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

The square also looses energy through a heat flux of 3 for all x when y=1. This is given by the following equation and is in the figure below for $$\displaystyle \Gamma_h$$.

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

$$\displaystyle q=\vec n\cdot\mathbf K(\vec\nabla T)|_{x,y=1}$$ B.C.3
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Where n is the surface normal. Lastly, for all y at x=-1 the boundary of the square looses heat through Newton's Law of Cooling. This is given in the equation and figure below for $$\displaystyle\Gamma_H$$.

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

$$\displaystyle q=H(u(x=-1,y)-u_\infty)$$ B.C.4
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Where H is the convective coefficient and $$\displaystyle u_\infty$$ is the temperature of the fluid outside the square far away.

Part 1: Develop Weak Form
Develop the Weak form of General 2-Dimensional Model:

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

$$\displaystyle \frac{\partial}{\partial x_i}\cdot K_{ij} \nabla_j u+f(x,y)=\rho c\frac{\partial u}{\partial t}$$ (6.10.1)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

with

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

$$\displaystyle \mathbf{K}=\mathbf{I}$$
 * style="width:10%; |
 * style="width:10%; |

$$\displaystyle f(\underline{x},t)=0$$

$$\displaystyle \frac{\partial u}{\partial t}=0$$

$$\displaystyle g=2,h=3,H=.5,u_\infty =0.1$$

(6.10.2)
 *

 * }

Solution
To develop the weak form one must first develop the strong form. Performing an energy balance on a differential portion of the square gives the following figure:

Performing a heat balance gives:

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

$$\displaystyle H_{in}-H_{out}+H_{generated}=H_{stored}$$
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

where

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

$$\displaystyle H_{in}=q_x|x+q_y|y$$
 * style="width:10%; |
 * style="width:10%; |

$$\displaystyle H_{out}=q_x|_{x+dx}+q_y|_{y+dy}$$

$$\displaystyle H_{generated}=f(x)$$

$$\displaystyle H_{stored}=\rho c \frac{\partial u}{\partial t}$$
 *

 * }

giving

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

$$\displaystyle q_x|_x+q_y|_y-q_x|_{x+dx}-q_y|{y+dy}+f(x)=\rho c\frac{\partial u}{\partial t}$$
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Using a Taylor Series:

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

$$\displaystyle q_x|_{x+dx}=q_x(x+dx)=q_x(x)+\frac{\partial q_x}{\partial x}+\dots$$
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

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

$$\displaystyle q_y|_{y+dy}=q_y(y+dy)=q_y(y)+\frac{\partial q_y}{\partial y}+\dots$$
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Plugging these equations and canceling out terms:

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

$$\displaystyle -\frac{\partial q_x}{\partial x}-\frac{\partial q_y}{\partial y}+f(x,y)=\rho c\frac{\partial u}{\partial t}$$
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Using the constitutive equation for energy transfer in this medium, Fourier's Law:

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

$$\displaystyle \vec q=-\mathbf K\nabla T$$
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Substituting in this equation, which is a vector equation and could be written $$\displaystyle q_{x_i}=-K_{i,j}\cdot\nabla_j T$$ gives:

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

$$\displaystyle \frac{\partial}{\partial x_i} K_{ij} \frac{\partial}{\partial x_j} u+f(x,y)=\rho c\frac{\partial u}{\partial t}$$
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Where we have written the solution in Einstein Notation and are using $$\displaystyle u$$ instead of $$\displaystyle T$$ for temperature. G2DM2.0/D1 gives the following where the Einstein Notation has been altered for the common terms in vector calculus divergence and gradient:

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

$$\displaystyle div(\underline{I}\cdot gradu)=0$$ (6.10.3)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

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

$$\displaystyle \vec n(\Gamma_h)\cdot\underline{I}\cdot gradu|_{\Gamma_h}=h=3$$ (6.10.4)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

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

$$\displaystyle u(\Gamma_g)=g=2$$ (6.10.5)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

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

$$\displaystyle H[u(\Gamma_H)-u_\infty ]=-\vec n(\Gamma_H)\cdot\underline{I}\cdot gradu(\Gamma_H)$$ (6.10.6)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Where $$\Gamma_g$$ represents the location on the boundary of the essential boundary condition, $$\Gamma_h$$ represents the location of the natural boundary condition, and $$\Gamma_H$$ represents the location of Newton's Law of Cooling boundary condition as stated previously.

Part 1: Weak Form
The first step to formulate the general weak form is to multiply by a weighting function $$w(x,y)$$ and integrate over the domain of the problem, represented by $$\Omega$$:

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

$$\displaystyle \int w(x,y)(\frac{\partial}{\partial x_i}\cdot K_{ij} \frac{\partial}{\partial x_j} u+f(x,y)-\rho c\frac{\partial u}{\partial t})d\Omega=0$$ (6.10.7)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Multiplying the natural and Newton's Law of Cooling Boundary conditions by the same weighting function:

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

$$\displaystyle w(\Gamma_h)[n_i\cdot(\Gamma_h)K_{ij}\frac{\partial}{\partial x_j}u|_{\Gamma_h}-3]=0$$ (6.10.8)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

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

$$\displaystyle w(\Gamma_H)[H[u(\Gamma_H)-u_\infty]+n_i(\Gamma_h)K_{ij}\frac{\partial}{\partial x_j}u|_{\Gamma_H}]=0$$ (6.10.9)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Using integration by parts ($$u=w$$ and $$dv=div(\underline{I}\cdot gradu)d\Omega$$) as used here and the fundamental theorem of calculus gives:

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

$$\displaystyle \int \frac{\partial}{\partial x_i}(wK_{ij}\frac{\partial u}{\partial x_j})d\Omega-\int\frac{\partial w}{\partial x_i}K_{ij}\frac{\partial u}{\partial x_j}d\Omega=0$$ (6.10.10)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Applying the divergence theorem on the first term:

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

$$\displaystyle \int \frac{\partial}{\partial x_i}(wK_{ij}\frac{\partial u}{\partial x_j})d\Omega=\int wn_iK_{ij}\frac{\partial }{\partial x_j}d(\partial\Omega)$$ (6.10.11)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Since $$\partial\Omega$$ is the boundary of the domain, it can be split up in

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

$$\displaystyle \int wn_iK_{ij}\frac{\partial }{\partial x_j}d(\partial\Omega)=\int_{\Gamma_g}wn_iK_{ij}\frac{\partial u}{\partial x_j}d\Gamma_g + \int_{\Gamma_h}wn_i K_{ij} \frac{\partial u}{\partial x_j}d\Gamma_h+\int_{\Gamma_H}w n_i K_{ij}\frac{\partial u}{\partial x_j}d\Gamma_H$$ (6.10.12)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

We will take advantage of Eq. 7.10.12 to remove the unknown flux on the essential boundary condition the arbitrary weighting function will be set to zero on the essential boundary, $$w(\Gamma_g)=0$$. For now, plugging everything in gives the general weak form:

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

$$ \displaystyle \int wn_iK_{ij}\frac{\partial }{\partial x_j}d(\partial\Omega)-\int_\Omega\frac{\partial w}{\partial x_i}K_{ij}\frac{\partial u}{\partial x_j}d\Omega-\int_\Omega w\rho c\frac{\partial u}{\partial t}d\Omega=0
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$     (Eq 6.10.13)
 *

 * }

Part 2: Semi Discrete Equations
To develop the semidiscrete version first form the following functions with the LIBF:

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

$$\displaystyle u\simeq u^h = \sum_{n=1}^{a*b}d_n\prod_{(j=1,j\neq n)}^a \prod_{(k=1,k\neq n)}^b\frac{(x-x_j)(y-y_k)}{(x_n-x_j)(y_n-y_k)}=\sum_{n=1}^a*d_n*N_n(x)$$ (6.10.14)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

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

$$\displaystyle w\simeq w^h = \sum_{m=1}^{a*b}c_m\prod_{(j=1,j\neq m)}^a \prod_{(k=1,k\neq m)}^b\frac{(x-x_j)(y-y_k)}{(x_m-x_j)(y_m-y_k)}=\sum_{m=1}^a*c_m*N_m(x)$$ (6.10.15)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Substituting Eq.'s 7.10.2, 7.10.12, 7.10.14, and 7.10.15 into the weak form:

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

$$\displaystyle -\int_{\Gamma_h}\sum_{m}c_mN_m(x,1)3d\Gamma_h-\int_{\Gamma_H}\sum_{m}c_mN_m(-1,y)0.5[\sum_{n}d_nN_n(-1,y)-.1]d\Gamma_H-\int_\Omega\frac{\partial \sum_{m}c_mN_m}{\partial x_i}I_{ij}\frac{\partial \sum_{n}d_nN_n}{\partial x_j}d\Omega=0 $$ (6.10.16)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

Factoring out the constant $$c_m$$ and recognizing subscripts m,n,i, and j mean to sum over:

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

$$\displaystyle c_m\Bigg [-\int_{\Gamma_h}N_m(x,1)3d\Gamma_h-\int_{\Gamma_H}N_m(-1,y)0.5[d_nN_n(-1,y)-.1]d\Gamma_H-\int_\Omega\frac{\partial N_m}{\partial x_i}I_{ij}\frac{\partial d_nN_n}{\partial x_j}d\Omega\Bigg]=0   $$ (6.10.17)
 * style="width:10%; |
 * style="width:10%; |
 *

 * }

We will now move to put this into matrix form similarly to this previous work. In order to do this the numbering of the nodes must begin on the essential boundary. This allows partitioning of the matrices as before, with essential and free components. By numbering the nodes on the essential boundary first, this means the basis functions, $$\displaystyle N_m,N_n$$ first functions will be on the essential boundary (m and n correspond to the node number). With this recognition we can form the following:

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

$$\displaystyle \vec C=\begin{bmatrix} \vec C_E \\ \vec C_F\end{bmatrix}=\begin{bmatrix} 0 \\ \vec C_F\end{bmatrix} $$     (6.10.18)
 *

 * }

The zero in equation comes from satisfying the constraint breaking solution where the coefficient at the essential boundary condition for the only nonzero basis function coefficient on the weighting function is zero. Hence the subscript E in equation 5.3.11 is for essential. Looking at the coefficients for uh:

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

$$\displaystyle \vec d=\begin{bmatrix} \vec d_E \\ \vec d_F\end{bmatrix}=\begin{bmatrix} g \\ \vec d_F\end{bmatrix} $$     (6.10.19)
 *

 * }

Likewise here, g is the value at the essential boundary condition. Since there is only one nonzero basis function here, the coefficient (since the basis function is one) must be equal to the essential boundary condition value. The subscript F can be thought of as free and will be determined.

Now defining the following

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

$$\displaystyle \tilde k(N_m,N_n)=\int_0^1\int_0^1(\frac{d}{dx}N_m)I_{ij}(\frac{d}{dx}N_n)dxdy+\int_{\Gamma_H}.5Nn(-1,y)d\Gamma_H $$     (6.10.20)
 *

 * }

This shows how the far left integral term of equation 5.3.10 can be written in components of a matrix. Note the additional term that contains the constants $$\displaystyle d_n$$ which is what the K matrix will multiply. Additionally defining:

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

$$\displaystyle K_{EE}=[\tilde k(N_{m=1\rightarrow a+b-1},N_{n=1\rightarrow a+b-1})] $$     (6.10.21)
 *

 * }

Where m increments from 1 to a+b-1 (a is number of x nodes and b is number of y nodes) and n increments from 1 to a+b-1. These end values (a+b-1) correspond to the number of nodes that are found on the essential boundary, hence being part of the $$\displaystyle K_{EE}$$ matrix. Continuing on with the defining:

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

$$\displaystyle K_{EF}=[\tilde k(N_{m=1\to a+b-1},N_{n=a+b\to a*b}] $$     (6.10.22)
 *

 * }

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

$$\displaystyle K_{FE}=[\tilde k(N_{m=a+b\to a*b},N_{n=1\to a+b-1})] $$     (6.10.23)
 *

 * }

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

$$\displaystyle K_{FF}=[\tilde k(N_{m=a+b\to a*b},N_{n=a+b\to a*b})] $$     (6.10.24)
 *

 * }

Giving:

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

$$\displaystyle \mathbf K=\begin{bmatrix} \mathbf K_{EE} \mathbf K_{EF} \\ \mathbf K_{FE} \mathbf K_{FF}\end{bmatrix} $$     (6.10.25)
 *

 * }

It can be seen that the subscripts E correspond the the essential nodes, and consequently the basis functions on the essential nodes (1 to a+b-1 where minus one takes into account double counting of the point shared with the connecting essential boundary borders). The F subscript corresponds to the free nodes (a+b to a*b) or free basis functions where the coefficient must still be determined.

Lastly, combining the constants into:
 * {| style="width:100%" border="0"

$$\displaystyle \mathbf F = [\tilde f_m] $$     (6.10.26)
 *

 * }

Where,

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

$$\displaystyle \tilde f_m=-\int_{\Gamma_h}3N_m(x,1)d\Gamma_h + \int_{\Gamma_H}.5N_m(0,y)d\Gamma_H $$     (6.10.27)
 *

 * }

Putting it all together:

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

$$\displaystyle \vec C \cdot \mathbf K\mathbf d=\begin{bmatrix}0\\\vec C_F\end{bmatrix}\cdot \Big (\begin{bmatrix} \mathbf K_{EE} \mathbf K_{EF} \\ \mathbf K_{FE} \mathbf K_{FF}\end{bmatrix}\begin{bmatrix}\mathbf g\\\vec d_F\end{bmatrix}-\begin{bmatrix}\mathbf F_E\\ \mathbf F_F\end{bmatrix} \Big )=0 $$     (6.10.28) Now, our weighting function is arbitrary so what is in parenthesis corresponding to CF must equal zero. Writing this out explicitly as a simple example for four Nodes gives the following:
 *

 * }

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

$$\displaystyle \begin{bmatrix}K_{11}K_{12}K_{13}K_{14}\\K_{21}K_{22}K_{23}K_{24}\\K_{31}K_{32}K_{33}K_{34}\\K_{41}K_{42}K_{43}K_{44}\end{bmatrix}\begin{bmatrix}g\\g\\g\\d_4\end{bmatrix}=\begin{bmatrix}f_1+r_1\\f_2+r_2\\f_3+r_3\\f_4\end{bmatrix} $$     (6.10.29)
 *

 * }

The first 3 g's correspond to the essential boundary, so in order to solve for our unknown d4, in terms of our knowns, we rearrange to get the following:

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

$$\displaystyle \begin{bmatrix}K_{44}\end{bmatrix}\begin{bmatrix}d_4\end{bmatrix}=\begin{bmatrix}f_4-K_{41}g-K_{42}g-K_{43}g\end{bmatrix} $$     (6.10.30)
 *

 * }

Part 3: MATLAB Solution and Plot
The following is the MATLAB code used to find the coefficients using Part 2's steps:

This gave the following plot with 2 nodes in the x and y direction (4 nodes total) and a Temperature at x=y=.5 of 2.6437:

For 3 nodes in the x and y direction (9 nodes total) and a Temperature at x=y=.5 of 4.7215:

Lastly, 4 nodes in the x and y direction (16 nodes total) and a Temperature at x=y=.5 of 7.0725:

Comments
Obtained a plot qualitatively similar to previous homework, but quantitatively closer to the correct solution (assuming a large negative temperature would be incorrect). Qualitatively the graph behaves as it should. The boundary for Newton's Law of Cooling, or convective heat loss, shows a slower decrease in temperature when compared to the conductive boundary. This makes sense since the convective coefficient of the fluid is .5 and serves to slightly insulate the heat loss in comparison to the conductive boundary condition. However, At the edge where the convective and conductive boundary conditions are applied the temperture still drops below zero, which is very possibly incorrect since the temperature is specified on the convective boundary as .1. To drop below the temperature of the surrounding fluid when two of the objects sides are held at a higher temperature than the surrounding fluid seems highly unlikely. But, we were unable to run down the error in the code.

title: User:Eml5526.s11.team7/HW7


text: File:Eml5526.s11.team7.hw7.2a n4sim.jpg

title: File:Eml5526.s11.team7.hw7.2a n4sim.jpg


title: File:Fe1.s11.team3.hw7.1.b.gif


title: File:Fe1.s11.team3.hw7.1.a.gif


title: File:Fe1.s11.team3.hw7.2.a.gif


title: File:Fe1.s11.team3.hw7.2.b.gif