User:Eml5526.s11.team01/Hwk7

=Problem 7.1: Solving on Bi-unit Square=

Given
$$\displaystyle \kappa =I$$, $$\displaystyle ~\rho c=3$$, $$\displaystyle g=2$$ on $$\displaystyle \partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$\displaystyle f =1$$, steady state.

Find
Solve the differential equation


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

$$\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |


 * $$\displaystyle (Eq. 1.1) $$
 * }
 * }

using LIBF

Solution The weak form:
 * {| style="width:100%" border="0"

$$\mathop{\iint }^{}w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{\kappa }_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)+f-\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }=0$$ Let’s calculate the derivative of the following factor
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.2) $$
 * }
 * }
 * {| style="width:100%" border="0"

$$\frac{\partial }\left( w{{K}_{ij}}\frac{\partial u} \right)=\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+w\frac{\partial }\left( {{K}_{ij}}\frac{\partial u} \right)$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.3) $$
 * }
 * }

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

$$w\frac{\partial }\left( {{K}_{ij}}\frac{\partial u} \right)=\frac{\partial }\left( w{{K}_{ij}}\frac{\partial u} \right)-\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.4) $$
 * }
 * }

Replacing in the integral we obtain
 * {| style="width:100%" border="0"

$$\mathop{\iint }^{}\left( \frac{\partial }\left( w{{K}_{ij}}\frac{\partial u} \right)-\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+wf\left( x,t \right)-w\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.5) $$
 * }
 * }

Appling the Gauss theorem on the first term we obtain
 * {| style="width:100%" border="0"

$$\mathop{\int }^{}w{{K}_{ij}}\frac{\partial u}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }+\mathop{\iint }^{}\left( -\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+wf\left( x,t \right)-w\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }=0$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.6) $$
 * }
 * }

Recognizing that the heat flux is given by


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

$${{q}_{i}}=-{{\text{ }\!\!\Kappa\!\!\text{ }}_{ij}}\frac{\partial u}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.7) $$
 * }
 * }

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

$$\mathop{\int }^{}-w{{q}_{i}}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }+\mathop{\iint }^{}\left( -\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)+wf\left( x,t \right)-w\rho c\frac{\partial u}{\partial t} \right)d\text{ }\!\!\Omega\!\!\text{ }=0$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.8) $$
 * }
 * }

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

$$\mathop{\iint }^{}\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}-\mathop{\int }^{}wH(u-{{u}_{\infty }})d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.9) $$
 * }
 * }


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

$$\mathop{\iint }^{}\frac{\partial w}\left( {{K}_{ij}}\frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}wHud{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}+\mathop{\int }^{}wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.10) $$
 * }
 * }

For the case when the conductivity is the identity matrix we have
 * {| style="width:100%" border="0"

$$\mathop{\iint }^{}\frac{\partial w}\left( {{\delta }_{ij}}\frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}wHud{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}+\mathop{\int }^{}wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.11) $$
 * }
 * }


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

$$\mathop{\iint }^{}\frac{\partial w}\left( \frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}wHud{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-\mathop{\int }^{}whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}+\mathop{\int }^{}wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.12) $$
 * }
 * }

But in this case we do not have natural boundary condition, therefore we have


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

$$\mathop{\iint }^{}\frac{\partial w}\left( \frac{\partial u} \right)d\text{ }\!\!\Omega\!\!\text{ }+\mathop{\int }^{}w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=\mathop{\iint }^{}wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.13) $$
 * }
 * }

