Nonlinear finite elements/Homework 2/Solutions



Problem 1: Classification of PDEs

Given:

$\begin{array}{ccc}\text{(1)}{u}_{,xx}& =& \alpha u_{,t}\\ \text{(2)}{u}_{,xxxx}& =& \alpha u_{,tt}\end{array}$

Solution

Part 1

Equation (1) models the one dimensional heat problem where $u(x,t)$ represents temperature at $x$ at time $t$ and $\alpha$ represents the thermal diffusivity.

Equation (2) models the transverse vibrations of a beam. The variable $u(x,t)$ represents the displacement of the point on the beam corresponding to position $x$ at time $t$. $\alpha$ is $-\frac{\rho A}{EI}$ where $E$ is Young's modulus, $I$ is area moment of inertia, $\rho$ is mass density, and $A$ is area of cross section.

Part 2

$\begin{array}{cc}{u}_{,xx}& =\alpha u_{,t}\equiv \frac{{\partial }^{2}u}{\partial x^2}=\alpha \frac{\partial u}{\partial t}\\ u_{,xxxx}& =\alpha u_{,tt}\equiv \frac{{\partial }^{4}u}{\partial x^4}=\alpha \frac{{\partial }^{2}u}{\partial t^2}\end{array}$

Part 3

Equation (1) is parabolic for all values of $\alpha$.

We can see why by writing it out in the form given in the notes on Partial differential equations.

$(1)\frac{{\partial }^{2}u}{\partial x^2}+(0)\frac{{\partial }^{2}u}{\partial x\partial t}+(0)\frac{{\partial }^{2}u}{\partial t^2}+(0)\frac{\partial u}{\partial x}+(-\alpha )\frac{\partial u}{\partial t}+(0)u=0$

Hence, we have $a=1$, $b=0$, $c=0$, $d=0$, $e=-\alpha$, $f=0$, and $g=0$. Therefore,

$b^2-4ac=0-(1)(0)=0\Rightarrow \text{parabolic}$

Equation (2) is hyperbolic for positive values of $\alpha$, parabolic when $\alpha$ is zero, and elliptic when $\alpha$ is less than zero.

Problem 2: Models of Physical Problems

Given:

$\frac{d}{dx}(a(x)\frac{du(x)}{dx})+c(x)u(x)=f(x)\text{for }x\in (O,L)$

Solution

• Heat transfer: :$kA\frac{d^{2}T}{d{x}^2}+ap\beta T=f$
• Flow through porous medium: :$\mu \frac{d^2\varphi }{d{x}^2}=f$
• Flow through pipe: :$\frac{1}{R}\frac{d^{2}P}{d{x}^2}=0$
• Flow of viscous fluids: :$\mu \frac{d^2{v}_z}{d{x}^2}=-\frac{dP}{dx}$
• Elastic cables: :$T\frac{d^{2}u}{d{x}^2}=f$
• Elastic bars: :$EA\frac{d^{2}u}{d{x}^2}=f$
• Torsion of bars: :$GJ\frac{d^2\theta }{d{x}^2}=0$
• Electrostatics: :$\epsilon \frac{d^2\varphi }{d{x}^2}=\rho$

Problem 4: Least Squares Method

Given:

The residual over an element is

$R^e(x,c_1^e,c_2^e,\dots ,c_n^e):=\frac{d}{dx}\left(a\frac{d{u}_h^e}{dx}\right)+c{u}_h^e-f(x)$

where

$u_h^e(x)=\sum _{j=1}^n{c}_j^e{\phi }_j^e(x).$

In the least squares approach, the residual is minimized in the following sense

$\frac{\partial }{\partial c_i}\underset{x_a}{\overset{x_b}{\int }}[R^e(x,c_1^e,c_2^e,\dots ,c_n^e){]}^{2}dx=0,i=1,2,\dots ,n.$

Solution

Part 1

We have,

$\frac{\partial }{\partial c_i}\underset{x_a}{\overset{x_b}{\int }}[R^e{]}^{2}dx=\underset{x_a}{\overset{x_b}{\int }}2R^e\frac{\partial R^e}{\partial c_i}dx=0.$

Hence the weighting functions are of the form

$w_i(x)=\frac{\partial R^e}{\partial c_i}$

Part 2

To make the typing of the following easier, I have gotten rid of the subscript $h$ and the superscript $e$. Please note that these are implicit in what follows.

We start with the approximate solution

$u(x)=\sum _{j=1}^n{c}_j{\phi }_j(x).$

The residual is

$R=\frac{d}{dx}\left[a\left(\sum _j{c}_j\frac{d{\phi }_j}{dx}\right)\right]+c\left(\sum _j{c}_j{\phi }_j\right)-f(x).$

