Homework 7

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

## Given: Problem defined on Mtg 39 (a)

Strong Form

 $\displaystyle \underbrace{div(\boldsymbol{\kappa}grad u)}_{\frac{\partial }{\partial x_{i}}(\kappa _{ij}\frac{\partial u}{\partial x_{j}})}+f=\rho c\frac{\partial u}{\partial t}$ (Eq 1.1)

for the static case

 \begin{align} \text{1)}\quad \boldsymbol{\kappa}\left( {x,t} \right) = \mathbf{I} \\ f\left( {x,t} \right) = 1 \\ \frac{\partial u}{\partial t}\left( {x,t} \right) = 0 \quad \Omega = \square\\ \end{align} (Eq 1.2)

Initial condition for dynamic case

 \begin{align} \text{2)}\quad \displaystyle \quad u({\color{Red}}\underline{x},t=0)=xy \quad \forall {\color{Red}}\underline{x} \quad \in \Omega \\ \end{align} (Eq 1.3)

Boundary conditions for Dynamic case

 \begin{align} \text{2a)}\quad f({\color{Red}}\underline{x},t)=0 \quad \forall {\color{Red}}\underline{x} \quad \in \Omega \quad \forall \quad t>0 \\ g({\color{Red}}\underline{x},t)=2 \quad on \Gamma _{g}=\partial \Omega \\ \end{align} (Eq 1.4)
 \begin{align} \text{2b)}\quad f({\color{Red}}\underline{x},t)=1 \quad \forall {\color{Red}}\underline{x} \quad \in \Omega \quad \forall \quad t>0 \\ g({\color{Red}}\underline{x},t)=2 \quad on \Gamma _{g}=\partial \Omega \\ \end{align} (Eq 1.5)

## 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

## Solution

### Solution for static case

#### ii. Finding approximate solution

2DLIBF are given by

 \displaystyle \begin{align} &N_{I}(x,y)=L_{i,m}(x)L_{j,n}(y)\\ & where \quad I=i+(j-1)m\\ \end{align}

Eq 1.6

$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

 $L_{i,m}(x) = \prod_{k = 1}^{m}\frac{x-x_{k}}{x_{i}-x_{k}}\;\; k\neq i$; where $\displaystyle m$ is number of nodes.

Eq 1.7

 $L_{j,n}(y) = \prod_{k = 1}^{m}\frac{y-y_{k}}{y_{j}-y_{k}}\;\; k\neq j$; where $\displaystyle n$ is number of nodes.

Eq 1.8

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:

 $\displaystyle \tilde{\kappa}(\omega,u)=\int_{\Omega}\bigtriangledown \omega.\boldsymbol{\kappa }.\bigtriangledown ud\Omega$
 $\displaystyle \tilde{f}(\omega)=-\int_{\Gamma _{h}}\omega h d\Gamma _{h}+\int_{\Omega}\omega f d\Omega$
 $\displaystyle \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

 $\displaystyle \tilde{m}(\omega,u)=0$
 $\displaystyle \tilde{f}(\omega)= \int_{\Omega}\omega d\Omega$
 $\displaystyle \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 )$ :

 \displaystyle \begin{align} & u^{h}(x)=\sum_{j} d_{j} b_{j} \\ &\\ & w^{h}(x)=\sum_{i} c_{i} b_{i} \\ \end{align} (Eq 1.9)

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.

 $\sum_{i} c_i \left [ \sum_{j} \left \{ \tilde{M_{ij}} \ddot\tilde{d_j} \right \}+ \left \{ \tilde{K_{ij}} \tilde{d_j} \right \}- \tilde{F_i} \right ] = 0$ Eq 1.10

where,

