User:Eml5526.s11.team3/Homework 7

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:

MATLAB Code

Matlab Code: %----------------------Problem 7.1 Static Case Heat Problem for f=1----------------------------- clear all close all clc syms x y; n=4;%the same number of nodes in x and y direction % T=4;%tension constant f=1;%force resource L_x=sym('L_x');%preallocate space for the functions interpolated in the x direction L_y=sym('L_y');%preallocate space for the functions interpolated in the y direction N=sym('N');%preallocate space for the shape functions xcoord=(n);%xcoordinate of different nodes ycoord=(n);%ycoordinate of different nodes d=zeros(n*n,1);%preallocate space for the displacement vector flag=ones(n*n,1);%preallocate space for flag K=[n*n,n*n];%preallocate space for the stiffness matrix F=[n*n,1];%preallocate space for the force vector e_nd=4*n-4;% number of nodes on the essential boundary %get the x- and y-coordinate of the nodes for i=1:n xcoord(i)=-1+(i-1)*2/(n-1); ycoord(i)=-1+(i-1)*2/(n-1); end %obtain the lagrange interpolation basis function in both directions %------------------------------------------------------------------- for i=1:n %LIBF in x direction L_x(i)=1; for j=1:n if i~=j L_x(i)=L_x(i)*(x-xcoord(j))/(xcoord(i)-xcoord(j)); end end %LIBF in y direction L_y(i)=1; for j=1:n if i~=j L_y(i)=L_y(i)*(y-ycoord(j))/(ycoord(i)-ycoord(j)); end end end %nodal shape functions, begin from the left-bottom node to right-top node step=0; for i=1:n for j=1:n N(j+step)=simple(L_y(i)*L_x(j)); end step=step+n;%increase row by row end %-------------------------------------------------------------------- %essential boundary condition %-------------------------------------------------------------------- %obtain the indices of free nodes freeindex=(n*n-e_nd); step1=0; step2=0; for i=1:(n-2) freeindex(1+step1:n-2+step1)=n+2+step2:2*n-1+step2; step1=step1+n-2; step2=step2+n; end %flag the nodes on the free boundary for i=1:length(freeindex) flag(freeindex(i))=0; end ID=(n*n); count = 0; count1 = 0; for i = 1:n*n if flag(i)==1% check if essential B.C count=count + 1; ID(i)=count;% arrange essential B.C nodes first d(count)=0;% store reordered essential B.C else count1=count1+1; ID(i)=e_nd+count1; end end temp_N=N; %rearrange the shape functions for i=1:n*n N(ID(i))= temp_N(i); end %------------------------------------------------------------------- %stiffness and force matrix %------------------------------------------------------------------- for i=1:n*n for j=1:n*n % K(i,j)=int(int(diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'),'x',-1,1),'y',-1,1); K1=diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'); K1_fn=matlabFunction(K1); K(i,j)=dblquad(K1_fn,-1,1,-1,1); end F(i,1)=int(int(f*N(i),'x',-1,1),'y',-1,1); % F1=f*N(i); % F1_fn=matlabFunction(F1); % F(i)=dblquad(F1_fn,-1,1,-1,1); end %------------------------------------------------------------------- %extract the stiffness matrix and force matrix of free nodes %------------------------------------------------------------------- K_F=K(e_nd+1:n*n,e_nd+1:n*n); F_F=F(e_nd+1:n*n); K_FE=K(e_nd+1:n*n,1:e_nd); d_E=2+d(1:e_nd,1); %------------------------------------------------------------------ %solve for d_F d_F=K_F\(F_F'-K_FE*d_E); %------------------------------------------------------------------ %d matrix for i=1:n*n if i<=e_nd d(i)=d_E(i); else d(i)=d_F(i-e_nd); end end %----------------------------------------------------------------- %trial solution Uh=N*d; %----------------------------------------------------------------- %surf the figure [X,Y]=meshgrid(-1:0.1:1,-1:0.1:1); U=subs(Uh,{x,y},{X,Y}); surf(X,Y,U); colorbar; title('Analysis of Static Case of Heat problem f=1'); xlabel('x'); ylabel('y'); % zlim([1 3]) zlabel('Temperature');

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:

MATLAB Code

Matlab Code: %----------------------Problem 7.1 Dynamic Case of heat problem for f=0----------------------------- clear all close all clc syms x y; global M_F F_FF K_F n n=4;%the same number of nodes in x and y direction rhoc=3; f=0;%force resource L_x=sym('L_x');%preallocate space for the functions interpolated in the x direction L_y=sym('L_y');%preallocate space for the functions interpolated in the y direction N=sym('N');%preallocate space for the shape functions xcoord=(n);%xcoordinate of different nodes ycoord=(n);%ycoordinate of different nodes d=zeros(n*n,1);%preallocate space for the displacement vector flag=ones(n*n,1);%preallocate space for flag K=[n*n,n*n];%preallocate space for the stiffness matrix F=[n*n,1];%preallocate space for the force vector e_nd=4*n-4;% number of nodes on the essential boundary %get the x- and y-coordinate of the nodes for i=1:n xcoord(i)=-1+(i-1)*2/(n-1); ycoord(i)=-1+(i-1)*2/(n-1); end Y_initial=-1; for j=1:n for i=1:n I=i+(j-1)*n; x_T(I)=xcoord(i); y_T(I)=Y_initial; end Y_initial=Y_initial+2/(n-1); end %obtain the lagrange interpolation basis function in both directions %------------------------------------------------------------------- for i=1:n %LIBF in x direction L_x(i)=1; for j=1:n if i~=j L_x(i)=L_x(i)*(x-xcoord(j))/(xcoord(i)-xcoord(j)); end end %LIBF in y direction L_y(i)=1; for j=1:n if i~=j L_y(i)=L_y(i)*(y-ycoord(j))/(ycoord(i)-ycoord(j)); end end end %nodal shape functions, begin from the left-bottom node to right-top node step=0; for i=1:n for j=1:n N(j+step)=simple(L_y(i)*L_x(j)); end step=step+n;%increase row by row end %-------------------------------------------------------------------- %essential boundary condition %-------------------------------------------------------------------- %obtain the indices of free nodes freeindex=(n*n-e_nd); step1=0; step2=0; for i=1:(n-2) freeindex(1+step1:n-2+step1)=n+2+step2:2*n-1+step2; step1=step1+n-2; step2=step2+n; end %flag the nodes on the free boundary for i=1:length(freeindex) flag(freeindex(i))=0; end ID=(n*n); count = 0; count1 = 0; for i = 1:n*n if flag(i)==1% check if essential B.C count=count + 1; ID(i)=count;% arrange essential B.C nodes first d(count)=0;% store reordered essential B.C else count1=count1+1; ID(i)=e_nd+count1; end end temp_N=N; temp_x=x_T; temp_y=y_T; %rearrange the shape functions for i=1:n*n N(ID(i))= temp_N(i); x_T(ID(i))=temp_x(i); y_T(ID(i))=temp_y(i); end %------------------------------------------------------------------- %stiffness and force matrix %------------------------------------------------------------------- for i=1:n*n for j=1:n*n K1=diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'); K1_fn=matlabFunction(K1); K(i,j)=dblquad(K1_fn,-1,1,-1,1); M1=rhoc*N(i)*N(j); M1_fn=matlabFunction(M1); M(i,j)=dblquad(M1_fn,-1,1,-1,1); end F(i,1)=int(int(f*N(i),'x',-1,1),'y',-1,1); % F1=f*N(i); % F1_fn=matlabFunction(F1); % F(i)=dblquad(F1_fn,-1,1,-1,1); end %------------------------------------------------------------------- %extract the stiffness matrix and force matrix of free nodes %------------------------------------------------------------------- K_F=K(e_nd+1:n*n,e_nd+1:n*n); % Calulatiing the K matrix for free degrees of freedom M_F=M(e_nd+1:n*n,e_nd+1:n*n); % Calulatiing the M matrix for free degrees of freedom F_F=F(e_nd+1:n*n); % Calulatiing the F matrix for free degrees of freedom K_FE=K(e_nd+1:n*n,1:e_nd); M_FE=M(e_nd+1:n*n,1:e_nd); d_E=2+d(1:e_nd,1); % Given temperature on boundary(essential BC) N_F=N(e_nd+1:n*n); % Basis functions for free nodes x_F=x_T(e_nd+1:n*n); % x coordinates for free nodes y_F=y_T(e_nd+1:n*n); % y coordinates for free nodes %------------------------------------------------------------------ % Calculating G matrix to find the temperatures at initial condition %------------------------------------------------------------------ for i=1:n*n-(4*n-4) G1(i)=N_F(i)*rhoc*x_F(i)*y_F(i); G1_fn=matlabFunction(G1(i)); G(i)=dblquad(G1_fn,-1,1,-1,1); end %------------------------------------------------------------------ % Initial conditions do=M_F\G'; F_FF=F_F'-K_FE*d_E; %------------------------------------------------------------------ %solve for d_F % Using ODE45 to solve for the first order time derivative equation in % dynamic heat case [T Y]=ode45(@transientheat,[0:1:20],do); %------------------------------------------------------------------ %d matrix for j=1:length(Y) for i=1:n*n if i<=e_nd d(i)=d_E(i); else d(i)=Y(j,i-e_nd)'; end end %----------------------------------------------------------------- %trial solution Uh=N*d; %----------------------------------------------------------------- %surf the figure figure(1) [X,Z]=meshgrid(-1:0.1:1,-1:0.1:1); U=subs(Uh,{x,y},{X,Z}); surf(X,Z,U); colorbar; title('Analysis of Dynamic Case of Heat Problem f=0 at Initial Condition'); xlabel('x'); ylabel('y'); zlim([0 2.5]); zlabel('Temperature'); % view(0,0);%change the view angle hold off end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dd_F = transientheat(t,y) global M_F F_FF K_F n % d_F=zeros(n,1); dd_F=zeros(n,1); dd_F = M_F\(F_FF-K_F*y); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

Matlab Code for Vibrating Membrane (Static case)

Static Case for Vibrating Membrane using 2D LIBF

%----------------------Problem 7.2 Static Case----------------------------- % Coded by Hailong Chen, FEM1.S11.team3, 4/16/2011 %-------------------------------------------------------------------------- clear all close all clc syms x y; format longeng; n=3;%the same number of nodes in x and y direction T=4;%tension constant f=1;%force resource L_x=sym('L_x',[1 1]);%preallocate space for the functions interpolated in the x direction L_y=sym('L_y',[1 1]);%preallocate space for the functions interpolated in the y direction N=sym('N',(1));%preallocate space for the shape functions xcoord=(n);%xcoordinate of different nodes ycoord=(n);%ycoordinate of different nodes d=zeros(n*n,1);%preallocate space for the displacement vector flag=ones(n*n,1);%preallocate space for flag K=zeros(n*n,n*n);%preallocate space for the stiffness matrix F=zeros(n*n,1);%preallocate space for the force vector e_nd=4*n-4;% number of nodes on the essential boundary %get the x- and y-coordinate of the nodes for i=1:n xcoord(i)=-1+(i-1)*2/(n-1); ycoord(i)=-1+(i-1)*2/(n-1); end %obtain the lagrange interpolation basis function in both directions %------------------------------------------------------------------- for i=1:n %LIBF in x direction L_x(i)=1; for j=1:n if i~=j L_x(i)=L_x(i)*(x-xcoord(j))/(xcoord(i)-xcoord(j)); end end %LIBF in y direction L_y(i)=1; for j=1:n if i~=j L_y(i)=L_y(i)*(y-ycoord(j))/(ycoord(i)-ycoord(j)); end end end %nodal shape functions, begin from the left-bottom node to right-top node step=0; for i=1:n for j=1:n N(j+step)=L_y(i)*L_x(j); end step=step+n;%increase row by row end %-------------------------------------------------------------------- %essential boundary condition %-------------------------------------------------------------------- %obtain the indices of free nodes freeindex=(n*n-e_nd); step1=0; step2=0; for i=1:(n-2) freeindex(1+step1:n-2+step1)=n+2+step2:2*n-1+step2; step1=step1+n-2; step2=step2+n; end %flag the nodes on the free boundary for i=1:length(freeindex) flag(freeindex(i))=0; end ID=(n*n); count = 0; count1 = 0; for i = 1:n*n if flag(i)==1% check if essential B.C count=count + 1; ID(i)=count;% arrange essential B.C nodes first d(count)=0;% store reordered essential B.C else count1=count1+1; ID(i)=e_nd+count1; end end temp_N=N; %rearrange the shape functions for i=1:n*n N(ID(i))= temp_N(i); end %------------------------------------------------------------------- %stiffness and force matrix %------------------------------------------------------------------- for i=1:n*n for j=1:n*n K(i,j)=T*int(int(diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'),'x',-1,1),'y',-1,1); end F(i,1)=int(int(f*N(i),'x',-1,1),'y',-1,1); end %------------------------------------------------------------------- %extract the stiffness matrix and force matrix of free nodes %------------------------------------------------------------------- K_F=K(e_nd+1:n*n,e_nd+1:n*n); F_F=F(e_nd+1:n*n,1); K_FE=K(e_nd+1:n*n,1:e_nd); d_E=d(1:e_nd,1); %------------------------------------------------------------------ %solve for d_F d_F=vpa(K_F)\vpa(F_F-K_FE*d_E); %------------------------------------------------------------------ %d matrix for i=1:n*n if i<=e_nd d(i)=d_E(i); else d(i)=d_F(i-e_nd); end end %----------------------------------------------------------------- %trial solution Uh=N*d; %----------------------------------------------------------------- %surf the figure [X,Y]=meshgrid(-1:0.1:1,-1:0.1:1); U=subs(Uh,{x,y},{X,Y}); surf(X,Y,U); colorbar; title('Analysis of Static Case of Membrane (9x9)'); xlabel('x'); ylabel('y'); zlabel('displacement'); center_approx=subs(Uh,{x,y},{0,0});%calculate the displacement of trial solution at the center %text(0,0,center_approx+0.01,'(0,0,0.07367127770)');%label the crest of the figure %view(0,0);%change the view angle %---------------------------------------------------------------- %value of center point for different interpolation nodes %n=3x3 0.07812500000 error between different n %n=4x4 0.07812500000 0 %n=5x5 0.07373046875 0.00439453125 %n=6x6 0.07373046875 0 %n=7x7 0.07367173855 0.00005873020

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

Dynamic Case for Free Vibration using 2D LIBF

%----------------------Problem 7.2a Static Case----------------------------- clear all close all clc d_Fs=hw7_2_static(4); syms x y; global M_F F_FF K_F n n=4;%the same number of nodes in x and y direction rho=3; T=4;%tension constant f=0;%force resource L_x=sym('L_x');%preallocate space for the functions interpolated in the x direction L_y=sym('L_y');%preallocate space for the functions interpolated in the y direction N=sym('N');%preallocate space for the shape functions xcoord=(n);%xcoordinate of different nodes ycoord=(n);%ycoordinate of different nodes d=zeros(n*n,1);%preallocate space for the displacement vector flag=ones(n*n,1);%preallocate space for flag K=[n*n,n*n];%preallocate space for the stiffness matrix F=[n*n,1];%preallocate space for the force vector e_nd=4*n-4;% number of nodes on the essential boundary %get the x- and y-coordinate of the nodes for i=1:n xcoord(i)=-1+(i-1)*2/(n-1); ycoord(i)=-1+(i-1)*2/(n-1); end Y_initial=-1; for j=1:n for i=1:n I=i+(j-1)*n; x_T(I)=xcoord(i); y_T(I)=Y_initial; end Y_initial=Y_initial+2/(n-1); end %obtain the lagrange interpolation basis function in both directions %------------------------------------------------------------------- for i=1:n %LIBF in x direction L_x(i)=1; for j=1:n if i~=j L_x(i)=L_x(i)*(x-xcoord(j))/(xcoord(i)-xcoord(j)); end end %LIBF in y direction L_y(i)=1; for j=1:n if i~=j L_y(i)=L_y(i)*(y-ycoord(j))/(ycoord(i)-ycoord(j)); end end end %nodal shape functions, begin from the left-bottom node to right-top node step=0; for i=1:n for j=1:n N(j+step)=simple(L_y(i)*L_x(j)); end step=step+n;%increase row by row end %-------------------------------------------------------------------- %essential boundary condition %-------------------------------------------------------------------- %obtain the indices of free nodes freeindex=(n*n-e_nd); step1=0; step2=0; for i=1:(n-2) freeindex(1+step1:n-2+step1)=n+2+step2:2*n-1+step2; step1=step1+n-2; step2=step2+n; end %flag the nodes on the free boundary for i=1:length(freeindex) flag(freeindex(i))=0; end ID=(n*n); count = 0; count1 = 0; for i = 1:n*n if flag(i)==1% check if essential B.C count=count + 1; ID(i)=count;% arrange essential B.C nodes first d(count)=0;% store reordered essential B.C else count1=count1+1; ID(i)=e_nd+count1; end end temp_N=N; temp_x=x_T; temp_y=y_T; %rearrange the shape functions for i=1:n*n N(ID(i))= temp_N(i); x_T(ID(i))=temp_x(i); y_T(ID(i))=temp_y(i); end %------------------------------------------------------------------- %stiffness and force matrix %------------------------------------------------------------------- for i=1:n*n for j=1:n*n % K(i,j)=T*int(int(diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'),'x',-1,1),'y',-1,1); % M(i,j)=int(int(rho*N(i)*N(j),'x',-1,1),'y',-1,1); K1=T*diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'); K1_fn=matlabFunction(K1); K(i,j)=dblquad(K1_fn,-1,1,-1,1); M1=rho*N(i)*N(j); M1_fn=matlabFunction(M1); M(i,j)=dblquad(M1_fn,-1,1,-1,1); end F(i)=int(int(f*N(i),'x',-1,1),'y',-1,1); % F1=f*N(i); % F1_fn=matlabFunction(F1); % F(i)=dblquad(F1_fn,-1,1,-1,1); end %------------------------------------------------------------------- %extract the stiffness matrix and force matrix of free nodes %------------------------------------------------------------------- K_F=K(e_nd+1:n*n,e_nd+1:n*n); M_F=M(e_nd+1:n*n,e_nd+1:n*n); F_F=F(e_nd+1:n*n); K_FE=K(e_nd+1:n*n,1:e_nd); M_FE=M(e_nd+1:n*n,1:e_nd); d_E=d(1:e_nd,1); N_F=N(e_nd+1:n*n); x_F=x_T(e_nd+1:n*n); y_F=y_T(e_nd+1:n*n); for i=1:n*n-(4*n-4) % G(i)=int(int(rho*N_F(i)*x_F(i)*y_F(i),'x',-1,1),'y',-1,1); G1(i)=N_F(i)*rho*x_F(i)*y_F(i); G1_fn=matlabFunction(G1(i)); G(i)=dblquad(G1_fn,-1,1,-1,1); end %------------------------------------------------------------------ %solve for d_F [PHI,LAMDA]=eig(K_F,M_F); omega1=sqrt(LAMDA(1,1)); T1= (2*pi)/omega1; do=d_Fs; ddo=zeros(n*n-(4*n-4),1); F_FF=F_F'-K_FE*d_E; doo=[do ;ddo]; %------------------------------------------------------------------ %solve for d_F [T Y]=ode45('dynamicmembrane',[0:0.3:2*T1],doo); %------------------------------------------------------------------ %d matrix for j=1:length(Y) for i=1:n*n if i<=e_nd d(i)=d_E(i); else d(i)=Y(j,i-e_nd); end end %----------------------------------------------------------------- %trial solution Uh=N*d; %----------------------------------------------------------------- %surf the figure [X,Z]=meshgrid(-1:0.1:1,-1:0.1:1); U=subs(Uh,{x,y},{X,Z}); surf(X,Z,U); colorbar; title('Analysis of Static Case of Membrane'); zlim([-2 2]); xlabel('x'); ylabel('y'); zlabel('displacement'); M(j) = imgcf; hold off end % for i=1:length(Y) % str = strcat('D:\Graduate Courses\EML 5526 FEM\HW\hw7.2b\',num2str(i), '.jpg'); % A=imread(str); % [I,map]=rgb2ind(A,256); % if(i==1) % imwrite(I,map,'fe1.s11.team3.hw7.2.b.gif','DelayTime',0.01,'LoopCount',Inf) % else % imwrite(I,map,'fe1.s11.team3.hw7.2.b.gif','WriteMode','append','DelayTime',0.01) % end % end %-------------------------------------------------------------------------------------------- function dy = dynamicmembrane(t,y) global M_F F_FF K_F n % d_F=zeros(n,1); dy=zeros(8,1); % y1(1)=y(1); % y2(2)=y(2); dy(1:4) = y(5:8); dy(5:8) = M_F\(F_FF-K_F*y(1:4)); end

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

Dynamic Case for Free Vibration using 2D LIBF

%----------------------Problem 7.2b Dynamic Case----------------------------- clear all close all clc % d_Fs=hw7_2_static(3); syms x y; global M_F F_FF K_F n n=4;%the same number of nodes in x and y direction rho=3; T=4;%tension constant f=0;%force resource L_x=sym('L_x');%preallocate space for the functions interpolated in the x direction L_y=sym('L_y');%preallocate space for the functions interpolated in the y direction N=sym('N');%preallocate space for the shape functions xcoord=(n);%xcoordinate of different nodes ycoord=(n);%ycoordinate of different nodes d=zeros(n*n,1);%preallocate space for the displacement vector flag=ones(n*n,1);%preallocate space for flag K=[n*n,n*n];%preallocate space for the stiffness matrix F=[n*n,1];%preallocate space for the force vector e_nd=4*n-4;% number of nodes on the essential boundary %get the x- and y-coordinate of the nodes for i=1:n xcoord(i)=-1+(i-1)*2/(n-1); ycoord(i)=-1+(i-1)*2/(n-1); end Y_initial=-1; for j=1:n for i=1:n I=i+(j-1)*n; x_T(I)=xcoord(i); y_T(I)=Y_initial; end Y_initial=Y_initial+2/(n-1); end %obtain the lagrange interpolation basis function in both directions %------------------------------------------------------------------- for i=1:n %LIBF in x direction L_x(i)=1; for j=1:n if i~=j L_x(i)=L_x(i)*(x-xcoord(j))/(xcoord(i)-xcoord(j)); end end %LIBF in y direction L_y(i)=1; for j=1:n if i~=j L_y(i)=L_y(i)*(y-ycoord(j))/(ycoord(i)-ycoord(j)); end end end %nodal shape functions, begin from the left-bottom node to right-top node step=0; for i=1:n for j=1:n N(j+step)=simple(L_y(i)*L_x(j)); end step=step+n;%increase row by row end %-------------------------------------------------------------------- %essential boundary condition %-------------------------------------------------------------------- %obtain the indices of free nodes freeindex=(n*n-e_nd); step1=0; step2=0; for i=1:(n-2) freeindex(1+step1:n-2+step1)=n+2+step2:2*n-1+step2; step1=step1+n-2; step2=step2+n; end %flag the nodes on the free boundary for i=1:length(freeindex) flag(freeindex(i))=0; end ID=(n*n); count = 0; count1 = 0; for i = 1:n*n if flag(i)==1% check if essential B.C count=count + 1; ID(i)=count;% arrange essential B.C nodes first d(count)=0;% store reordered essential B.C else count1=count1+1; ID(i)=e_nd+count1; end end temp_N=N; temp_x=x_T; temp_y=y_T; %rearrange the shape functions for i=1:n*n N(ID(i))= temp_N(i); x_T(ID(i))=temp_x(i); y_T(ID(i))=temp_y(i); end %------------------------------------------------------------------- %stiffness and force matrix %------------------------------------------------------------------- for i=1:n*n for j=1:n*n % K(i,j)=T*int(int(diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'),'x',-1,1),'y',-1,1); % M(i,j)=int(int(rho*N(i)*N(j),'x',-1,1),'y',-1,1); K1=T*diff(N(i),'x')*diff(N(j),'x')+diff(N(i),'y')*diff(N(j),'y'); K1_fn=matlabFunction(K1); K(i,j)=dblquad(K1_fn,-1,1,-1,1); M1=rho*N(i)*N(j); M1_fn=matlabFunction(M1); M(i,j)=dblquad(M1_fn,-1,1,-1,1); end F(i)=int(int(f*N(i),'x',-1,1),'y',-1,1); % F1=f*N(i); % F1_fn=matlabFunction(F1); % F(i)=dblquad(F1_fn,-1,1,-1,1); end %------------------------------------------------------------------- %extract the stiffness matrix and force matrix of free nodes %------------------------------------------------------------------- K_F=K(e_nd+1:n*n,e_nd+1:n*n); M_F=M(e_nd+1:n*n,e_nd+1:n*n); F_F=F(e_nd+1:n*n); K_FE=K(e_nd+1:n*n,1:e_nd); M_FE=M(e_nd+1:n*n,1:e_nd); d_E=d(1:e_nd,1); N_F=N(e_nd+1:n*n); x_F=x_T(e_nd+1:n*n); y_F=y_T(e_nd+1:n*n); for i=1:n*n-(4*n-4) % G(i)=int(int(rho*N_F(i)*(x_F(i)+1)*(y_F(i)+1/2)*cos(pi*x_F(i)/2)*cos(pi*y_F(i)/2),'x',-1,1),'y',-1,1); G1(i)=N_F(i)*rho*(x_F(i)+1)*(y_F(i)+1/2)*cos(pi*x_F(i)/2)*cos(pi*y_F(i)/2); G1_fn=matlabFunction(G1(i)); G(i)=dblquad(G1_fn,-1,1,-1,1); end %------------------------------------------------------------------ %solve for d_F [PHI,LAMDA]=eig(K_F,M_F); omega1=sqrt(LAMDA(1,1)); T1= (2*pi)/omega1; do=((M_F\G')); ddo=zeros(n*n-(4*n-4),1); F_FF=F_F'-K_FE*d_E; doo=[do ;ddo]; %------------------------------------------------------------------ %solve for d_F [T Y]=ode45('dynamic2',[0:0.3:2*T1],doo); %------------------------------------------------------------------ %d matrix for j=1:length(Y) for i=1:n*n if i<=e_nd d(i)=d_E(i); else d(i)=Y(j,i-e_nd); end end %----------------------------------------------------------------- %trial solution Uh=N*d; %----------------------------------------------------------------- %surf the figure [X,Z]=meshgrid(-1:0.1:1,-1:0.1:1); U=subs(Uh,{x,y},{X,Z}); surf(X,Z,U); colorbar; title('Analysis of Dynamic Case of Membrane'); zlim([-2,2]); xlabel('x'); ylabel('y'); zlabel('displacement'); M(j) = imgcf; hold off end % for i=1:length(Y) % str = strcat('D:\Graduate Courses\EML 5526 FEM\HW\hw7.2b\',num2str(i), '.jpg'); % A=imread(str); % [I,map]=rgb2ind(A,256); % if(i==1) % imwrite(I,map,'fe1.s11.team3.hw7.2.b.gif','DelayTime',0.01,'LoopCount',Inf) % else % imwrite(I,map,'fe1.s11.team3.hw7.2.b.gif','WriteMode','append','DelayTime',0.01) % end % end

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

Quadrilateral elements

Node numbers for the triangular case

Triangular elements

MATLAB Code

Matlab Code: clear all close all clc syms x y LX LY N K_global zi1 zi2 J load data.mat figure(1) plot(nodes(:,2),nodes(:,3),'r.'); for I=1:61 for J=1:61 K_global(I,J)=0*(x/x); end end %%%% basis function N(1) = (zi1-1)*(zi2-1)/4; N(2) = -(zi1+1)*(zi2-1)/4; N(3) = (zi1+1)*(zi2+1)/4; N(4) = -(zi1-1)*(zi2+1)/4; %%%%%% for e=1:48 E=elements(e,2:5); X=[nodes(E(2),2) nodes(E(3),2) nodes(E(4),2) nodes(E(1),2) ]; Y=[nodes(E(2),3) nodes(E(3),3) nodes(E(4),3) nodes(E(1),3)]; for i=1:4 Nxi(i,1)=diff(N(i),zi1); Nxi(i,2)=diff(N(i),zi2); end %%% Mr.jacob J=0*[x/x x/x; x/x x/x]; for i=1:4 J(1,1)=J(1,1)+(N(i)*X(i)); J(2,1)=J(1,1); J(2,2)=J(2,2)+(N(i)*Y(i)); J(1,2)=J(2,2); end J(1,1)=diff(J(1,1),zi1); J(1,2)=diff(J(1,2),zi1); J(2,1)=diff(J(2,1),zi2); J(2,2)=diff(J(2,2),zi2); J=inv(J)'; %%%% % stifness for i=1:4 for j=1:4 K(e,i,j) = ((Nxi(i,:)*J)*4*(Nxi(j,:)*J)')*det(J); K_fn=matlabFunction(K(e,i,j)); K(e,i,j)=dblquad(K_fn,-1,1,-1,1); % I=i+(3)*(e-1); % J=j+(3)*(e-1); K_global(E(i),E(j))=K_global(E(i),E(j))+K(e,i,j); K_global(E(i),E(j))=4*K_global(E(i),E(j)); end end for i=1:4 FF_fn=matlabFunction((Nxi(i,1)+Nxi(i,2))*det(J)); FF(E(i))=dblquad(FF_fn,-1,1,-1,1); end end % essential condition EC=[1,4,6,8,10,34,43,52,61]; % for k=1:length(EC) for i=1:61 for j=1:61 if i~=EC(k) || j~=EC(k) K_ff(i,j)=K_global(i,j); end end end end % F-matrix for i=1:61 if i==EC FF(i)=0; end end FF=FF'; D=double(inv(K_ff)*FF); uh=0*x; for e=1:48 for i=1:2 for j=1:2 I=i+(j-1)*2; uh=uh+N(e,I); end end end figure(3) ezsurf(uh,[0,1]) axis([0 1 0 1 -2 2])

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:

MATLAB Code

Matlab Code: clear all clc for m=2:1:3; num_x=vpa(ones(1,m)); %% to initialise the numerator matrix%% dem_x=vpa(ones(1,m)); %% to initialise the denominator matrix %% k=2/(m-1); %% to generate uniform discretization %% for i =1:m syms x %% defining the function as function of x %% for j= 1:m if i==j num_x(i)= num_x(i)*1; dem_x(i)=dem_x(i)*1; else num_x(i)=num_x(i)*((x+1-k*(j-1))); dem_x(i)=dem_x(i)*(k*(i-1)-k*(j-1)); end end b_x(i)= num_x(i)/dem_x(i); %% calculating the basis function %% end end for n=2:1:3; num_y=vpa(ones(1,n)); %% to initialise the numerator matrix%% dem_y=vpa(ones(1,n)); %% to initialise the denominator matrix %% k=2/(n-1); %% to generate uniform discretization %% for i =1:n syms y %% defining the function as function of x %% for j= 1:n if i==j num_y(i)= num_y(i)*1; dem_y(i)=dem_y(i)*1; else num_y(i)=num_y(i)*((y+1-k*(j-1))); dem_y(i)=dem_y(i)*(k*(i-1)-k*(j-1)); end end b_y(i)= num_y(i)/dem_y(i); %% calculating the basis function %% end end for j=1:n for i=1:m N(j,i)=b_x(i)*b_y(j); I=i+(j-1)*m; single_diff_x(I)=diff(N(j,i),x,1); single_diff_y(I)=diff(N(j,i),y,1); N_1(I)=N(j,i); N_2(I)=subs(N_1(I),x,-1); N_3(I)=subs(N_1(I),y,1); end end for I=1:m*n grad(I,:)=single_diff_x(I)*single_diff_x + single_diff_y(I)*single_diff_y; end K=double(int(int(grad,y,-1,1),x,-1,1))+int(0.5*N_2'*N_2,y,-1,1); F=-int(3*N_3',x,-1,1)+int(0.1*0.5*N_2',y,-1,1); K(1,1)=1e9; K(2,2)=1e9; K(3,3)=1e9; K(6,6)=1e9; K(9,9)=1e9; F(1)=K(1,1)*2; F(2)=K(2,2)*2; F(3)=K(3,3)*2; F(6)=K(6,6)*2; F(9)=K(9,9)*2; d=K\F; u_x_y=N_1(1)*d(1)+N_1(2)*d(2)+N_1(3)*d(3)+N_1(4)*d(4)+N_1(5)*d(5)+N_1(6)*d(6)+N_1(7)*d(7)+N_1(8)*d(8)+N_1(9)*d(9); ezsurf(u_x_y,[-1,1,-1,1]);

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