The derivative of the residual is

$\frac{\partial R}{\partial c_i}=\frac{d}{dx}\left[a\left(\sum _j\frac{\partial c_j}{\partial c_i}\frac{d{\phi }_j}{dx}\right)\right]+c\left(\sum _j\frac{\partial c_j}{\partial c_i}{\phi }_j\right)-f(x).$

For simplicity, let us assume that $a$ is constant, and $c=0$. Then, we have

$R=a\left(\sum _j{c}_j\frac{d^2{\phi }_j}{d{x}^2}\right)-f,$

and

$\frac{\partial R}{\partial c_i}=a\left(\sum _j\frac{\partial c_j}{\partial c_i}\frac{d^2{\phi }_j}{d{x}^2}\right).$

Now, the coefficients $c_i$ are independent. Hence, their derivatives are

$\frac{\partial c_j}{\partial c_i}=\{\begin{array}{cc}0& \text{for}i\ne j\\ 1& \text{for}i=j\end{array}$

Hence, the weighting functions are

$\frac{\partial R}{\partial c_i}=a\frac{d^2{\phi }_i}{d{x}^2}.$

The product of these two quantities is

$R\frac{\partial R}{\partial c_i}=(a\sum _j{c}_j\frac{d^2{\phi }_j}{d{x}^2}-f)\left(a\frac{d^2{\phi }_i}{d{x}^2}\right)=\sum _j\left[\left(a^2\frac{d^2{\phi }_i}{d{x}^2}\frac{d^2{\phi }_j}{d{x}^2}\right)c_j\right]-af\frac{d^2{\phi }_i}{d{x}^2}.$

Therefore,

$\underset{x_a}{\overset{x_b}{\int }}R\frac{\partial R}{\partial c_i}dx=0⟹\sum _j\left[\left(\underset{x_a}{\overset{x_b}{\int }}{a}^2\frac{d^2{\phi }_i}{d{x}^2}\frac{d^2{\phi }_j}{d{x}^2}dx\right)c_j\right]-\underset{x_a}{\overset{x_b}{\int }}af\frac{d^2{\phi }_i}{d{x}^2}dx=0$

These equations can be written as

$K_{ij}{c}_j=f_i$

where

$K_{ij}=\underset{x_a}{\overset{x_b}{\int }}{a}^2\frac{d^2{\phi }_i}{d{x}^2}\frac{d^2{\phi }_j}{d{x}^2}dx$

and

$f_i=\underset{x_a}{\overset{x_b}{\int }}af\frac{d^2{\phi }_i}{d{x}^2}dx.$

Part 3

1. Template:Font

For the least squares method to be a true "finite element" type of method we have to use finite element shape functions. For instance, if each element has two nodes we could choose

${\phi }_1=\frac{x_b-x}{h}\text{and}{\phi }_2=\frac{x-x_a}{h}.$

Another possiblity is set set of polynomials such as

$\phi =\{1,x,x^2,x^3\}.$

Problem 5: Heat Transefer

Given:

$-\frac{d}{dx}\left(k\frac{dT}{dx}\right)=q\text{ for }x\in (0,L)$

where $T$ is the temperature, $k$ is the thermal conductivity, and $q$ is the heat generation. The Boundary conditions are

$\begin{array}{ccc}T(0)& =& T_0\\ k\frac{dT}{dx}+\beta (T-T_{\mathrm{\infty }})+\hat{q}{|}_{x=L}& =& 0\end{array}$

Solution

Part 1

First step is to write the weighted-residual statement

$0=\underset{x_a}{\overset{x_b}{\int }}{w}_i^e[-\frac{d}{dx}\left(k\frac{d{T}_h^e}{dx}\right)-q]dx$

Second step is to trade differentiation from $T_h^e$ to $w_i^e$, using integration by parts. We obtain

$0=\underset{x_a}{\overset{x_b}{\int }}(k\frac{w_i^e}{dx}\frac{d{T}_h^e}{dx}-w_i^{e}q)dx-{\left[w_i^{e}k\frac{d{T}_h^e}{dx}\right]}_{x_a}^{x_b}\text{(3)}$

Rewriting (3) in using the primary variable $T$ and secondary variable $k\frac{dT}{dx}\equiv Q$ as described in the text book to obtain the weak form

$0=\underset{x_a}{\overset{x_b}{\int }}(k\frac{w_i^e}{dx}\frac{d{T}_h^e}{dx}-w_i^{e}q)dx-w_i^e(x_a)Q_a^e-w_i^e(x_b)Q_b^e\text{(4)}$

From (4) we have the following

$\begin{array}{ccc}{K}_{ij}^e& =& \underset{x_a}{\overset{x_b}{\int }}k\frac{N_i^e}{dx}\frac{d{N}_j^e}{dx}dx\\ f_i^e& =& \underset{x_a}{\overset{x_b}{\int }}q{N}_i^{e}dx\\ K_{ij}^e{T}_j^e& =& f_i^e+N_i^e(x_a)Q_a^e+N_i^e(x_b)Q_b^e\end{array}$

Part 2

The following inputs are used to solve for the solution of this problem:

  /prep7
  pi = 4*atan(1) !compute value for pi
  ET,1,LINK34!convection element
  ET,2,LINK32!conduction element
  R,1,1!AREA = 1
  MP,KXX,1,0.01!conductivity
  MP,HF,1,25 !convectivity
  emax = 2
  length = 0.1
  *do,nnumber,1,emax+1 !begin loop to creat nodes
     n,nnumber,length/emax*(nnumber-1),0
     n,emax+2,length,0!extra node for the convetion element
  *enddo
  type,2 !switch to conduction element
  *do,enumber,1,emax
     e,enumber,enumber+1
  *enddo
  TYPE,1 !switch to convection element
  e,emax+1,emax+2!creat zero length element
  D,1,TEMP,50
  D,emax+2,TEMP,5
  FINISH

  /solu
  solve
  fini

  /post1
  /output,p5,txt
  prnsol,temp
  prnld,heat
  /output,,
  fini

ANSYS gives the following results

Node	Temperature
1	50.000
2	27.590
3	5.1793

Node	Heat Flux
1	-4.4821
4	4.4821

Part 3

The following Maple code can be used to solve the problem.

> restart;
> with(linalg):
> L:=0.1:kxx:=0.01:beta:=25:T0:=50:Tinf:=5:
> # Number of elements
> element:=10:
> # Division length
> ndiv:=L/element:
> # Define node location
> for i from 1 to element+1 do
>   x[i]:=ndiv*(i-1);
> od:
> # Define shape function
> for j from 1 to element do
>   N[j][1]:=(x[j+1]-x)/ndiv;
>   N[j][2]:=(x-x[j])/ndiv; 
> od:
> # Create element stiffness matrix
> for k from 1 to element do
>   for i from 1 to 2 do
>     for j from 1 to 2 do
>       Kde[k][i,j]:=int(kxx*diff(N[k][i],x)*diff(N[k][j],x),x=x[k]..x[k+1]); 
>     od;
>   od;
> od;
> Kde[0][2,2]:=0:Kde[element+1][1,1]:=0:
> # Create global matrix for conduction Kdg and convection Khg
> Kdg:=Matrix(1..element+1,1..element+1,shape=symmetric):
> Khg:=Matrix(1..element+1,1..element+1,shape=symmetric):
> for k from 1 to element+1 do
>   Kdg[k,k]:=Kde[k-1][2,2]+Kde[k][1,1];
>   if (k <= element) then
>     Kdg[k,k+1]:=Kde[k][1,2];
>   end if;
> od:
> Khg[element+1,element+1]:=beta:
> # Create global matrix Kg = Kdg + Khg
> Kg:=matadd(Kdg,Khg):
> # Create T and F vectors
> Tvect:=vector(element+1):Fvect:=vector(element+1):
> for i from 1 to element do
>   Tvect[i+1]:=T[i+1]:
>   Fvect[i]:=f[i]:
> od:
> Tvect[1]:=T0:
> Fvect[element+1]:=Tinf*beta:
> # Solve the system Ku=f
> Ku:=multiply(Kg,Tvect):
> # Create new K matrix without first row and column from Kg
> newK:=Matrix(1..element,1..element):
> for i from 1 to element do
>   for j from 1 to element do
>     newK[i,j]:=coeff(Ku[i+1],Tvect[j+1]);
>   od;
> od;
> # Create new f matrix without f1
> newf:=vector(element,0):
> newf[1]:=-Ku[2]+coeff(Ku[2],T[2])*T[2]+coeff(Ku[2],T[3])*T[3]:
> newf[element]:=Tinf*beta:
> # Solve the system newK*T=newf
> solution:=multiply(inverse(newK),newf);
> exact := 50-448.2071713*x;
> with(plots):
> # Warning, the name changecoords has been redefined
> data := [[ x[n+1], solution[n]] <math>n=1..element]:
> FE:=plot(data, x=0..L, style=point,symbol=circle,legend="FE"):
> ELAS:=plot(exact,x=0..L,legend="exact",color=black):
> display({FE,ELAS},title="T(x) vs x",labels=["x","T(x)"]);