User:Eml5526.s11.team6/hwk7

=Problem 7.1 Solve G2DM1.0/D1 (Static and Transient) using 2D LIBF and LLEBF=

Given
The G2Dm1.0/D1 PDE

where: $$\left[ \right]$$ is an identity matrix The essential boundary condition is $$g(x,t) = 2\ \ \ \ on\ \ \ \ \partial \Omega $$.

Find
1a) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for steady state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

1b1) Solve G2DM1.0/D1 using 2D LLEBF with uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for steady state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

1b2) Solve G2DM1.0/D1 using 2D LLEBF with non-uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for steady state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

2a) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 0\ \ \ \ in\ \ \ \ \Omega = \square $$.

2b1) Solve G2DM1.0/D1 using 2D LIBF until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

2b2a) Solve G2DM1.0/D1 using 2D LLEBF with uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

2b2a) Solve G2DM1.0/D1 using 2D LLEBF with non-uniform mesh until accuracy $${10^{ - 6}}$$ at center (x,y)=(0,0) for transient state and $$f(x) = 1\ \ \ \ in\ \ \ \ \Omega = \square $$.

Solution
For Eq.6.5.1, we can write the weak form

1a) 2D LIBF
For 2D LIBF, we choose Lagrange Interpolation Basis Function for each element. For example, we choose n=m=2


 * $$N_{I}(x,y)=L_{i,m}(x)\cdot L_{j,n}(y),\;\; where \;\;I=i+(j-1)m$$

After substituting Eq.7.1.3 into Eq.7.1.2 in which Lagrange Interpolation Basis Function are both trial function and weight function, we integrate over the whole domain to obtain the capacitance, conductance, and flux matrices. These matrices are used to achieve final results for the temperatures at each node using the following Matlab codes.

1b) 2D LLEBF
For the 2D LLEBF case, we choose Linear Lagrange elemental basis functions for each element, i.e.,

After substituting Eq.7.1.4 into Eq.7.1.2 in which Linear Lagrange elemental basis function are both trial function and weight function, we integrate over the whole domain to obtain the capacitance, conductance, and flux matrices. These matrices are used to achieve final results for the temperatures at each node using the following Matlab codes.

1b1) Uniform Mesh
The first solution is calculated using a uniform mesh (a mesh where the nodes are equidistant from each other).



One sign that the solution was correctly coded is that the essential boundary conditions are met in these figures. Notice that along the edges of the temperature distribution graph the temperature is equal to 2. Additionally, because there is a positive heat source within the boundaries of our model we would expect the temperature to rise from the boundaries of the model, which is also apparent. These are good indicators that the solution provided is correct, or at least we are heading in the correct direction.

1b2) Non-Uniform Mesh
The second solution is calculated using a non-uniform mesh (an irregularly shaped mesh that displays the abilities of the robust code adopted from the textbook by Fish and Belytschko).



Similarly, the same boundary conditions are evident, and the rise in temperature in the center of the object is obvious.

2b2b) Non-Uniform Mesh
=Problem 7.2: Using 2D LIBF to Solve For Transverse Displacement of an Elastic Membrane=

Find
Solve the given problem using 2D LIBF until the relative error when evaluated at the point $$(x,y)=(0,0)$$ converges to O(10-6).

For the dynamic cases prepare a movie to show $$u^h(x,t)$$

Solution
2D LIBF take the following form:

where $$L_{i,m}(x)$$ and $$L_{j,n}(y)$$ are 1D LIBF with m and n nodes in the x and y directions, respectively.

Static Case
For the static case the strong form reduces to

Introducing a trial solution and weight function as follows

the discrete weak form becomes





Comments
The code still had unresolved errors at the time of the hw submission deadline, as such the desired convergence (or any for that matter) was not achieved. However the essential boundary condition was satisfied. This code is very computationally expensive as it depends on symbolic integration and therefore higher order basis functions were not able to be calculated due to the code timing out.

=Problem 7.3 Static Solution for Unit Circle with quadratic and triangular elements=

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

$$ \displaystyle T = 4 $$

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

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

Find
For both quad and triangle elements.

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

Plot deformed shape in 3-D.

Solution
By using MATLAB, quarter circle is created with 9 elements and the datapoints are also located and calculated using MATLAB. Following figure shows the mesing of the unit circle.



The co-ordinates are obtained using matlab code and the nodes are numbered from 1 to 16 for 9 elements.

The boundary conditions are defined for harizontal and vertical axes above as

Hence, we can say that nodes on the arc of the circle are fixed and don't have any displacement and remaining nodes can be moved. The MATLAB code is as follows:

The result of the MATLAB code in the form of Deformed 3D shape is shown below.



As we can see, the boundary of the circle is not deflected but other nodes show deflection.

The displacement contour is as shown below



=Problem 7.4-Solution to a general 2-D problem (Heat Transfer)with Lagrange Interpolation Basis Functions(LIBF) =

Given
Problem statement initiated under problem 6.10 (Re-Done here) in mtgs - [[media:Fe1.s11.mtg37.djvu|37-2]] and [[media:Fe1.s11.mtg37.djvu|38-2]]

2D space with boundary convection
Consider a heat conduction problem in 2-dimensional space with boundary convection as shown in the below Figure 6.10.1.



Boundary Conditions
Above 2D space has temperature $$ T(x,y)$$ as space unknown to solve for. We can state it as a boundary value problem with following boundary conditions specified at 3 different kind of boundaries as : $${\Gamma _g}$$, $${\Gamma _h}$$ and $${\Gamma _H}$$

Essential BC on $${\Gamma _g}$$
$$\displaystyle (Eq. 7.4.1) $$

Note : $$T (x,y)$$ can also can be defined as a varying quantity along $${\Gamma _g}$$

Natural BC on $${\Gamma _h}$$
Refer to the Figure 7.4.2 below. We specify the heat flux on boundary $${\Gamma _h}$$ as $${q_n} = {q_0}$$.

This is as we have : $$ {q_n} = \overline q .\overline n $$ with $$\overline n $$ as the unit normal vector on $${\Gamma _h}$$ Hence we can also write : $${q_n} = ( - k.\overline \nabla T).\overline n  =  - k(\overline \nabla  T.\overline n )$$ OR we then have :

$$\displaystyle (Eq. 6.10.2) $$ Note- $$ \nabla T$$ is positive as explained in the - Assumption- subsection below

We are to consider only the normal component of $$\overline q $$, i.e. $$ {q_n}$$, for heat transactions as that is the only component causing heat flow in or out of the element. The tangential component of $$\overline q $$, i.e. $$ {q_t}$$, is just circulating over the boundary.

Assumption
$$\overline n  $$ is unit vector pointing in the outward direction. This implies that $$ \nabla T$$ is positive and temperature increases as we go out away from the surface towards inside of the body. Also we note that $$\overline n  $$ is a vector quantity normal to the $$ T(x,y)$$ contours in 2D space. As $$ T(x,y)$$ increases from inside to outside, heat flux $${q_n}$$ goes into the volume from outside of the boundary $${\Gamma _h}$$ This implies that when $${q_0}$$ is positive, heat comes into the system

Convection BC on $${\Gamma _H}$$
We have as the convection boundary where heat transfer occurs via heat convection Hence we can say that the normal component $${q_n}$$ of the heat flux $$\overline q $$ equals the heat transfer due to heat convection at $${\Gamma _H}$$. Hence we write-

$${q_n} = h(T - Ta)$$ where, we have - $$T$$ : Temperature at the surface of $${\Gamma _H}$$ $${T_a} = {T_\infty } =$$ Say, Ambient temperature $$h$$ :Coefficient of heat convection at $${\Gamma _H}$$ Thus :

$$\displaystyle (Eq. 7.4.3) $$

Given domain with 2-dimensional model data set 2
The following data set needs to be solved using 2D LIBFs -

Solution
Referenced from Section 6.3 (pg 142 )and Section 8.1 (pg 189 ) of text book - ' A first course in finite element method' by authors Fish & Belytschco

Deriving the Equilibrium Equation for 2D- heat transfer
We first derive the equilibrium equation using elemental plane surface as shown in Figure- 7.4.3 below. The elemental surface has length as 'dx' and height as 'dy'. The inflow and outflow of heat is shown symbolically with arrows at various sections.

Note- Heat flux ' $$ \underline q $$ ' is a vector quantity in real 2D space defined per unit surface area and also $$f = $$ Heat generation per unit volume for the given elemental surface We assume unit dimension along the direction perpendicular to the screen or paper

Now the equilibrium of heat can be written in following manner -

Now we use | Taylor Series Expansion technique and can then write -

$$ {({q_x})_{x + dx}} = {({q_x})_x} + {\left( {\frac} \right)_x}dx + higher\_order\_terms$$ OR we can say with the higher order terms neglected :

Similarly we can write for y- direction -

And now $$ Eq.7.4.4 $$ can be modified as below - Canceling out $$dxdy $$ throughout as we have used it for elemental surface and it has to be true over the entire 2D space,

