User:Fem.tm7/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)
 * <p style="text-align:right">
 * }


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


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

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

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 3x3 Nodes using adaptive Simpson quadrature



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

Results with 4x4 Nodes using adaptive Simpson quadrature



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

Results with 5x5 Nodes using adaptive Simpson quadrature



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

This method shows better convergence than the above method using symbolic integration.

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

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

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

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

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

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

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

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

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


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

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

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

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

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


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


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

$$\displaystyle u(\Gamma_g)=g=2$$ (6.10.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


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

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

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


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

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

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

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

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

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


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

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

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

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

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

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

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

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


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


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

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

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

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

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


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

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

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.