User:Eml5526.s11.team2.greatsemesterHW7

=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 $$

2) Dynamic Case(Transient)
$$\rho c=3$$

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

(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 $$

(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 $$

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%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle Md^{(s)}+Kd=F $$     (7.1.14) For static case
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| 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" |
 * <p style="text-align:right">
 * }


 * {| 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" |
 * <p style="text-align:right">
 * }

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)
 * <p style="text-align:right">
 * }

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

=Problem 7.2: 2-D Vibrating Membrane Using Static and Dynamic LIBF Analysis=

Given: The 2-D Vibrating Membrane Data Set
Given the data set provided below.


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

$$ \begin{align} \displaystyle T \cdot div(\nabla \mathbf{u}) + \mathbf{f} &= \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} \\ T &:= \mathbf{Tension} \left( \frac{force}{length} \right) \\ \mathbf{u}(\mathbf{x},t) &:= \mathbf{Transverse \quad Displacement} \\ \mathbf{f}(\mathbf{x},t) &:= \mathbf{Distribute \quad Transverse \quad Force} \\ \rho(\mathbf{x}) &:= \mathbf{Mass \quad Density}\left( \frac{mass}{area}\right) \\ \Omega &:= \mathbf{Membrane \quad Boundary}\\ d\Omega &:= \mathbf{Displacement \quad on \quad the \quad Membrane \quad Boundary}\\ \end{align} $$     (7.2.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

1) Static (Steady State):

 * Find the solution for the static vibrating membrane $$ \displaystyle \mathbf{u}_s^{h} $$ to $$ 10^{-6} $$ accuracy for the problem formulated below.


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

$$ \begin{align} \displaystyle T \cdot div(\nabla \mathbf{u}) + \mathbf{f} &= 0 \\ T &= 4 \\ \mathbf{f}(\mathbf{x},t) &= 1\quad in \quad \Omega=\square \\ g &= 0\quad on \quad d\Omega \implies \frac{d^{2}\mathbf{u}}{dt^{2}}= 0 \\ \end{align} $$     (7.2.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

2) Dynamic (Transient): $$ \rho=3 $$

 * Find the solution for the dynamic vibrating membrane $$ \displaystyle \mathbf{u}_d^{h} $$ to $$ 10^{-6} $$ accuracy for the problem formulated below.


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

$$ \begin{align} \displaystyle T \cdot div(\nabla \mathbf{u}) + \mathbf{f} &= \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} \\ T &= 4 \\ \rho &= 3 \\ \mathbf{f}(\mathbf{x},t) &= 0\quad in \quad \Omega=\square \\ \mathbf{u}(\mathbf{x},t=0) &= 0 \quad \forall \mathbf{x} \in \Omega\\ \mathbf{\dot{u}}(\mathbf{x},t=0) &= 0 \quad \forall \mathbf{x} \in \Omega\\ \end{align} $$     (7.2.3)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

2a) Free Vibration: Solve General Eigenvalue Problem (EVP)

 * The general EVP is given by,
 * {| style="width:100%" border="0"

$$ \displaystyle \mathbf{K \Phi} = \lambda \mathbf{M \Phi} $$     (7.2.4)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * Where each of the component to the EVP are defined here,
 * {| style="width:100%" border="0"

$$ \begin{align} \displaystyle \quad \mathbf{\Phi} &:= \mathbf{Eigenvector} \\ \mathbf{\lambda}&:= \mathbf{Eigenvalue} \\ \mathbf{M}&:= \mathbf{Mass(Inertial)} \ \mathbf{Matrix} \\ \mathbf{K}&:= \mathbf{Stiffness} \ \mathbf{Matrix} \\ \end{align} $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * Determine the first 3 eigenpairs $$\displaystyle \left( \lambda_{i},\mathbf{\phi _i} \right) $$ where the fundamental frequency $$\displaystyle \omega_1 $$ and fundamental period $$\displaystyle T_1 $$ are given by


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

$$ \begin{align} \displaystyle \omega _1 &= \left( \lambda _1 \right)^{\frac{1}{2}} \quad s.t. \quad {f_1} = \frac \\ and \quad {T_1} &= \frac{1} \end{align} $$     (7.2.5)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * Plot $$\displaystyle \mathbf{u}^{h}(\mathbf{0},t)$$ vs. t for $$\displaystyle t \in [0,2T_1] $$.
 * Produce a movie of the vibrating membrane.

3) Dynamic (Transient): Similar to Part 2 but Using Asymmetric $$ u_0 $$

 * Find the solution for the dynamic vibrating membrane $$ \displaystyle \mathbf{u}_d^{h} $$ to $$ 10^{-6} $$ accuracy for the problem formulated below.


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

$$ \begin{align} \displaystyle T \cdot div(\nabla \mathbf{u}) + \mathbf{f} &= \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} \\ T &= 4 \\ \rho &= 3 \\ \mathbf{f}(\mathbf{x},t) &= 0\quad in \quad \Omega=\square \\ \mathbf{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\\ \mathbf{\dot{u}}(\mathbf{x},t=0) &= 0 \quad \forall \mathbf{x} \in \Omega\\ \end{align} $$     (7.2.6)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

3a) Free Vibration: Solve General Eigenvalue Problem (EVP)

 * The general EVP given in (7.2.4).
 * Determine the first 3 eigenpairs $$\displaystyle \left( \lambda_{i},\mathbf{\phi _i} \right) $$ where the fundamental frequency $$\displaystyle \omega_1 $$ and fundamental period $$\displaystyle T_1 $$ are given in (7.2.5).
 * Plot $$\displaystyle \mathbf{u}^{h}(\mathbf{0},t)$$ vs. t for $$\displaystyle t \in [0,2T_1] $$.
 * Produce a movie of the vibrating membrane.