Capacitance matrix

 $\displaystyle \overset{\tilde{\ }}{\mathop{m}}\,(w^h,u^h)=\int\limits_{\Omega }{w^h\rho c \frac{\partial u^h}{\partial t}d\Omega }$ $\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}\left( \int\limits_{\Omega }{{{b}_{i}}\rho c {{b}_{j}}}d\Omega \right)}}\dot{d_j}$ Eq 1.11

 $\displaystyle \overset{\tilde{\ }}{\mathop{M}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{{{b}_{i}}\rho c {{b}_{j}}}d\Omega$ Eq 1.12

Conductivity matrix

 $\displaystyle k(w^h,u^h)=\int\limits_{\Omega }{\nabla w^h.\kappa .\nabla u^hd\Omega}$ $\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}}\int\limits_{\Omega }{\left( \frac{\partial {{b}_{i}}}{\partial {{x}_{i}}}\kappa \frac{\partial {{b}_{j}}}{\partial {{x}_{j}}}d\Omega \right)}}{{d}_{j}}$ Eq 1.13

 $\displaystyle \overset{\tilde{\ }}{\mathop{K}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{b_{i}^{'}\kappa b_{j}^{'}d\Omega}$ Eq 1.14

Force matrix

 $\displaystyle \overset{\tilde{\ }}{\mathop{f}}\,(w^h)=-\int\limits_{\Omega }{w^hf(x,t)d\Omega }$ $\displaystyle =\sum\limits_{i}{{{c}_{i}}}\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }$ Eq 1.15

 $\displaystyle \overset{\tilde{\ }}{\mathop{F}}\,(b_i)=-\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }$ Eq 1.16

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

 $\mathbf{\tilde{M}}= \begin{bmatrix} \mathbf{M_{EE}} \quad\vline& \mathbf{M_{EF}}\\ \hline \mathbf{M_{FE}} \quad\vline& \mathbf{M_{FF}}\\ \end{bmatrix}$

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

 $\mathbf{\tilde{\kappa}}= \begin{bmatrix} \mathbf{\kappa_{EE}} \quad\vline& \mathbf{\kappa_{EF}}\\ \hline \mathbf{\kappa_{FE}} \quad\vline& \mathbf{\kappa_{FF}}\\ \end{bmatrix}$

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

 $\mathbf{\tilde{F}}= \begin{bmatrix} \mathbf{F_{EE}} \\ \hline \mathbf{F_{FF}} \\ \end{bmatrix}$

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

 $\displaystyle \mathbf{M_{FF}}\mathbf\dot{{d_{F}}}+ \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{{d_{F}}}=0$ and $\displaystyle \mathbf{g}$ = constant. Thus $\displaystyle \mathbf\dot{g} = 0.$

Taking this considerations, Eq 1.17 reduces to :-

 $\displaystyle \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 :-

 $\displaystyle \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:

Analysis of static case for heat problem f=1

Analysis of static case for heat problem f=0

### Solution for dynamic case

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

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

 $\displaystyle \mathbf{M_{FF}}\mathbf\dot{{d_{F}}}+ \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{{d_{F}}})$ $\displaystyle \mathbf{d_{F}}=\mathbf{K_{FF}^{-1}}(\mathbf{F_{F}}-(\mathbf{K_{FE}}\mathbf{g}+\mathbf{M_{FF}}\mathbf\dot{{d_{F}}}))$ (Eq 1.20 )

the initial conditions obtained is:

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

where

 $\displaystyle G={{\langle {{b}_{i}},{{u}_{0}}\rangle }_{M}}=\int\limits_{\Omega }{{{b}_{i}}}\rho c{{u}_{o}}d\Omega$ $\displaystyle \mathbf{u}=\mathbf{xy}$ (Eq 1.21)

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

Analysis of Dynamics case of heat problem(f=0)
Analysis of Dynamics case of heat problem(f=1)

The Matlab code to generate the plots and LIBF is:

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

## Given: Problem defined on Mtg 40 (c)

Refer to Lecture notes Lecture 40-5 Lecture 41-1,41-2 for more detailed description.

Strong Form

 $T div (grad u) + f = \rho \ \frac{\partial^2u}{\partial t^2} \qquad\text{on}\quad \Omega=\square$ (Eq 2.1)