Performing the discretization on the dominium we have $$u=\left[ N \right]{{u}_{e}}$$ and $$w=\left[ N \right]{{w}_{e}}$$, where $$\left[ N \right]$$is the vector of interpolation functions (Lagrange polynomial) and $$~{{u}_{e}}$$ and $${{w}_{e}}$$ are the nodal values o u and w respectively. The derivatives with respect to $${{x}_{1}}$$ and $${{x}_{2}}$$ are
 * {| style="width:100%" border="0"

$$~\left[ \begin{matrix} \frac{\partial u}{\partial {{x}_{1}}} \\ \frac{\partial u}{\partial {{x}_{2}}} \\ \end{matrix} \right]=\left[ \begin{matrix} \frac{\partial N}{\partial {{x}_{1}}} \\ \frac{\partial N}{\partial {{x}_{2}}} \\ \end{matrix} \right]{{u}_{e}}=\left[ B \right]{{u}_{e}}$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.14) $$
 * }
 * }

Replacing on the semi-discrete weak form we obtain


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

$${{w}^{eT}}\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }{{u}^{e}}+{{w}^{eT}}\mathop{\int }^{}\rho c{{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }{{\dot{u}}_{e}}={{w}^{eT}}\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.15) $$
 * }
 * }

From the last equation we can obtain the conductivity matrix, the heat source vector and the capacitance matrix. The conductivity matrix is given by


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

$$\overset{}{\mathop{k}}\,=\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle (Eq. 1.16) $$
 * }
 * }

The capacitance matrix is given by:


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

$$\overset{}{\mathop{M}}\,=~\mathop{\int }^{}\rho c{{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.17) $$
 * }
 * }

And the heat source is


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

$$\overset{}{\mathop{f}}\,=\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.18) $$
 * }
 * }

To solve the equation:


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

$$Ku+M\dot{u}=f$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.19) $$
 * }
 * }

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

$$\dot{u}={{M}^{-1}}\left( f-Ku \right)$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.20) $$
 * }
 * }

To solve to u from u ̇ we used the Crank-Nicolson scheme where:
 * {| style="width:100%" border="0"

$${{u}_{i+1}}={{u}_{i}}+\text{ }\!\!\Delta\!\!\text{ }t{{\dot{u}}_{i+\alpha }}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.21) $$
 * }
 * }

And:

For this scheme $$\alpha =\frac{1}{2}$$, therefore solving we have


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

$${{u}_{i+1}}={{u}_{i}}+\frac{\text{ }\!\!\Delta\!\!\text{ }t}{2}\left( {{{\dot{u}}}_{i}}+{{{\dot{u}}}_{i+1}} \right)$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.23) $$
 * }
 * }

Part 1: Steady state
$$\displaystyle \frac{\partial u}{\partial t}=0$$ and $$f\left( x \right)=1$$ using LIBF. To solve the problem using LIBF we divide the x axis and y axis in 2 parts initially. In that way the approximated solution using LIBF is given by:


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

$$u=\frac{\left( x-{{x}_{2}} \right)\left( x-{{x}_{3}} \right)}{\left( {{x}_{1}}-{{x}_{2}} \right)\left( {{x}_{1}}-{{x}_{3}} \right)}\frac{\left( y-{{y}_{2}} \right)\left( y-{{y}_{3}} \right)}{\left( {{y}_{1}}-{{y}_{2}} \right)\left( {{y}_{1}}-{{y}_{3}} \right)}{{u}_{1}}+...+\frac{\left( x-{{x}_{1}} \right)\left( x-{{x}_{2}} \right)}{\left( {{x}_{3}}-{{x}_{1}} \right)\left( {{x}_{3}}-{{x}_{2}} \right)}\frac{\left( y-{{y}_{1}} \right)\left( y-{{y}_{2}} \right)}{\left( {{y}_{3}}-{{y}_{1}} \right)\left( {{y}_{3}}-{{y}_{2}} \right)}{{u}_{9}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.24) $$
 * }
 * }

Where $${{u}_{i}}$$ is the nodal value at the node i. the numbering is given in the figure. The last equation can be written as $$u=\left[ N \right]$$.

And [N] can be replaced in the weak form to obtain the solution.

The obtained solution is shown in the next figure (obtained with the routine 1.m) Because lack of memory, the tolerance was only 10^-2

