Nonlinear finite elements/Homework 2/Solutions

Problem 1: Classification of PDEs edit



Solution edit

Part 1 edit

What physical processes do the two PDEs model?

Equation (1) models the one dimensional heat problem where   represents temperature at   at time   and   represents the thermal diffusivity.

Equation (2) models the transverse vibrations of a beam. The variable   represents the displacement of the point on the beam corresponding to position   at time  .   is   where   is Young's modulus,   is area moment of inertia,   is mass density, and   is area of cross section.

Part 2 edit

Write out the PDEs in elementary partial differential notation


Part 3 edit

Determine whether the PDEs are elliptic, hyperbolic, or parabolic. Show how you arrived at your classification.

Equation (1) is parabolic for all values of  .

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


Hence, we have  ,  ,  ,  ,  ,  , and  . Therefore,


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

Problem 2: Models of Physical Problems edit



Solution edit

  • Heat transfer: : 
  • Flow through porous medium: : 
  • Flow through pipe: : 
  • Flow of viscous fluids: : 
  • Elastic cables: : 
  • Elastic bars: : 
  • Torsion of bars: : 
  • Electrostatics: : 

Problem 4: Least Squares Method edit


The residual over an element is




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


Solution edit

Part 1 edit

What is the weighting function used in the least squares approach?

We have,


Hence the weighting functions are of the form


Part 2 edit

Develop a least-squares finite element model for the problem.

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

We start with the approximate solution


The residual is


The derivative of the residual is


For simplicity, let us assume that   is constant, and  . Then, we have




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


Hence, the weighting functions are


The product of these two quantities is




These equations can be written as






Part 3 edit

  1. Discuss some functions   that can be used to approximate  

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


Another possiblity is set set of polynomials such as


Problem 5: Heat Transefer edit



where   is the temperature,   is the thermal conductivity, and   is the heat generation. The Boundary conditions are


Solution edit

Part 1 edit

Formulate the finite element model for this problem over one element

First step is to write the weighted-residual statement


Second step is to trade differentiation from   to  , using integration by parts. We obtain


Rewriting (3) in using the primary variable   and secondary variable   as described in the text book to obtain the weak form


From (4) we have the following


Part 2 edit

Use ANSYS to solve the problem of heat conduction in an insulated rod. Compare the solution with the exact solution. Try at least three refinements of the mesh to confirm that your solution converges.

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

  pi = 4*atan(1) !compute value for pi
  ET,1,LINK34!convection element
  ET,2,LINK32!conduction element
  R,1,1!AREA = 1
  MP,HF,1,25 !convectivity
  emax = 2
  length = 0.1
  *do,nnumber,1,emax+1 !begin loop to creat nodes
     n,emax+2,length,0!extra node for the convetion element
  type,2 !switch to conduction element
  TYPE,1 !switch to convection element
  e,emax+1,emax+2!creat zero length element



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 edit

Write you own code to solve the problem in part 2.

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)"]);