Nonlinear finite elements/Homework 2/Solutions



Problem 1: Classification of PDEsEdit

Given:

 

SolutionEdit

Part 1Edit

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 2Edit

Write out the PDEs in elementary partial differential notation

 

Part 3Edit

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 ProblemsEdit

Given:

 

SolutionEdit

  • 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 MethodEdit

Given:

The residual over an element is

 

where

 

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

 

SolutionEdit

Part 1Edit

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

We have,

 

Hence the weighting functions are of the form

 

Part 2Edit

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

 

and

 

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

 

Hence, the weighting functions are

 

The product of these two quantities is

 

Therefore,

 

These equations can be written as

 

where

 

and

 

Part 3Edit

  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 TranseferEdit

Given:

 

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

 

SolutionEdit

Part 1Edit

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 2Edit

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:

  /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 3Edit

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