Part 2: Dynamic (Transient)
Solve the differential equation,


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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.25) $$
 * }
 * }

$$\displaystyle \kappa =I$$, $$\displaystyle ~\rho c=3$$, $$\displaystyle  g=2$$ on $$\displaystyle  \partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\displaystyle  \text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where  $$\displaystyle  f =1$$, steady state.

Similar results to the first part were obtained as shown in the figure. The routine 2.m was used to solve the problems. This code can be found in problem 1 appendix.

Part 3: Solving
Solve the differential equation,


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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.26) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, steady state, using LLEBF with NON-uniform mesh.

In this case the analysis was made using 100 elements and no convergence analysis was made. The problem was solved by using the matlab routine 3.m (can be found in the appendix) Inside that routine is included the file Nonunoformmesh.m which contains the connectivity, the coordinates of the nodes and the boundary conditions. The results are shown in the figure 1.4; it can be observed that the results are similar to those obtained with uniform mesh.



Part 4
Solve the differential equation


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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.27) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, transient; $$\text{u}\left( x,0 \right)\text{ }=\text{ xy}$$  using LIBF

Because of memory problems only 16 elements and a tolerance of 10^(-3) was used. Figure 1.5 shows the temperature of the center as function of time obtained with the Crank-Nicholson scheme. It is observed that the temperature tend to the value at the boundaries (2). It is observed also that the temperature decreased at the beginning; this error is probably because of the time step size which probably is too long.



5) Solve the differential equation


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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.28) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, transient; $$\text{u}\left( x,0 \right)\text{ }=\text{ xy}$$   using LLEBF; uniform mesh The initial condition is shown in figure 1.6



We used the Crank-Nicolson scheme before mentioned to solve the transient part of the problem. Initially we used a time step of 0.01 and the solution present some errors as can be observed in the following two figures.



Then we decreased the time step to 0.001 and better results were obtained as shown in the next 3 figures.





The next figure shows the temperature as function of time for the center of the domain point (0,0); it is observed that the temperature converge to the value of 2.07 which is the value at steady state.

The next figure shows the error before the convergence to the steady state. It was found that the error increases at the beginning. This instability could be probably eliminated decreasing more the time step; however that increase the computation time.

Part 6: Solving
Solve the differential equation,


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

$$div\left( \kappa .grad(u) \right)+f=\rho c\frac{\partial u}{\partial t}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 1.29) $$
 * }
 * }

With $$\kappa =I$$, $$~\rho c=3$$, $$g=2$$ on $$\partial \text{ }\!\!\Omega\!\!\text{ }$$, and $$\text{ }\!\!\Omega\!\!\text{ }$$ of a bi-unit square, where $$f =1$$, transient; $$\text{u}\left( x,0 \right)\text{ }=\text{ xy}$$  using LLEBF; NON-uniform mesh Solution The next figure shows the mesh used

The initial condition x*y is shown in the next figure



Because the value at steady state with f = 1 is close to the value at the borders (2.07) we increase the source term to f=20; in that way the steady state solution was as shown in the next figure

The next figure shows the evolution of the temperature at the center



The code 6.m, as found in the appendix, was used to solve this problem

Appendix
=Problem 2: Membrane=

Given
Given the membrane square; $$\displaystyle T = 4$$; $$\displaystyle f = 1$$; $$\displaystyle \rho =3$$; $$\displaystyle g = 0$$ on all the boundaries; bi-unit square; using LIBF

Find
Find the steady state and dynamic solution

Part 1

 * The general PDE is:


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

$$Tdiv\left( grad\left( u \right) \right)+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.1) $$
 * }
 * }


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

$$T\nabla .\left( \nabla u \right)+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.2) $$
 * }
 * }


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

$$T{{\nabla }^{2}}(u)+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.3) $$
 * }
 * }


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

