User:Eml5526.s11.team7.jin/HW6

=Problem 6.2 Solving the General 1D problem with Quadratic Lagrangian element basis function (QLEBF)=

Problem Description
Similar to HW5.{2,4,8}, solve G1DM1.0/D1b using Quadratic Lagrangian Element Basis Function (QLEBF) with uniform discretization (equidistant element nodes) until convergence of $${{u}^{h}}(0.5)$$ to $$O({{10}^{-6}})$$. ( $$nel=2,4,6,8,...$$).


 * 1. For $$nel=2$$, compute $$\mathbf{\tilde{K}}=\sum\limits_{e=1}^{2}$$, with $${{\mathbf{\tilde{K}}}^{e}}$$ by (1)p.30-5, display $${{\mathbf{\tilde{K}}}^{e}},e=1,2$$.


 * 2. Compute $${{\mathbf{k}}^{e}}$$, $${{\mathbf{L}}^{e}}$$, for $$e=1,2$$.


 * 3. Compute $${{\tilde{\Kappa }}^{e}}={{\mathbf{L}}^{e}}^{T}{{\mathbf{k}}^{e}}{{\mathbf{L}}^{e}}$$, for $$e=1,2$$, compare to the result by (1).


 * 4. Plot all QLEBF for $$nel=3$$.


 * 5. Plot $$u_^{h}$$ vs. $$u$$, $$\left[ u_^{h}(0.5)-u(0.5) \right]$$ vs. $$\tilde{n}$$

Part 1
For element 1, the 3 basis functions are:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_1^1 = \frac \cdot \frac = 8{x^2} - 6x + 1

$$     (6.2.1)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_2^1 = \frac \cdot \frac = - 16{x^2} + 8x

$$     (6.2.2)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_3^1 = \frac \cdot \frac = 8{x^2} - 2x

$$     (6.2.3) By (1)p.30-5,
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{ij}^e = \int\limits_ {b{{_i^e}^\prime }(x){a_2}(x)b{{_j^e}^\prime }(x)dx}

$$     (6.2.4) Thus, we have
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{11}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( {16x - 6} \right)dx} = \frac{6}

$$     (6.2.5)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{12}^1 = K_{21}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( { - 32x + 8} \right)dx} =  - \frac{3}

$$     (6.2.6)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{13}^1 = K_{31}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} = \frac{6}

$$     (6.2.7)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{22}^1 = \int\limits_0^{0.5} {\left( { - 32x + 8} \right)\left( {2 + 3x} \right)\left( { - 32x + 8} \right)dx} = \frac{3}

$$     (6.2.8)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{23}^1 = K_{32}^1 = \int\limits_0^{0.5} {\left( { - 32x + 8} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} =  - \frac{3}

$$     (6.2.9)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{33}^1 = \int\limits_0^{0.5} {\left( {16x - 2} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} = \frac{6}

$$     (6.2.10) Then
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde K}}^1} = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{89/6}&0&0 \\   0&0&0&0&0 \\   0&0&0&0&0 \end{array}} \right]

$$     (6.2.11) Similar to element 1, when  e=2
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde K}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0 \\  0&0&0&0&0 \\   0&0&{107/6}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.12) Thus,
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{\mathbf{\tilde K}} = \sum\limits_{e = 1}^2  = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{98/3}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.13)
 * 
 * }

Part 2
By (1)p.31-3,
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

k_{ij}^e = \int\limits_ {b{{_i^e}^\prime }(x){a_2}(x)b{{_j^e}^\prime }(x)dx}

$$     (6.2.14) For element 1
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{k}}^1} = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6} \\  { - 38/3}&{88/3}&{ - 50/3} \\   {11/6}&{ - 50/3}&{89/6} \end{array}} \right]

$$     (6.2.15)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{L}}^1} = \left[ {\begin{array}{ccccccccccccccc} 1&0&0&0&0 \\  0&1&0&0&0 \\   0&0&1&0&0 \end{array}} \right]

$$     (6.2.16) For element 2
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{k}}^2} = \left[ {\begin{array}{ccccccccccccccc} {107/6}&{ - 62/3}&{17/6} \\  { - 62/3}&{136/3}&{ - 74/3} \\   {17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.17)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{L}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&1&0&0 \\  0&0&0&1&0 \\   0&0&0&0&1 \end{array}} \right]

$$     (6.2.18)
 * <p style="text-align:right">
 * }