Essential boundary conditions,

 $\displaystyle g = 0 \qquad \text {on} \quad {\partial \Omega }$ (Eq 2.2)

## 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

## Solution

### Constructing Weak Form and Discrete Weak Form

#### i. Weak Form

 $T div (grad u) + f = \rho \ \frac{\partial^2u}{\partial t^2} \qquad\text{on}\quad \Omega=\square$ (Eq 2.1)

From the PDE the weak form can be written as:

 $\displaystyle \int{w\left( T\frac{\partial }{\partial {{x}_{i}}}\left( \frac{\partial u}{\partial {{x}_{i}}} \right)+f\left( x,t \right)-\rho \frac{\partial^2u}{\partial t^2} \right)}d\text{ }\!\!\Omega\!\!\text{ =0}$ (Eq 2.3)

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

 $\displaystyle \int{w\left( T\frac{\partial }{\partial {{x}_{i}}}\left( \frac{\partial u}{\partial {{x}_{i}}} \right) \right)}d\text{ }\!\!\Omega\!\!\text{ =}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial }{\partial {{x}_{i}}}\left( wT\frac{\partial u}{\partial {{x}_{i}}} \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}T\frac{\partial u}{\partial {{x}_{i}}}d}\text{ }\!\!\Omega\!\!\text{ }$ (Eq 2.4)

Hence Eq. 2.1 can be now written as

 $\displaystyle \int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial }{\partial {{x}_{i}}}\left( wT\frac{\partial u}{\partial {{x}_{i}}} \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}T\frac{\partial u}{\partial {{x}_{i}}}d}\text{ }\!\!\Omega\!\!\text{ } { +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho \frac{\partial^2u}{\partial t^2}}d\text{ }\!\!\Omega\!\!\text{ }=0$ (Eq 2.5)

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,

 $\displaystyle \int\limits_{\text{ }\!\!\Omega\!\!\text{ }}wT\frac{\partial u}{\partial{{x}_{i}}}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}T\frac{\partial u}{\partial {{x}_{i}}}d}\text{ }\!\!\Omega\!\!\text{ } { +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho \frac{\partial^2u}{\partial t^2}}d\text{ }\!\!\Omega\!\!\text{ }=0$ (Eq 2.6)

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

 $\displaystyle \int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho \frac{\partial^2u}{\partial t^2}}d\text{ }\!\!\Omega\!\!\text{ }+\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}T\frac{\partial u}{\partial {{x}_{i}}}d}\text{ }\!\!\Omega\!\!\text{ }=-\int\limits_{{{\Gamma }_{h}}}wT\frac{\partial u}{\partial{{x}_{i}}}{{n}_{i}}d{{\Gamma }_{h}}-\int\limits_{\Omega }{wf(x,t)}d\Omega$ (Eq 2.7)

This Equation is of the form:-

 $\displaystyle \tilde{m}(w,u)+\tilde{k}(w,u)=\tilde{f}(w) \quad \forall w$ such that $\displaystyle w=0 \quad on \quad {{\Gamma }_{g}}$ (Eq 2.8)

where,

 $\displaystyle \tilde{m}(w,u)=\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho \frac{\partial^2u}{\partial t^2}d\text{ }\!\!\Omega\!\!\text{ }}$ $\displaystyle\tilde{k}(w,u)=\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}T\frac{\partial u}{\partial {{x}_{i}}}d}\text{ }\!\!\Omega\!\!\text{ }$ $\displaystyle\tilde{f}(w)=-\int\limits_{\Omega }{wf(x,t)}d\Omega$ (Eq 2.9)

#### ii. Finding approximate solution

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

 \displaystyle \begin{align} & u^{h}(x)=\sum_{j} d_{j} b_{j} \\ &\\ & w^{h}(x)=\sum_{i} c_{i} b_{i} \\ \end{align} (Eq 2.10)

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

 \displaystyle \begin{align} &N_{I}(x,y)=L_{i,m}(x)L_{j,n}(y)\\ & where \quad I=i+(j-1)m\\ \end{align}