$$T\frac{{{\partial }^{2}}u}{\partial x_{i}^{2}}+f=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.4) $$
 * }
 * }

The weak form is similar to the latter problem. The obtained result is


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

$${{w}^{eT}}\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }{{u}^{e}}+{{w}^{eT}}\mathop{\int }^{}\left( \frac{\rho }{T} \right){{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }{{\dot{u}}_{e}}={{w}^{eT}}\frac{1}{T}\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.5) $$
 * }
 * }

Therefore, the stiffness matrix is given by:


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

$$\overset{}{\mathop{k}}\,=\mathop{\iint }^{}{{\left[ B \right]}^{T}}\left[ B \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.6) $$
 * }
 * }

The mass matrix is:


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

$$\overset{}{\mathop{M}}\,=~\mathop{\int }^{}\left( \frac{\rho }{T} \right){{\left[ N \right]}^{T}}\left[ N \right]d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.7) $$
 * }
 * }

And the heat source is


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

$$\overset{}{\mathop{f}}\,=\frac{1}{T}\mathop{\iint }^{}{{\left[ N \right]}^{T}}f\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.8) $$
 * }
 * }

In the static case, the equation is the same as that for heat transfer in steady state; therefore we used the same routine 1.m for solve this problem. In the next figures is shown the solution and the error Convergence of 10^-6 was achieved with 1600 elements after 20 iterations



Part 2

 * Free vibration eigenvalue problem. Find the three first fundamental eigenpairs. For free vibration the factor f is zero, therefore the PDE is reduced to:
 * {| style="width:100%" border="0"

$$\rho \frac{{{\partial }^{2}}u}{\partial x_{i}^{2}}=\rho \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.9) $$
 * }
 * }

In 2D we can write:


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

$$\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}u}{\partial {{y}^{2}}}=\frac{\rho }{T}\frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.10) $$
 * }
 * }

This is case where analytical solution can be found; thererfor, we did not use any numerical method; instead, using separation of variables we can write the solution as:
 * {| style="width:100%" border="0"

$$u\left( x,y,t \right)=F\left( x \right)G\left( y \right)H(t)$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.11) $$
 * }
 * }

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

$$F(x)G\left( y \right)H\left( t \right)+F\left( x \right){{G}^{\left( y \right)}}H\left( t \right)=\frac{\rho }{T}F\left( x \right)G\left( y \right)H''(t)$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.12) $$
 * }
 * }

Where ‘’ means second derivative with respect to x, y or t depending if the derivative is of F, of G or of H respectively. If the latter equation is divided by u we obtain
 * {| style="width:100%" border="0"

$${{F}^{}}(x)G\left( y \right)H\left( t \right)+F\left( x \right){{G}^{}}(y)H\left( t \right)=\frac{\rho }{T}F\left( x \right)G\left( y \right)H''(t)$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.13) $$
 * }
 * }


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

$$\frac{F}{F}+\frac{G}{G}={{c}^{2}}\frac{H''}{H}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.14) $$
 * }
 * }

As every fraction in the latter equation depends only on x, y or t respectively they should be equal to a constant. Let’s call that constant w^2
 * {| style="width:100%" border="0"

$$\frac{F}{F}+\frac{G}{G}={{c}^{2}}\frac{H''}{H}=-{{w}^{2}}$$ Now we have the equations:
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.15) $$
 * }
 * }
 * {| style="width:100%" border="0"

$${{c}^{2}}\frac{H''}{H}=-{{w}^{2}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.16) $$
 * }
 * {| style="width:100%" border="0"
 * {| style="width:100%" border="0"

$$\frac{F}{F}=-\frac{G}{G}-{{w}^{2}}$$ In the last equation we again observe that the left hand side depends only on x and the left hand side depends on y therefore should be constant. Let’s call $$k_{1}^{2}$$and $$k_{2}^{2}$$ then
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.17) $$
 * }
 * }
 * {| style="width:100%" border="0"

$$\frac{F''}{F}=-k_{1}^{2}$$ And
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.18) $$
 * }
 * }
 * {| style="width:100%" border="0"

$$\frac{G''}{G}=~-k_{2}^{2}$$ Replacing equations
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.19) $$
 * }
 * }
 * {| style="width:100%" border="0"

$${{w}^{2}}=k_{1}^{2}+k_{2}^{2}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.20) $$
 * }
 * }