$$ Eq. 7.4.8 $$ can be written with different notation in the form of Divergence of $$\overline q $$ as - Note that $$ Eq. 7.4.9 $$ is valid for 2D as well as 3D heat transfer with $$\overline q $$ as 2D or 3D vector respectively in real 2D or 3D space. Here we are solving for temperature $$T(x,y)$$ so we use 2D space with $$T(x,y)$$

Constitutive Equation for 2D- heat transfer
Here we use the Fourier Law for 2D space which is stated as : $$\overline q  =  - \left[ k \right]\overline \nabla  T $$ It can be written as - Note- Usually, $$=0 $$ and $$=0 $$ for most of the materials as we do not see flow of heat perpendicular to the direction of the temperature gradient line. Also, we have : $$= $$ for isotropic materials and $$ \ne  $$ for anisotropic materials We consider isotropic materials here hence we can assume : $$= $$. Hence $$ Eq. 7.4.10 $$ now becomes -

Hence we write then :

Governing Equation for 2D- heat transfer
Using $$ Eq.7.4.9 $$ and $$ Eq. 7.4.10 $$ we can say that :

Thus the strong form for 2D - Heat Transfer can be written as follows-

$$\displaystyle (Eq. 7.4.11) $$

Deriving the Weak Form for 2D- heat transfer with boundary convection
Let us select the weighting function $$ w(x,y) $$ such that $$ w(x,y)=0$$ at boundary $$ {{\Gamma }_{g}} $$ where the essential boundary condition is applied. This is actually obvious to say as the essential boundary condition is specified at $$ {{\Gamma }_{g}} $$, the value of $$ T (x,y)$$ will be constant there and will no longer require any weighting function multiplication to approximate it in the given domain. Also this will eliminate the unknown flux at the boundary $${{\Gamma }_{g}} $$

Integrating over volume (for general 3D - case) after multiplying with the weighting function w(x), in the form of weighted residual, we get -

Note- For 2D- case : Put $$ dz $$ = $$1$$ which will give us the Area Integral Now we use Green's Theorem to perform the integration by parts. Our aim to get rid of the derivative on $$ \overline{\nabla }.({{K}_{ij}}\overline{\nabla }T)$$ and a derivative should appear on $$ w(x) $$ We state the general Green's Theorem for scalar fields $$ \varphi $$, T - is as follows -

Here, $${{\nabla }^{2}} $$ serves as Laplacian or Divergence of the gradient. Also, $$\overline{n} $$ is the unit normal at boundary and is variable at each point on the boundary. Thus ,we have a volume integral converted to a volume plus a surface integral Hence, $$ Eq. 6.10.12 $$ can now be modified as follows -

This is general equation for 3D as well as 2D case and can be converted to specific case depending upon whether the gradient is 2D or 3D Now again we note that in the weak form the natural boundary conditions get pulled in automatically, for example- the surface integral in above equation will be split into integral over the 3 boundaries - $${{\Gamma }_{g}} $$, $${{\Gamma }_{H}} $$ , and$${{\Gamma }_{h}} $$ , as follows :

Note- In the above Equation, the first term on the RHS will vanish as the weighting function $$ w (x) $$ vanishes on $$ {{\Gamma }_{g}}$$ as explained. Also here we flipped the signs for convection boundary condition on boundary on $${{\Gamma }_{g}} $$ Now, we put $$ Eq. 7.4.15 $$ into $$ eq. 6.10.14$$ and keep the terms required to get the stiffness matrix $$ \left[ K \right]$$ along with the terms involving temperature on LHS and all the terms with known values on the RHS.

$$\displaystyle (Eq. 7.4.16) $$ $$ Eq. 7.4.16 $$ is the weak form for the given 2D Heat Transfer domain with convection boundary specified as $$ {{\Gamma }_{H}}$$ This equation is of the following form -

It is the continuous form of the weak equation and we need to convert it into a discrete form using the method of Desensitization in FEM. Also note that - 1. In $$ Eq. 7.4.16 $$, temperature is a continuous function of field parameters (x,y). 2. This equation is good for 2D as well as 3D domains. 3. We can handle anisotropic materials by converting conductivity 'k' into matrix form as explained.

Developing the semi-discrete equations
Now the temperature distribution in 2D space can be written as follows with subscript 'e' denoting to an element i.e. temperature at a point (x,y) in the 2D space is a sum over that of the elemental nodes. Also here $$ N_i $$ denotes shape function and $$ T_i $$ denotes nodal temperature values.

This can be written in matrix form as follows where the nodal values matrix is shown by $$\left\{ {{X}_{e}} \right\} $$