Eq 2.11

$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

 $L_{i,m}(x) = \prod_{k = 1}^{m}\frac{x-x_{k}}{x_{i}-x_{k}}\;\; k\neq i$; where $\displaystyle m$ is number of nodes.

Eq 2.12

 $L_{j,n}(y) = \prod_{k = 1}^{m}\frac{y-y_{k}}{y_{j}-y_{k}}\;\; k\neq j$; where $\displaystyle n$ is number of nodes.

Eq 2.13

 $\sum_{i} c_i \left [ \sum_{j} \left \{ \tilde{M_{ij}} \ddot\tilde{d_j} \right \}+ \left \{ \tilde{K_{ij}} \tilde{d_j} \right \}- \tilde{F_i} \right ] = 0$ Eq 2.14

where,

Capacitance matrix

 $\displaystyle \overset{\tilde{\ }}{\mathop{m}}\,(w^h,u^h)=\int\limits_{\Omega }{w^h\rho \frac{\partial^2 u^h}{\partial t^2}d\Omega }$ $\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}\left( \int\limits_{\Omega }{{{b}_{i}}\rho {{b}_{j}}}d\Omega \right)}}\ddot{d_j}$ Eq 2.15

 $\displaystyle \overset{\tilde{\ }}{\mathop{M}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{{{b}_{i}}\rho {{b}_{j}}}d\Omega$ Eq 2.16

Conductivity matrix

 $\displaystyle k(w^h,u^h)=\int\limits_{\Omega }{\nabla w^h.T .\nabla u^hd\Omega}$ $\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}}\int\limits_{\Omega }{\left( \frac{\partial {{b}_{i}}}{\partial {{x}_{i}}}T \frac{\partial {{b}_{j}}}{\partial {{x}_{j}}}d\Omega \right)}}{{d}_{j}}$ Eq 2.17

 $\displaystyle \overset{\tilde{\ }}{\mathop{K}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{b_{i}^{'}T b_{j}^{'}d\Omega}$ Eq 2.18

Force matrix

 $\displaystyle \overset{\tilde{\ }}{\mathop{f}}\,(w^h)=-\int\limits_{\Omega }{w^hf(x,t)d\Omega }$ $\displaystyle =\sum\limits_{i}{{{c}_{i}}}\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }$ Eq 2.19

 $\displaystyle \overset{\tilde{\ }}{\mathop{F}}\,(b_i)=-\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }$ Eq 2.20

### 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}$

 $\sum_{i} c_i \left [ \sum_{j} \left \{ \tilde{K_{ij}} \tilde{d_j} \right \}- \tilde{F_i} \right ] = 0$ Eq 2.21

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

 $4 \left ( \frac{\partial^2u}{\partial x_{1}^2} + \frac{\partial^2u}{\partial x_{2}^2}\right ) + 1 = 0$ Eq 2.22

Membrane Analysis using 3x3 nodes
Membrane Analysis using 4x4 nodes
Membrane Analysis using 5x5 nodes
Membrane Analysis using 6x6 nodes
Membrane Analysis using 7x7 nodes
Convergence of Membrane Analysis

### 2. Dynamic solution

$\quad \rho \ = \ 3$

#### 2a. Free Vibration Analysis

The equation is,

 $\mathbf{K} \mathbf{\phi} = \lambda \mathbf{M}\mathbf{\phi}$ Eq 2.23
 $\mathbf{\phi}\left[ \mathbf{K}-\lambda \mathbf{M} \right]=0$ Eq 2.24
 $\lambda = \mathbf{M^{-1}} \mathbf{K}$ Eq 2.25

Where,
$\quad \lambda - Eigenvalues$

$\quad \mathbf{\phi} - Eigenvectors$