From above we have
 * {| style="width:100%" border="0"

$${{F}^{''}}+k_{1}^{2}F=0$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.21) $$
 * }
 * }

This equation has the solution $$F\left( x \right)=Acos\left( {{k}_{1}}x \right)+Bsin\left( {{k}_{1}}x \right)$$ To apply easily the boundary conditions we place the origin of coordinates at the lower left corner (Figure)

From the figure we observe that F(0) = 0, and F(L) = 0 With the first condition we have:
 * {| style="width:100%" border="0"

$$F\left( 0 \right)=Acos\left( 0 \right)+Bsin\left( 0 \right)=0$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.22) $$
 * }
 * }

Therefore A = 0. Also:
 * {| style="width:100%" border="0"

$$F\left( L \right)=Bsin\left( {{k}_{1}}L \right)=0$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.23) $$
 * }
 * }

In the last equation B cannot be zero because if B is zero, as A is zero F(x) would be zero and therefore u would be zero; therefore the $$\sin \left( {{k}_{1}}L \right)$$ has to be zero, which means that $${{k}_{1}}L=n\pi $$ where $$n=1,~2,~3,~\ldots $$ As in our case the length of the domain is 1 (L = 1) we have $${{k}_{1}}=n\pi $$ In similar way we can obtain that $${{k}_{2}}=m\pi $$ From above we have
 * {| style="width:100%" border="0"

$${{H}^{''}}\left( t \right)+\fracH\left( t \right)=0$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.24) $$
 * }
 * }

Therefore the solution is:
 * {| style="width:100%" border="0"

$$H\left( t \right)=Dcos\left( \frac{w}{c}t \right)+Esin\left( \frac{w}{c}t \right)$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.25) $$
 * }
 * }

From the latter equation we observe that the natural frequency is
 * {| style="width:100%" border="0"

$${{w}_{n}}=\frac{w}{c}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.26) $$
 * }
 * }

Replacing we have


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

$${{w}_{mn}}=\frac{1}{c}\sqrt{k_{1}^{2}+k_{2}^{2}}=\sqrt{\frac{T}{\rho }\left( {{\left( \frac{n\pi }{L} \right)}^{2}}+{{\left( \frac{m\pi }{L} \right)}^{2}} \right)}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.27) $$
 * }
 * }

The first 3 natural frequencies are when n = m=1; n = 1 and m = 2; n = 2 and m = 2

With T = 4, $$\rho =3$$, and L = 1 we obtain: w1 = 5.13; w2 = 8.11;  w3 = 10.26

Part 3

 * Transient analysis with f = 0; u(x,0) = u-static

In the dynamic case, the Newmark scheme was used to solve the equation:


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

$$\left[ M \right]\ddot{u}+\left[ K \right]u=F$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.28) $$
 * }
 * }

According to the Newmark scheme the function u and its first derivative are approximated by
 * {| style="width:100%" border="0"

$${{u}_{i+1}}={{u}_{i}}+\Delta t{{\dot{u}}_{i}}+\frac{1}{2}\Delta {{t}^{2}}{{\ddot{u}}_{i+\gamma }}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.29) $$
 * }
 * }
 * {| style="width:100%" border="0"

$${{\dot{u}}_{i+1}}={{\dot{u}}_{i}}+\Delta t{{\ddot{u}}_{i+\alpha }}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.30) $$
 * }
 * }

Where


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