1) Solution to the Static Case
First we will recall that the governing equation in (7.2.2) is,

Our governing PDE can be rewritten as the Weak Form shown here, Next we will integrate by parts, rearrange terms, apply Gauss' Theorem and arrive at our final formulation.

It is customary to represent the final expression in (7.2.8) as a more concise formulation shown here,
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} \tilde{m}(w,u) + \tilde{k}(w,u) &= \tilde{f}(w) \quad \forall \ w=0 \ on \ \Gamma_g \\ \tilde{m}(w,u) &= \int\limits_\Omega {w\rho \frac} d\Omega =0 \\ \tilde{k}(w,u) &= \int\limits_\Omega {T\nabla w\nabla u} d\Omega \\ \tilde{f}(w) &= - \int_{\Gamma h} {whd{\Gamma _h}}  + \int\limits_\Omega  {wf} d\Omega \\ \end{align} $$     (7.2.9)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }



2) Solution to Eigenvalue Problem: Separation of Variables
The governing PDE is a function of $$ \displaystyle x, y,$$ and $$ \displaystyle t $$ we can rewrite $$ \displaystyle u(x,y,t)$$ as a product, where we seek a solutions of the form,


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

$$ \displaystyle u (x,y,t) = {X}(x) {Y}(y) {T}(t) $$     (7.2.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Our governing PDE is for the static membrane is shown here,
 * {| style="width:100%" border="0"

$$ \displaystyle T\left( {\frac + \frac} \right) = \rho \frac $$     (7.2.13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Upon rearrangement we obtain the following form which is more suitable for applying Separation of Variables,
 * {| style="width:100%" border="0"

$$ \displaystyle \left( {\frac + \frac} \right) = \frac{\rho }{T}\frac $$     (7.2.14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Substituting the appropriate derivatives of (7.2.12) into (7.2.14) gives,
 * {| style="width:100%" border="0"

$$ \displaystyle \left( {XYT + YXT} \right) = \frac{\rho }{T}T''XY $$ (7.2.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Dividing the LHS and RHS by the product $$ \displaystyle X(x)Y(y)T(t) $$ and setting the equation equal to a constant $$ \displaystyle -\omega^2 $$, (In this case the constant is set up to work with the natural frequency calculation):
 * {| style="width:100%" border="0"

$$ \displaystyle \frac{X} + \frac{Y} = \frac{\rho }{T}\frac{T} = - {\omega ^2} $$     (7.2.16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Letting $${X} $$ and $$ {Y} $$ be equal to the arbitrary constants $$\displaystyle k_1 $$ and $$\displaystyle k_2 $$ we get the following ODEs shown in (7.2.17) and (7.2.18). One can consider these ODEs as being constrained by equation (7.2.19).
 * {| style="width:100%" border="0"

$$ \displaystyle {\frac{X}} = - {k_1}^2 $$     (7.2.17)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle {\frac{Y}} = - {k_2}^2 $$    (7.2.18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$ \displaystyle {\omega ^2} = {k_1}^2 + {k_2}^2 $$    (7.2.19)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Upon rearrangement of (7.2.17) and (7.2.18) $$ \displaystyle X(x) $$ and $$ \displaystyle Y(y) $$ are represented in the two second order homogeneous ODEs.
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{align} {X''} + {k_1}^2{X} &= 0 \\ {Y''} + {k_2}^2{Y} &= 0 \\ \end{align} $$    (7.2.20)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

$$ \displaystyle X(x) $$ and $$ \displaystyle Y(y) $$ have the same general solutions,
 * {| style="width:100%" border="0"

$$ \displaystyle {X}=Acos({k_1}x)+Bsin({k_2}x) $$     (7.2.21)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For the simplicity of imposing the BCs let us assume the domain is $$ \displaystyle {x,y}\in [0,L] $$. When we impose the BC by evaluating at 0, the sin term drops out therefore making $$ \displaystyle A=0 $$. Next observe at L, $$\displaystyle k_{1}L = n\pi $$ and $$\displaystyle k_{2}L = m\pi $$.

Similarly, the second order homogeneous ODE for $$ \displaystyle T(t) $$ is,
 * {| style="width:100%" border="0"

$$ \displaystyle {T''} + \frac{T}{\rho }{\omega ^2}{T} = 0 $$     (7.2.22)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Which has the following general solution,
 * {| style="width:100%" border="0"

$$ \displaystyle {T} = C\cos \left( {\omega\sqrt {\frac{T}{\rho }} t} \right) + D\sin \left( {\omega\sqrt {\frac{T}{\rho }} t} \right) $$     (7.2.23)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Therefore, the natural frequency is $$ \displaystyle \omega\sqrt {\frac{T}{\rho }} $$. Expressing the natural frequeny in terms of equation (7.2.19) gives us,
 * {| style="width:100%" border="0"

$$ \displaystyle \omega_{m,n}=\sqrt {\frac{T}{\rho }}\sqrt {k_1^2+k_2^2}=\sqrt {\frac{T}{\rho }\left[ {{{\left( {\frac{{n\pi }}{L}} \right)}^2} + {{\left( {\frac{{m\pi }}{L}} \right)}^2}} \right]} $$     (7.2.24)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Substituting the numerical values for density, temperature, and BC's from the problem statement and letting $$ \displaystyle m=n=1 $$, $$ \displaystyle m=2 $$ and $$ \displaystyle n=1 $$, and $$ \displaystyle m=n=2 $$,


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

$$ \displaystyle \begin{align} {\omega _{1,1}} &= 5.13 \\ {\omega _{2,1}} &= 8.11 \\ {\omega _{2,2}} &= 10.26 \\ \end{align} $$   (7.2.25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

=Problem 7.3: 2-D Static Vibrating Membrane Analysis Using Quadratic and Rectangular Meshes=

Given: The 2-D Vibrating Membrane Data Set
Given the data set provided below.


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

$$ \begin{align} \displaystyle T \cdot div(\nabla \mathbf{u}) + \mathbf{f} &= \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} \\ T &:= \mathbf{Tension} \left( \frac{force}{length} \right) \\ \mathbf{u}(\mathbf{x},t) &:= \mathbf{Transverse \quad Displacement} \\ \mathbf{f}(\mathbf{x},t) &:= \mathbf{Distribute \quad Transverse \quad Force} \\ \rho(\mathbf{x}) &:= \mathbf{Mass \quad Density}\left( \frac{mass}{area}\right) \\ \Omega &:= \mathbf{Membrane \quad Boundary}\\ d\Omega &:= \mathbf{Displacement \quad on \quad the \quad Membrane \quad Boundary}\\ \end{align} $$     (7.3.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

1) Static (Steady State):

 * Let $$ \displaystyle \Omega = \bigcirc $$ be a circle of unit radians.
 * Find the solution for the static vibrating membrane $$ \displaystyle \mathbf{u}_s^{h} $$ to $$ 10^{-6} $$ accuracy at the center for the problem formulated below.


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

$$ \begin{align} \displaystyle T \cdot div(\nabla \mathbf{u}) + \mathbf{f} &= 0 \\ T &= 4 \\ \mathbf{f}(\mathbf{x},t) &= 1\quad in \quad \Omega= \bigcirc \\ g &= 0\quad on \quad \Gamma_g=d\Omega \implies \frac{d^{2}\mathbf{u}}{dt^{2}}= 0 \\ \Omega &= \bigcirc \\ \end{align} $$     (7.3.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Solution: Matlab
After modifying code to no avail, I was fortunate enough to use the abilities of wikiversity to find a similar code and then modified it. The code used to solve this problem was referenced from Deshpande. I made modifications to the code to allow for the use of a greater number of elements. Here I used 16 elements instead of the 9 used by team 6. In modifying the code I also used code from Fish and Belytschko, A First Course in Finite Elements. The code was found through their website and was in chapter 9 [|F&B]. The code used is as follows

The graphs produced by the code are



=Problem 7.4=

Given: The Governing Equation and Boundary Conditions
Governing equation was derived in meeting 32-4.
 * $$\displaystyle Q_{1}$$:Heat flow into w throgh $$\displaystyle \partial w$$
 * $$\displaystyle Q_{2}$$:Heat generated in w by heat source.
 * $$\displaystyle Q_{3}$$:Heat in w due to change in term u.


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

$$  \displaystyle Q_{1}=-\int_{\partial w}\mathbf{q}\mathbf{n}d(\partial w)=\int_{w}div\mathbf{q}dw $$     (7.4.1)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle Q_{2}=\int_{w}f(\mathbf{x},t)dw $$     (7.4.2)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle Q_{3}=\int_{w}\rho c\frac{\partial u}{\partial t}dw $$     (7.4.3)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle Q_{1}+Q_{2}=Q_{3} $$     (7.4.4)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \mathbf{q}=-\mathbf{K}grad(u) $$     (7.4.5)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle div(\mathbf{K}grad(u))+f=\rho c\frac{\partial u}{\partial t} $$ (7.4.6)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Boundry Conditions

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

$$  \displaystyle u(\mathbf{x},t)]_{\Gamma _{g}}=g(\mathbf{x},t) $$     (7.4.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle (q.n)]_{\Gamma _{h}}=h(\mathbf{x},t) $$     (7.4.8)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle (q.n)]_{\Gamma _{H}}=H(u-u_{\infty }) $$     (7.4.9)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

WRF

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

$$  \displaystyle \int_{\Omega }w[div(Kgrad(u))+f-\rho c\frac{\partial u}{\partial t}]d\Omega =0 $$     (7.4.10)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} K=  & \int_{\Omega }wdiv(Kgrad(u))d\Omega =\frac{\partial }{\partial x_{i}}(K_{ij}\frac{\partial u}{\partial x_{j}})\\ & \int_{\Omega }n_{i}wK_{ij}\frac{\partial u}{\partial x_{j}}d(\partial\Omega)-\int_{\Omega }\frac{\partial w}{\partial x_{i}}K_{ij}\frac{\partial u}{\partial x_{j}}d(\Omega) \end{align} $$     (7.4.11)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle K=-\int_{\Gamma _{g}}wnqd\Gamma _{g}-\int_{\Gamma _{h}}wnqd\Gamma _{h}-\int_{\Gamma _{H}}wnqd\Gamma _{H}+\int_{\Omega }grad(w).K.grad(u)d\Omega $$     (7.4.12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

We deliberately chose w i.e. i vanishes on essential boundry. So the first term vanishes too. In the first term we have unknown flux that we can not deal with right now. But other fluxes are defined.We can write our continuous form of components as:


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

$$  \displaystyle \tilde{m}(w,u)=\int_{\Omega}w\rho c\frac{\partial u}{\partial t}d\Omega $$     (7.4.13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \tilde{k} (w,u)=\int_{\Omega}grad(w)Kdiv(u)d\Omega+\int_{\Gamma_{H}}wHud\Gamma_{H} $$     (7.4.14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{f}(w,u)=\int_{\Omega}wfd\Omega-\int_{\Gamma_{h}}whd\Gamma_{h}+\int_{\Gamma_{H}}wHu_{\infty }d\Gamma_{H} $$     (7.4.15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

!Note this time f depeneds on both w and u.

Discrete WF

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

$$  \displaystyle w^{app}=\sum_{I=1}^{n}N(x,y)_{I}c_{I} $$     (7.4.16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle u^{app}=\sum_{J=1}^{n}N(x,y)_{J}d_{J} $$     (7.4.17)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle grad(u)=grad(N(x,y)_{I})d_{I}=B(x,y)_{I}d_{I} $$     (7.4.18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{m}(w,u)+\tilde{k}(w,u)=\tilde{f}(w,u) $$     (7.4.19) Now lets substitute everything in the equation above.
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \sum_{I=1}^{n}c^{T}\left [\sum_{J=1}^{n}\left \{ \int_{\Omega }\rho cN_{I}N_{J}d_{I}^{s}d\Omega+\int_{\Omega}B_{I}KB_{J}d_{I}d\Omega -\int_{\Omega}N_{I}fd\Omega+\int_{\Gamma_{h}}N_{I}hd\Gamma_{h}+\int_{\Omega}N_{I}N_{J}(\Gamma_{H})H(d(\Gamma_{H})-d_{\infty })\right \} \right ]=0 $$     (7.4.20)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{d}}=\begin{bmatrix} \bar{de}\\ df
 * style="width:95%" |
 * style="width:95%" |

\end{bmatrix} $$     (7.4.21)
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{w}}=\begin{bmatrix} 0\\ wf
 * style="width:95%" |
 * style="width:95%" |

\end{bmatrix} $$     (7.4.22)
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{K}}=\begin{bmatrix} K_{E} &K_{EF} \\ K_{EF}^{T} & K_{F} \end{bmatrix} $$     (7.4.23)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Example of calculation of elements of K & M matrices when only one node lie on essential boundary.
 * {| style="width:100%" border="0"

$$  \displaystyle\begin{align} & K_{E}=\int_{\Omega }B_{1}.K.B_{1}d\Omega+\int\limits_ {{N_1}H{N_1}d{\Gamma _H}} &\\ & K_{EF}=\int_{\Omega }B_{1}.K.B_{J}d\Omega+\int\limits_ {{N_1}H{N_j}d{\Gamma _H}} & J=2,n\\ & K_{FF}=\int_{\Omega }B_{I}.K.B_{J}d\Omega+\int\limits_ {{N_i}H{N_j}d{\Gamma _H}} & I,J=2,n
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$     (7.4.24)
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \tilde{\mathbf{M}}=\begin{bmatrix} M_{E} &M_{EF} \\ M_{EF}^{T} & M_{F} \end{bmatrix} $$     (7.4.25)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle\begin{align} & M_{E}=\int_{\Omega }N_{1}\rho cN_{1}d\Omega&\\ & M_{EF}=\int_{\Omega }N_{1}\rho cN_{J}d\Omega& j=2,n\\ & M_{FF}=\int_{\Omega }N_{I}\rho cN_{J}d\Omega& i,j=2,n
 * style="width:95%" |
 * style="width:95%" |

\end{align} $$     (7.4.26)
 * <p style="text-align:right">
 * }


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

$$ \underline {\tilde c} = \left[ {\begin{array}{ccccccccccccccc} {\underline } \\ \hline \end{array}} \right] $$     (7.4.27)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \underline {\tilde c} = \left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right] $$     (7.4.28)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Structure of d matrix
 * {| style="width:100%" border="0"

$$  \displaystyle \underline {\tilde d} = \left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right] $$     (7.4.29)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Structure of K matrix for this type of problem
 * {| style="width:100%" border="0"

$$  \displaystyle {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T{\left[ {\begin{array}{ccccccccccccccc} &\vline & && \\ \hline &\vline & && \\ &\vline & && \\ &\vline & && \end{array}} \right]_{\tilde n \times \tilde n}}{\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{\tilde n \times 1}}}+{\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T{\left[ {\begin{array}{ccccccccccccccc} {}&\vline & {}&{}&{} \\ \hline {}&\vline & {}&{}&{} \\ {}&\vline & {}&{\underline K _{HH}^H}&{} \\ {}&\vline & {}&{}&{} \end{array}} \right]_{\tilde n \times \tilde n}}{\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{\tilde n \times 1}}} $$     (7.4.30)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Components of force matrix
 * {| style="width:100%" border="0"

$$  \displaystyle {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T\left[ {\begin{array}{ccccccccccccccc} {\underline F _E^f} \\ \hline {\underline F _h^f} \\ {\underline F _H^f} \\ {\underline F _R^f} \end{array}} \right] + {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T\left[ {\begin{array}{ccccccccccccccc} {} \\ \hline {\underline F _h^h} \\ {} \\  {} \end{array}} \right] + {\left[ {\begin{array}{ccccccccccccccc} \\ \hline \\   \\ \end{array}} \right]_{_{1 \times \tilde n}}}^T\left[ {\begin{array}{ccccccccccccccc} {} \\ \hline {} \\  {\underline F _H^H} \\ {} \end{array}} \right] $$     (7.4.31)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


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

$$  \displaystyle \begin{matrix} Temperature & at & global & nodes \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }




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

$$  \displaystyle \begin{matrix} Temperature & Contours \end{matrix} $$
 * style="width:95%" |
 * style="width:95%" |
 * }



=Contributing Members & Referenced Lecture=