We calculate eigenvalues and then put them into the equation below to find frequencies.

The fundamental frequency is $\omega _{1}=\lambda_{1}^{1/2}$, and $f_{1}=\frac{\omega _{1}}{2\Pi }$.

The period T is $T_{1}=\frac{1}{f_{1}}$

#### 2b. Plotting displacement and the movie of the vibrating membrane for symmetric

 $T = 4, \qquad f(\mathbf{x},t) = 0 \quad \ in \ \Omega$ Eq 2.26

Initial conditions:

 $u(\mathbf{x},t=0) = u^h_s(x) \quad \ \forall \ \mathbf{x} \ \epsilon \ \Omega$ Eq 2.27
 $\dot{u} (\mathbf{x},t=0) = 0 \ \quad \forall \ \mathbf{x} \ \epsilon \ \Omega$ Eq 2.28

Note: To see the movie click on the picture twice and wait for a second for movie to start.

Analysis of Membrane problem with symmetric initial condition

#### 2c. Plotting displacement and the movie of the vibrating membrane for non-symmetric

 $T = 4, \qquad f(\mathbf{x},t) = 0 \quad \ in \ \Omega$ Eq 2.29

Initial conditions:

 $u(\mathbf{x},t=0) = (x+1)(y+\frac{1}{2}) cos (\frac{\pi}{2}x) cos (\frac{\pi}{2}y) \quad \ \forall \ \mathbf{x} \ \epsilon \ \Omega$ Eq 2.30
 $\dot{u} (\mathbf{x},t=0) = 0 \ \quad \forall \ \mathbf{x} \ \epsilon \ \Omega$ Eq 2.31

Note: To see the movie click on the picture twice and wait for a second for movie to start.

Analysis of Membrane problem with non-symmetric initial condition

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

Node numbers for the quadrilateral case
Node numbers for the triangular case
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 had to terminate the program.

# Problem 7.4

## Problem statement

Please refer to lecture note 37-2 and 38 for more details.

Solve the heat conduction problem in 2D with boundary convection.

1> Construct the weak form for heat conduction in 2D with boundary convection.

2> Develop semi-discrete equations (ODEs) similar to (3) 23-3 Give detailed expression of all quantities.

3>Solve G2DM2.0/D1 using 2D LIBF till $\displaystyle 10^{-6}$ accuracy at the center. Plot solution in 3D with contour.

Problem figure

## Solution

### Developing the Weak Form of Newton's law of Cooling

i> From Lecture 33-2, the PDE of the problem is given by:-

 $\displaystyle \frac{\partial }{{\partial{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{{\partial{x}_{j}}} \right)+f\left( x,t \right)=\rho c\frac{\partial u}{\partial t}$ (Eq )<1>

From the PDE the weak form can be written as:-

 $\displaystyle \int{w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)+f\left( x,t \right)-\rho c\frac{\partial u}{\partial t} \right)}d\text{ }\!\!\Omega\!\!\text{ =0}$ (Eq )<2>

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

 $\displaystyle \int{w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right) \right)}d\text{ }\!\!\Omega\!\!\text{ =}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial }{\partial {{x}_{i}}}\left( w{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}}d}\text{ }\!\!\Omega\!\!\text{ }$