$${{\ddot{u}}_{i+\gamma }}=\left( 1-\alpha \right){{\ddot{u}}_{i}}+\alpha {{\ddot{u}}_{i+1}}$$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right;">$$\displaystyle (Eq. 2.31) $$
 * }
 * }

The values of $$\alpha =\gamma =\frac{1}{2}$$where used. The problem was solved by using the routine 3.m in the appendix of problem 1. The next figure shows the displacement of the center of the membrane as function of time.

Part 4

 * Transient analysis with f = 0;

The next figure is the displacement vs time of the center.



'''Above are the Images from the Movie of the VIBRATION MEMBRANE. To see the full 1:51 minutes Movie go to the Following link.'''

=Problem 7.3 Static Solution for Unit Circle=

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

$$ T = 4 $$

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

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

Find
For both quad and triangle elements.

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

Plot deformed shape in 3-D.

Quad)
Using SolidWorks software. A unit circle was created, with a mesh of 9 elements, all quadrilateral. It was from the data points that the mesh nodes were created by hand using the measurement tool. The picture below shows the mesh construction on SolidWorks.



After obtaining the coordinates from the unit circle and numbering the nodes from 1 to 16 since there are 9 elements, the boundary conditions were defined.

The boundary on the bottom is defined as

and for the left hand side of the circle, the boundary condition is as follows:

This means that these boundaries are "free" and the nodes are open to deflection. However the arc of the circle is not "free" these boundary conditions are locked and the nodes stay in one place.

After obtaining this information, a Matlab code was created to re-create the nodes and boundary conditions for the structure. The code is listed below.

The following code then calls the mesh file and makes a 3D plot of the unit circle.

Here is the result of the 3D plot. The arc of the circle does not deflect, but the free ends of the circle are clearly displaced.



Triangular)
The procedure is done the same way except with a different Matlab code, shown below for the mesh and plot.



=Problem 7.4: Newtons Law of Colling=

Given
The partial Differential Equation is

1.Develop Weak form similar to (3)-(6) page 34-3.
The weak form will be as below,

Lets calculate the derivative of the following factor, Therefore, Replacing in the integral we obtain, Applying the Gauss Theorem on the first term we obtain, We can recognize that the Heat Flux is Given by, Eventually We have,

Organizing We get, Simplifying we get, For the Case when Conductivity if the Identity Matrix, we have following equations, Also We can Write as,

2) Develop Semidiscrete Equations (ODE's) similar to (3) page 23-3. Give Detailed Description for all the quantities.
Performing the Discretization on the Dominium we have,

Where N is the vector of the Interpolation (Lagrange Polynomial) and $$\displaystyle u_{e},w_{e}$$ are the Nodal Values of u & w Respectively. The Derivative With respect to $$\displaystyle x_{1}$$, $$\displaystyle x_{2}$$ are Replacing on the Semi Discrete Weak form we obtain, From the Last Equation we can Obtain the Conductivity Matrix, the Heat source vector and the Capacitance matrix. The Conductivity Matrix is Given by

Further the Capacitance Matrix is Given by,

And the Heat Source is

3.Solve G2DM2.0/D1 Using 2D LIBF (similar to HW 6.6) till $$10^{-6}$$ accuracy at center. Plot Solution in 3D w/Contour.
To solve G2DM2.0/D1 Following Routine was developed. The following figure shows the result for 16 points.

To solve the problem using LIBF we will devide the x-axis and Y-axis in 2 parts initially. In that way the Approximated solution by LIBF is given by:

We Have developed a COMSOL code to plot this routine. The COMSOL code is as follows

Appendix
The Plots for the Routine Developed along with the intended 3-D contour is as below.





Comment:

We can see from the above figures that there also exists a negative value, which is the same with our results from HW6.10.

Comment: The Contour is plotted for $$\displaystyle n_{el}=4$$. The above MATLAB code is developed for larger values of $$n_{el}$$.
=Contributing Members=