# Nonlinear finite elements/Homework 4/Solutions
## Problem 1: Axially Loaded Bar with Nonlinear Modulus

Given:
Bar of uniform cross-section under axial body load and compressive axial end load ($R$). Nonlinear material with constitutive relation
    $\sigma =E_{0}(1-\varepsilon )\varepsilon ,\qquad \varepsilon ={\cfrac {du}{dx}}~.$
Finite element model uses quadratic elements with shape functions
    ${\begin{aligned}N_{1}(x)&={\cfrac {2(x_{2}-x)(x_{3}-x)}{h^{2}}}\\N_{2}(x)&={\cfrac {4(x-x_{1})(x_{3}-x)}{h^{2}}}\\N_{3}(x)&={\cfrac {2(x-x_{1})(x-x_{2})}{h^{2}}}\end{aligned}}$
where $h=x_{3}-x_{1}$.

### Solution
#### Part 1

Plot the shape functions for an element. Show that they add up to 1 at all points within the element.
A plot of the shape functions for an element is shown in Figure 1(a). The shape functions for a mesh are shown in Figure 1(b). If we add the shape functions we get
    ${\begin{aligned}\sum _{i}N_{i}(x)&={\cfrac {2(x_{2}-x)(x_{3}-x)}{h^{2}}}+{\cfrac {4(x-x_{1})(x_{3}-x)}{h^{2}}}+{\cfrac {2(x-x_{1})(x-x_{2})}{h^{2}}}\\&={\cfrac {2}{h^{2}}}\left[(x_{2}-x)(x_{3}-x)+(x-x_{1})(x_{3}-x)+(x-x_{1})(x_{3}-x)+(x-x_{1})(x-x_{2})\right]\\&={\cfrac {2}{h^{2}}}\left[(x_{2}-x_{1})(x_{3}-x)+(x-x_{1})(x_{3}-x_{2})\right]\\&=\left({\cfrac {2}{h^{2}}}\right)\left({\cfrac {h}{2}}\right)(x_{3}-x_{1})=1\qquad \square \end{aligned}}$
| Figure 1 (a). Shape functions for an element. | Figure 1 (b). Shape functions for a mesh. |

#### Part 2

Plot the displacement, stress, strain, and the error norm when you use 10 elements to discretize the bar. Discuss and comment on your observations.
If we use an isoparametric form of the equations to simplify the integration of the stiffness matrix and load vector terms, the shape functions become
    ${\begin{aligned}N_{1}(\xi )\,&=-{\frac {1}{2}}\left(1-\xi \right)\xi \\N_{2}(\xi )\,&=\,\left(1+\xi \right)\left(1-\xi \right)\\N_{3}(\xi )\,&={\frac {1}{2}}\left(1+\xi \right)\xi ~.\end{aligned}}$
The map between $x$ and $\xi$ is
    $x(\xi )=x_{1}N_{1}(\xi )+x_{2}N_{2}(\xi )+x_{3}N_{3}(\xi );~\qquad x_{2}={\cfrac {x_{1}+x_{3}}{2}}~.$
Therefore, the Jacobian of the deformation is
    $J={\cfrac {dx}{d\xi }}={\cfrac {h}{2}}~.$
The terms of the stiffness matrix are given by
    ${\begin{aligned}K_{ij}&=\int _{x_{1}}^{x_{3}}AE_{0}\left(1-\sum _{k=1}^{3}u_{k}{\cfrac {dN_{k}(x)}{dx}}\right){\cfrac {dN_{i}(x)}{dx}}{\cfrac {dN_{j}(x)}{dx}}~dx\\&=\int _{-1}^{1}AE_{0}\left[1-\sum _{k=1}^{3}u_{k}\left(J^{-1}{\cfrac {dN_{k}(\xi )}{d\xi }}\right)\right]\left(J^{-1}{\cfrac {dN_{i}(\xi )}{d\xi }}\right)\left(J^{-1}{\cfrac {dN_{j}(\xi )}{d\xi }}\right)~J~d\xi ~.\end{aligned}}$
The element stiffness matrix is
    $\mathbf {K} ={\cfrac {AE_{0}}{3h^{2}}}{\begin{bmatrix}7h+u_{3}+15u_{1}-16u_{2}&&-8(h+2u_{1}-2u_{2})&&h-u_{3}+u_{1}\\-8(h+2u_{1}-2u_{2})&&16(h-u_{3}+u_{1})&&-8(h-2u_{3}+2u_{2})\\h-u_{3}+u_{1}&&-8(h-2u_{3}+2u_{2})&&7h-15u_{3}-u_{1}+16u_{2}\end{bmatrix}}$
The terms of the tangent stiffness matrix are given by
    $T_{ij}=K_{ij}+\sum _{m=1}^{3}{\frac {\partial K_{im}}{\partial u_{j}}}u_{m}~.$
Therefore, the element tangent stiffness matrix is
    $\mathbf {T} ={\begin{bmatrix}7h+2u_{3}+30u_{1}-32u_{2}&&-8(h+4u_{1}-4u_{2})&&h-2u_{3}+2u_{1}\\-8(h+4u_{1}-4u_{2})&&16(h-2u_{3}+2u_{1})&&-8(h-4u_{3}+4u_{2})\\h-2u_{3}+2u_{1}&&-8(h-4u_{3}+4u_{2})&&7h-30u_{3}-2u_{1}+32u_{2}\end{bmatrix}}$
