User:Anand ST14

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