Hence equation 2 can be now written as:-

 $\displaystyle \int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial }{\partial {{x}_{i}}}\left( w{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0$ (Eq )<3>

Now we can write $\displaystyle d\Omega$ as $\displaystyle d\Omega ={{\Gamma }_{g}}\cup {{\Gamma }_{h}}\cup {{\Gamma }_{H}}$, where
$\displaystyle {{\Gamma }_{g}}$ = essential boundary condition.
$\displaystyle {{\Gamma }_{h}}$ = natural boundary condition.
$\displaystyle {{\Gamma }_{H}}$ = Mixed boundary condition.

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

 $\displaystyle \int\limits_{\text{ }\!\!\Omega\!\!\text{ }}w{{K}_{ij}}\frac{\partial u}{\partial{{x}_{j}}}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0$ (Eq )<4>

Now the heat flux is given by:-

 $\displaystyle {{q}_{i}}=-{{\text{ }\!\!\Kappa\!\!\text{ }}_{ij}}\frac{\partial u}{\partial {{x}_{j}}}$ (Eq )<5>

So by rearranging Equation 5, we can write

 $\displaystyle \int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{-w{{q}_{i}}{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0$ (Eq )<6>

We can write the first term as-

$\displaystyle \int\limits_{\Gamma }{-w{{q}_{i}}{{n}_{i}}d\Gamma =}\int\limits_{{{\Gamma }_{g}}}{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{g}}}-\int\limits_{{{\Gamma }_{h}}}{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{h}}}-\int\limits_{{{\Gamma }_{H}}}{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{H}}}$

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

On $\displaystyle {{\Gamma }_{h}},\quad q(x,t)n(x)=h(x,t) ,\forall x\in {{\Gamma }_{h}}$

$\displaystyle {{\Gamma }_{H}},\quad q(x,t)n(x)=H[u(x,t)-{u}_{\infty}] ,\forall x\in {{\Gamma }_{H}}$

So we can write Equation 6 as :-

 $\displaystyle -\int\limits_{{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}{whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}-\int\limits_{{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}{wH(u-{{u}_{\infty }})d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)d\text{ }\!\!\Omega\!\!\text{ }}+\int\limits_{\Omega }{wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-}\int\limits_{\Omega }{w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=0}$ (Eq )<7>

This can be re-arranged in the form:-

 $\displaystyle \int\limits_{\Omega }{w\rho c\frac{\partial u}{\partial t}d\Omega +\int\limits_{\Omega }{\nabla w.\kappa .\nabla ud\Omega +}}\int\limits_{{{\Gamma }_{H}}}{wHu}d{{\Gamma }_{H}}=-\int\limits_{{{\Gamma }_{h}}}{wh}d{{\Gamma }_{h}}+\int\limits_{{{\Gamma }_{H}}}{wH{{u}_{\infty }}}d{{\Gamma }_{H}}+\int\limits_{\Omega }{wf(x,t)}d\Omega$ (Eq )<8>

This Equation is of the form:-

 $\displaystyle \tilde{m}(w,u)+\tilde{k}(w,u)=\tilde{f}(w) \quad \forall w$ such that $\displaystyle w=0 \quad on \quad {{\Gamma }_{g}}$

where,

 $\displaystyle \tilde{m}(w,u)=\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }}$ $\displaystyle\tilde{k}(w,u)=\int\limits_{\Omega }{\nabla w.\kappa .\nabla ud\Omega }+\int\limits_{{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}{wHu}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$ $\displaystyle\tilde{f}(w)=-\int\limits_{{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}{whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}+\int\limits_{{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}{wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}+\int\limits_{\Omega }{wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }}$ (Eq )<9>

### Developing Semi-Discrete equations

ii>

To form the semi-discrete equations, we take $\displaystyle w \quad and \quad u$ as:-

 $\displaystyle w^h=\sum\limits_{i}{{{c}_{i}}{{b}_{i}}} \quad and \quad u^h=\sum\limits_{j}{{{d}_{j}}{{b}_{j}}}$ (Eq )<10>

Accordingly the mass, capacitance and heat source equations can be written as:-

Capacitance matrix

 $\displaystyle \overset{\tilde{\ }}{\mathop{m}}\,(w^h,u^h)=\int\limits_{\Omega }{w^h\rho c\frac{\partial u^h}{\partial t}d\Omega }$ $\displaystyle=\int\limits_{\Omega }{\sum{{{c}_{i}}{{b}_{i}}\rho c\sum{{{d}_{j}}{{b}_{j}}d\Omega }}}$ $\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}\left( \int\limits_{\Omega }{{{b}_{i}}\rho c{{b}_{j}}}d\Omega \right)}}d_{j}^{'}$

 $\displaystyle \overset{\tilde{\ }}{\mathop{M}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{{{b}_{i}}\rho c{{b}_{j}}}d\Omega$ (Eq )<11>

Conductivity matrix

 $\displaystyle k(w^h,u^h)=\int\limits_{\Omega }{\nabla w^h.\kappa .\nabla u^hd\Omega +\int\limits_{{{\Gamma }_{H}}}{w^hHu^hd{{\Gamma }_{H}}}}$ $\displaystyle=\int\limits_{\Omega }{\frac{\partial \sum{{{c}_{i}}{{b}_{i}}}}{\partial {{x}_{i}}}}.\kappa .\frac{\partial \sum{{{d}_{j}}{{b}_{j}}}}{\partial {{x}_{j}}}d\Omega +\int\limits_{{{\Gamma }_{\Eta }}}{\sum{{{c}_{i}}{{b}_{i}}}H\sum{{{d}_{j}}{{b}_{j}}}d{{\Gamma }_{\Eta }}}$ $\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}}\int\limits_{\Omega }{\left( \frac{\partial {{b}_{i}}}{\partial {{x}_{i}}}.\kappa \frac{\partial {{b}_{j}}}{\partial {{x}_{j}}}d\Omega \right)}}{{d}_{j}}+\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}}\int\limits_{{{\Gamma }_{\Eta }}}{\left( {{b}_{i}}H{{b}_{j}}d{{\Gamma }_{\Eta }} \right)}}{{d}_{j}}$

 $\displaystyle \overset{\tilde{\ }}{\mathop{K}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{b_{i}^{'}\kappa b_{j}^{'}d\Omega +\int\limits_{{{\Gamma }_{\Eta }}}{{{b}_{i}}H{{b}_{j}}d}}{{\Gamma }_{\Eta }}$ (Eq )<12>

Force matrix

 $\displaystyle \overset{\tilde{\ }}{\mathop{f}}\,(w^h)=-\int\limits_{{{\Gamma }_{h}}}{w^hhd{{\Gamma }_{h}}}+\int\limits_{{{\Gamma }_{H}}}{w^hH{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{w^hf(x,t)d\Omega }$ $\displaystyle =-\int\limits_{{{\Gamma }_{h}}}{\sum\limits_{i}{{{c}_{i}}{{b}_{i}}}hd{{\Gamma }_{h}}}+\int\limits_{{{\Gamma }_{H}}}{\sum\limits_{i}{{{c}_{i}}{{b}_{i}}}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{\sum\limits_{i}{{{c}_{i}}{{b}_{i}}}f(x,t)d\Omega }$ $\displaystyle =-\sum\limits_{i}{{{c}_{i}}}\int\limits_{{{\Gamma }_{h}}}{{{b}_{i}}hd{{\Gamma }_{h}}}+\sum\limits_{i}{{{c}_{i}}}\int\limits_{{{\Gamma }_{H}}}{{{b}_{i}}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\sum\limits_{i}{{{c}_{i}}}\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }$

 $\displaystyle \overset{\tilde{\ }}{\mathop{F}}\,(b_i)=-\int\limits_{{{\Gamma }_{h}}}{{{b}_{i}}hd{{\Gamma }_{h}}}+\int\limits_{{{\Gamma }_{H}}}{{{b}_{i}}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }$ (Eq )<13>

### 3D Contour plots

iii>

consider the case of n=m=3 the number of nodes is 9 The trial solution is

 $\displaystyle u_{p}^{h}(x)={{d}_{1}}*N_{1}+{{d}_{2}}N_{2}+{{d}_{3}}N_{3}+\cdots +{{d}_{9}}N_{9},_{{}}^{{}}j=0,1,2,\cdots ,9$

Since the trial solution must satisfy the essential boundary condition,, we should have ${{d}_{1}}=2,{{d}_{2}}=2,{{d}_{3}}=2,{{d}_{6}}=2,{{d}_{9}}=2$., because it is specified that through out the boundary the temperature is 2 i.e for all the nodes on the boundary temperature is 2

In order to solve for the rest coefficient ${{d}_{j}},_{{}}^{{}}j=1,2,\cdots$, we construct the Conductance matrix, Heat Source matrix and Capacitance matrix as given by Equations 12,13 and 14.

Using Matlab to construct above matrices until satisfying the convergence criterion, we finally have

 $\displaystyle {{K}}=\left( \begin{matrix} 151/200 & -1/5 & -1/30 & -133/1000 & -16/45 & 1/9 & -33/500 & 1/9 & -1/45\\ -1/5 & 88/45 & -1/5 & -16/45 & -16/15 & -16/45 & 1/9 & 0 & 1/9\\ -1/30 & -1/5 & 28/45 & 1/9 & -16/45 & -1/5 & -1/45 & 1/9 & -1/30\\ -133/100 & -16/45 & 1/9 & 112/45 & -16/15 & 0 & -133/1000 & -16/45 & 1/9\\ -16/45 & -16/15 & -16/45 & -16/15 & 256/45 & -16/15 & -16/45 & -16/15 & -16/45\\ 1/9 & -16/45 & -1/5 & 0 & -16/15 & 88/45 & 1/9 & -16/45 & -1/5\\ -33/500 & 1/9 & -1/45 & -133/1000 & -16/45 & 1/9 & 1511/2000 & -1/5 & -1/30\\ 1/9 & 0 & 1/9 & -16/45 & -16/15 & -16/45 & -1/5 & 88/45 & -1/5\\ -1/45 & 1/9 & -1/30 & 1/9 & -16/45 & -1/5 & -1/30 & -1/5 & 28/45\\ \end{matrix} \right)$

 $\displaystyle {F}=\left( \begin{matrix} 166/1000\\ 0\\ 0\\ 33/500\\ 0\\ 0\\ -59/60\\ -4\\ -1\\ \end{matrix} \right)$

 $\displaystyle {d}=\left( \begin{matrix} 2.000\\ 2.000\\ 2.000\\ 0.335\\ 0.847\\ 2.000\\ -1.460\\ -1.331\\ 2.000\\ \end{matrix} \right)$

The trial solution using 2D LIBF obtained is as follows:-

 \displaystyle \begin{align} &U_{F}^{h}=(x)0.50*x*y*(x - 1)*(y - 1) - 1*x*(1.0*x + 1.0)*(1.0*y + 1.0)*(y - 1)\\& \quad + 0.665*y*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1) - 0.1678*x*(1.0*y + 1.0)*(x - 1)*(y - 1)\\& \quad - 1.000*y*(1.0*x + 1.0)*(x - 1)*(y - 1) + 0.847*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1)*(y - 1) \\& \quad + 0.500*x*y*(1.0*x + 1.0)*(1.0*y + 1.0) + 0.500*x*y*(1.0*x + 1.0)*(y - 1) - 0.365*x*y*(1.0*y + 1.0)*(x - 1)\\ \end{align} (Eq )<14>

The Matlab code to generate the plots and LIBF is:

The 3D contour plot generated from MATLAB is dislpayed below.

3D plot for n=m=3

# Contributing Members & Referenced Lecture

 Team 3 Contribution Table Problem Number Lecture Assigned To Solved By Typed By Proofread By 7.1 Lecture 39-1 Lecture 40-4 Avinash Avinash Avinash ushnish Hailong Chen 7.2 Lecture 40-5 Lecture 41-1,41-2 Hailong Chen sahin Hailong Chen sahin Hailong Chen sahin Avinash 7.3 Lecture 41-3 vnarayanan risher risher vnarayanan vnarayanan User:Eml5526.s11.team3 7.4 Lecture 41-4 ushnish ushnish ushnish User:Eml5526.s11.team3