The terms of the element load vector are given by
    ${\begin{aligned}f_{i}&=\int _{x_{1}}^{x_{3}}axN_{i}(x)~dx+\left.N_{i}\left(AE{\cfrac {du}{dx}}\right)\right|_{x_{1}}^{x_{3}}\\&=\int _{-1}^{1}a\left(\sum _{k=1}^{3}x_{k}N_{k}(\xi )\right)N_{i}(\xi )~J~d\xi +\left.N_{i}(\xi )R_{i}\right|_{-1}^{1}~.\end{aligned}}$
Therefore, the element load vector is
    $\mathbf {f} ={\cfrac {ah}{6}}{\begin{bmatrix}x_{1}\\2(x_{1}+x_{3})\\x_{3}\end{bmatrix}}+{\begin{bmatrix}-R_{1}\\0\\R_{3}\end{bmatrix}}$
The stress-strain curve in compression is shown in Figure 2.
* Figure 2. Stress-strain curve in compression.
If we use a mesh containing $2$ equally sized elements, we get the solutions shown in Figure 3.
| Figure 3(a). Displacement | Figure 3(b). Stress |
| Figure 3(a). Strain | Figure 3(b). Error norm |
If we use a mesh containing $10$ equally sized elements, we get the solutions shown in Figure 4.
| Figure 4(a). Displacement | Figure 4(b). Stress |
| Figure 4(a). Strain | Figure 4(b). Error norm |
Observations:
* The two-node and three-node elements give the same solution.
* The solution using two elements shows that the stresses and strains are interpolated linearly within each element and are non longer constant.
* There is no jump in the stress/strain at the nodes.
* The solution converges quite rapidly - four iterations. This is due to the shape of the stress-strain curve in compression.

#### Part 2

Do the same for 100 elements. What are your observations on the effect of mesh refinement on the solution?
If we use a mesh containing $100$ equally sized elements, we get the solutions shown in Figure 5.
| Figure 5(a). Displacement | Figure 5(b). Stress |
| Figure 5(a). Strain | Figure 5(b). Error norm |
Observations:
* The solution is approximated quite accurately with 10 elements. Using 100 elements does not lead to any significant improvement in the solution.
* The solution converges in four iterations unlike the case where a tensile load was applied.

#### Part 4

Now, increase the applied load. At what load does the solution fail to converge? Explain why the solution fails to converge at this load.
The solution converges at all loads unlike in the tension problem. If the load is increased to $R=-20.1$, the elements invert and we get a negative Jacobian. Even though we get a solution, that solution is unphysical.
The Matlab code use to compute the results is shown below.
You could alternatively write your own Maple code or use ANSYS for these calculations depending on your interests.

## Problem 2: Nonlinear Thermal Conduction

Consider the model chip system shown in Figure 6.
* Figure 6. Schematic of a model chip.
Assume that heat is being generated in the silicon carbide region of the chip at the rate of 100 W/cm$^{2}$. The cooling channels contain a fluid at a constant 250 K and (supposedly) can dissipate heat at the rate of 800 W/cm$2$. The surroundings are at room temperature (300 K).
The thermal conductivity of the silicon carbide and the aluminum nitride are also given in Figure 6. The thermal conductivity of the epoxy resin is 1.5 W/m-K.
* Figure 7. Meshed model

#### Part 1

First run a linear simulation of the problem using room temperature thermal conductivities of the materials. This will act as a check of the model and setup and will help you determine any additional constraints that you will need. Do the cooling channels really dissipate heat at 800 W/cm$2$?
The following ANSYS input is used for this simulation. Due to symmetry, only half of the model will be meshed Figure 7. A constant temperature of 300K is applied around the edge Epoxy resin. Heat is being generated from areas made up by the Silicon Carbide at the rate of 100 W/cm$^{2}$. Each line forming the cooling channels has constant temperature of 250K.
Shown in Figure 8 and 9 below are the contour plots of temperature and the vector sum of heat flux. From Figure 9 we see that the cooling channels do not dissipate heat at the rate suggested in the problem statement.
| Figure 8. Temperature of the linear thermal conductivity model (K). | Figure 9. Flux of the linear thermal conductivity model (W/m-K). |

#### Part 2

Fit curves to the thermal conductivity versus temperature data for Silicon Carbide and Aluminum Nitride.
Below are the curve-fits for each material.
* Silicon Carbide:$k=2(10^{-8})T^{4}-2(10^{-5})T^{3}+0.0057T^{2}-0.7605T+87.873$
* Aluminum Nitride:$k=-1(10^{-9})T^{4}+1(10^{-7})T^{3}+0.0003T^{2}-0.0941T+28.634$

#### Part 3

Perform a nonlinear simulation using these temperature dependent properties to find the equilibrium temperature inside the silicon carbide region. State all simplifying assumptions that you have made.
The same code that was used in the previous problem is used here. However, the thermal conductivities are tabulated according to their temperature.
Shown in Figure 10 and 11 below are the contour plots of temperature and the vector sum of heat flux.
| Figure 10. Temperature of the nonlinear thermal conductivity model (K). | Figure 11. Flux of the nonlinear thermal conductivity model (W/m-K). |
There are many approaches to create the mesh for this problem. You will have to find the way the work best for you. These two codes listed are just two ways that the meshes can be manipulated.