Part 3
For element 1
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde {\rm K}}}^1} = {{\mathbf{L}}^1}^T{{\mathbf{k}}^1}{{\mathbf{L}}^1} = \left[ {\begin{array}{ccccccccccccccc} 1&0&0 \\  0&1&0 \\   0&0&1 \\   0&0&0 \\   0&0&0 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6} \\  { - 38/3}&{88/3}&{ - 50/3} \\   {11/6}&{ - 50/3}&{89/6} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 1&0&0&0&0 \\  0&1&0&0&0 \\   0&0&1&0&0 \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{89/6}&0&0 \\   0&0&0&0&0 \\   0&0&0&0&0 \end{array}} \right]

$$     (6.2.19) For element 2
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde {\rm K}}}^2} = {{\mathbf{L}}^1}^T{{\mathbf{k}}^1}{{\mathbf{L}}^1} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0 \\  0&0&0 \\   1&0&0 \\   0&1&0 \\   0&0&1 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {107/6}&{ - 62/3}&{17/6} \\  { - 62/3}&{136/3}&{ - 74/3} \\   {17/6}&{ - 74/3}&{131/6} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0&0&1&0&0 \\  0&0&0&1&0 \\   0&0&0&0&1 \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0 \\  0&0&0&0&0 \\   0&0&{107/6}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.20) The result is exactly same as achieved in part 1.
 * <p style="text-align:right">
 * }

Part 4
The plot of QLEBF for nel = 3:

Part 5
The plots of $$u_^{h}$$ vs. $$u$$,



The plots of $$\left[ u_^{h}(0.5)-u(0.5) \right]$$ vs. $$\tilde{n}$$

=Problem 6.9 Finite Element Formulation for Multidimensional Scalar Field Problems=

Problem Description
Consider the heat conduction problem depicted in Figure 6.9.1. The coordinates are given in meters. The conductivity is isotropic, with $$\mathbf{D}=k\left[ \begin{matrix} 1 & 0 \\   0 & 1  \\ \end{matrix} \right]$$, and $$k=5\text{W}{}^{\text{o}}{{\text{C}}^{\text{-1}}}$$. The temperature $$T=0$$ is prescribed along edges AB and AD. The heat fluxes $$\bar{q}=0$$ and $$\bar{q}=20\text{W}{{\text{m}}^{\text{-1}}}$$ are prescribed on edges BC and CD, respectively. A constant heat source $$s=6\text{W}{{\text{m}}^{\text{-2}}}$$ is applied over the plate. Consider this problem modeled with 16 quadrilateral finite elements as shown in Figure 6.9.2.

Solve this problem using the finite element code given in Section 12.5.



Solution
By using the matlab codes downloaded from download, The post-processing results for temperature and flux are shown in Figure 6.9.3 and Figure 6.9.4.

All the matlab codes for this problem can be found in the file folder 'ch08'.

The main program is given in heat2d.m file. It calls the following functions: include_flags, [K,f,d] = preprocessor, [ke, fe] = heat2Delem(e), [K,f] = assembly(K,f,e,ke,fe), f = src_and_flux(f), [d,f_E] = solvedr(K,f,d) and postprocessor(d).

In the function of [K,f,d] = preprocessor, input file, which defines material properties, mesh data, vector and matrices initializations, essential and natural conditions and mesh generation, is read.

In the function of [ke, fe] = heat2Delem(e), the element conductance matrix and nodal source vector are integrated using Gauss quadrature. And then assembled by the function of [K,f] = assembly(K,f,e,ke,fe).

The f = src_and_flux(f) function calculates the nodal boundary flux vector.

In the postprocessor(d) function, the results are plotted in Figure 6.9.3 and Figure 6.9.4.

Problem Description
Refer to p205-207 of F&B for detailed problem description.

Problem Description
Consider a chimney constructed of two isotropic materials: dense concrete ($$k=2.0W{}^{\circ }{{C}^{-1}}$$) and bricks ($$k=0.9W{}^{\circ }{{C}^{-1}}$$). The temperature of the hot gases on the inside surface of the chimney is $$140{}^{\circ }C$$, whereas the outside is exposed to the surrounding air, which is at $$T=10{}^{\circ }C$$. The dimensions of the chimney (in meters) are shown below. For the analysis, exploit the symmetry and consider 1/8 of the chimney cross sectional area. Consider a mesh of eight elements as shown below. Determine the temperature and flux in the two materials. Analyze the problem with 2×2, 4×4 and 8×8 quadrilateral elements for 1/8 of the problem domain. A2×2 finite element mesh is shown in Figure below. Symmetry implies insulated boundary conditions on edges AD and BC. Note that element boundaries have to coincide with the interface between the concrete and bricks.