Also applying the same for the weighting function, we can write : This can be written in matrix form as follows -

We note that we have got gradient terms in the weak form as seen in $$Eq. 7.4.16 $$. The $$\overline{\nabla }T $$ can be written as :

Here $$ \left[ B \right]$$ denotes the matrix of derivatives of the shape functions as shown in below $$ Eq. 7.4.22 $$.

Similar equation can be written for the weighting function as follows -

Now the LHS of the weak form in $$ Eq. 7.4.16 $$ can be written as for a single element -

Hence, LHS becomes for $$ Eq. 7.4.16$$ as :

Depending upon various boundary conditions the RHS of $$ Eq. 7.4.16 $$ can be written as $$\left[ {{F}_{e}} \right]{{\left\{ \partial {{X}_{e}} \right\}}^{T}} $$ Hence the discrete form Equation for first term of the stiffness matrix becomes -

where, $$ \left[ {{K}_{e1}} \right]=\int\limits_{V}{{{\left[ B \right]}^{T}}\left[ K \right]}\left[ B \right]dv$$

The second term of the stiffness can be written in discrete form as follows-

Now the mass matrix can be written in discretized form from the first term of the continuous weak form of $$ Eq. 7.4.16 $$ as follows -

Now the first term on RHS can be discretized as -

Now the second term on RHS can be discretized as -

And the last term on RHS of the weak form can be discretized as -

Capacitance matrix

Conductivity matrix

Heat source matrix

Semi- discrete equations form the given data set G2DM2.0/D1
Referring to mtg - [[media:Fe1.s11.mtg37.djvu|41-5]], the discrete form equations for the given data set as explained in the 'Given' section above, can be written as follows -

Note- In all the matrices below, we have overlap at Node 13 of the domain as shown in Figure 7.4.5 below -

Now the conduction contribution from $${{\Gamma }_{H}} $$ can be written as follows -

Solve the weak form
The following Matlab codes are used to solve the weak form given above using LLEBF.

compare results for different numbers of elements
We use uniform elements and the results of different numbers of elements are as following









From the above results, we can find these results are identical to each other, which means that our codes are correct to achieve the same results even with different numbers of elements.

We can find "negative temperature" at the intersection corner of natural boundary condition and mixed boundary condition. Although it's quite unphysical to have "negative temperature", we try to adjust the boundary condition to find out the reason why the "negative temperature" exists.

compare results for different essential boundary conditions
First we compare the results of different essential boundary conditions with g=2,4,10,300









From the above results, we can find that we won't have "negative temperature" if g>4. From the 2D temperature distribution, we can find that whatever the essential boundary is, the distribution is similar with each other but the magnitude of maximum and minimum are different. However, the distance between maximum and minimum temperature is identical to each other, which means that the essential boundary can "shift" the results. Therefore maybe the "negative temperature" exists because the essential boundary is a little low for this case.

compare results for different natural boundary conditions
Then we compare the results of different natural boundary conditions with h=6,3,2,1









From the above results, we can find that the lower outward heat flux is, the minimum temperature is, because with lower outward heat flux, the temperature will dissipate slowly. When h = 1, the "negative temperature" disappears because less and less temperature will flow out of the area.

Why do we have "negative temperature"?
It's true that the "negative temperature" exist without any physical meaning, because in the universe the minimum temperature is zero K where "kelvin" is a unit of measurement for temperature. But if we suppose the unit of the temperature is Celsius scale, which is another unit scale for thermal temperature, we can have negative temperature just like what we have in winter. In other word, the actual governing equation we solved in this problem is a Laplace equation which can not guarantee that all the results on all the nodes are not negative, unless we can give it the "true" boundary condition which are physical in real world.

And also from the above procedures that we compare results from different boundary conditions, we can find that the boundary conditions do dominate the temperature distribution and "negative temperature" disappears with certain boundary conditions.

Therefore, We believe that we can have physical temperature distribution with physical boundary conditions, even sometimes "negative temperature" can exist with using different aspects of measure unit.

=Contributing Members=

Eml5526.s11.team6.vork 17:52, 21 April 2011 (UTC)

Eml5526.s11.team6.deshpande 19:27, 21 April 2011 (UTC)

Eml5526.s11.team6.kurth 20:31, 21 April 2011 (UTC)

Eml5526.s11.team6.gravois 20:49, 21 April 2011 (UTC)

Eml5526.s11.team6.tupsakhare 20:54, 21 April 2011 (UTC)

eml5526.s11.team6.lee 20:58, 21 April 2011 (UTC)

Eml5526.s11.team6.joglekar 23:01, 21 April 2011 (UTC)