Steady-state Heat Diffusion in a Slice of Reactor Dome.

The Problem.
 We want to solve the heat transfer problem on a 2-D slice of reactor dome constructed from two materials: steel and concrete. The self-adjoint form of the heat equation models the temperature distribution :
(K(x,y)* Tx )x + (K(x,y) * Ty )y = 0
where K(x,y) is the thermal conductivity of the materials. The PDE problem with domain and boundary conditions is shown to the right. Click here to see a physical description of the problem.
We will generate a uniform triangular mesh, and we select the bi-linear FEM discretizer with the Jacobi CG iterative linear system solver. We want to store the solution T and its derivatives Tx and Ty.

Use PDELab to ...

Why use PDELab?
 In this case study, we use PDELab's finite element methods to solve a 2D heat transfer problem. PDELab's bi-linear FEM discretizes a self-adjoint elliptic partial differential equation on an arbitrary 2-D domain for which a triangular or quadrilaterl finite element mesh can be generated. The system of linear equations produced by bi-linear FEM can be solved using any of the PDELab/ITPACK linear system solvers.

PDELab provides the following editors to assist in the process of specifying and solving the problem :

The purpose of this case study is to demonstrate how easily PDELab solves interesting problems.

Define the Problem
Step 1:When PDELab top level window appears, select 2-D and FEM and click on New File. The PDELab 2-D Finite Element Method Session appears. The Session toolkit is to the right of the Session window.

Step 2: Click on the Framework Editor in the Session toolkit.

Step 3: Choose PELLPACK from the Framework menu. The Framework Editor allows users to select the framework for defining the problem. In this case, selecting PELLPACK tells the symbolic processor that only linear, elliptic equations may be entered without special processing. The framework selection also signals the fortran code generator that the PELLPACK problem template must be generated. We select Self-Adjoint as the Equation Form.
Since PELLPACK modules can discretize only a single PDE, the Number of Equations is equal to 1, and the field cannot be edited. Click on PDE System. Enter the steady-state equation in the PDE System editor as shown here. Since we chose the self-adjoint form, we only enter the functions P, Q, and the right-hand-side. The function K(x,y) describing the thermal conductivity will be a user-defined [fortran] function which we will place in the subprograms segment of the session. Click on Done when finished.

Step 4: The equation has been specified, so click on Generate Framework. The Maxima symbolic processor is accessed to check for errors, and the GenCray code generator generates the appropriate PELLPACK problem specification. Any input errors are printed out in the Maxima window of the Framework Editor.

When the input has been successfully processed, select Quit from the File menu. The PDELab language program appears in the Session window. The generated program is shown below. The default rectangular domain and boundary conditions appear in the boundary segment. The default solution scheme includes the bi-linear fem discretizer and Jacobi CG linear system solver.

PDELab language `.e' file 

OPTIONS.
   clockwise = .true.
   time
EQUTION.
   ((K(X,Y))*ux)x + ((K(X,Y))*uy)y = 0
BOUNDARY.
   U=0  on x = 0
   U=0  on x = 1
   U=0  on y = 0
   U=0  on y = 1
DISCRETIZATION.
   Bi-Linear FEM
INDEXING.
   As Is
SOLUTION.
   Jacobi CG
END.

Step 5: We will now specify the domain and boundary conditions using the 2D BoundaryTool. Click on the 2D BoundaryTool.

Step 6: We could draw the domain using the 2D BoundaryTool, but decide instead to write our own parameterization of the boundary of the domain directly in the Session window. Click on the bar at the top of the session window : Click here to edit session. Change the boundary segment from 

boundary.
   u=0  on x = 0
   u=0  on x = 1
   u=0  on y = 0
   u=0  on y = 1
to 
boundary.
      (.7/.79)*u+ux=56./.79  on X=3.0, Y=.7-t for t=0.0 to .7
      uy=0  on X=3.0-t, Y=0.0 for t=0.0 to 2.0
      u=450 on X=cos(t), Y=sin(t) for t=0.0 to 1.5708
      ux=0  on X=0.0, Y=1.0+t for t=0.0 to .75
      (.7/30.62)*u +(x/sqrt(x*x+y*y))*ux+(y/sqrt(x*x+y*y))*uy =56./30.62 &
            on X=1.75*cos(pi/2.-t),  Y=1.75*sin(pi/2.-t)  &
              for t=0.0 to 1.1592795
      (.7/.79)*u + uy = 56./.79  on X=1.6039015+t,Y=.7 &
              for t=0.0 to 1.3960985

Step 7: Click on the 2D BoundaryTool. The parameterized boundary we just defined is dynamically loaded into PDELab (assuming you had no typing errors in the boundary definition). If you click on Set Conditions, you will see that the conditions have been set according to the specification entered in the session. The boundary conditions specify which sides of the reactor are insulated, which side is subjected to constant temperature, and which side has a heat flux condition between the interior of the domain and the surrounding air. Click on Cancel to dismiss the Boundary Conditions Editor.

Step 8: The final step in the problem specification is the definition of the function describing the thermal conductivity of the materials in the dome. Click on the bar above the session where it states Click here to edit session. Use the scroll bar to find the end of the .e file. Just before the end. segment, type in the following 

SUBPROGRAMS.  +both
*
* material properties : thermal conductivity for steel and concrete on
* which make up the reactor dome
*
          real function k(x,y)
          k=0.
          if (x.le.1..and.y.le.sqrt(3.0625-x**2)
     a          .and.y.ge.sqrt(1.-x**2)) k=30.6
          if (x.le.1.75.and.y.le.sqrt(3.0625-x**2)
     a       .and.x.gt.1.) k=30.62
              if (x.lt.1.75.and.x.ge.1.6039015.and.
     a        y.gt.sqrt(3.0625-x**2)) k=.79
          if (x.ge.1.75
     a           .and.x.le.3..and.y.le..7) k=.79
          return
          end
Click on the edit bar Click here to edit session when you have finished typing. The problem specification is complete.

Specify the Solution
Step 9: We will now generate the mesh. Do not close the BoundaryTool; it must be up when the Mesh Generator is invoked. Click on the MeshGenerator. The domain boundary appears in the window with the piece numbers identified.

Step 10: The Triangular mesh is the default selection. Enter 0.03 as the Edge Constraint, and use the default Element Constraint d < edge < 2d indicates that all generated element edge lengths will be between 0.06 and 2*0.06. Second, 30 deg < angle < 120 deg indicates that the angles of each mesh element will be no less 30 deg and no more than 120 deg, These two controls create a very uniform triangular mesh). Now click on Generate Mesh. The mesh to the right is generated and displayed. Finally, click on Save. A Filename Selection dialog requests a name for the file where the mesh data will be stored. Enter a filename, and click on Continue. The mesh data is written to the file, and a new mesh segment is written to the Session. It is placed immediately after the `boundary' segment and indicates that a file containing the mesh now exists and will be used in the problem solution process. 

MESH.
   read fem from file &
       (filename='reactor.mesh', &
        i1ndim=2, i1mnpt=2031, i1melm=3834, i1mtyp=1, i1mnen=4, &
        i3nint=4, i3ntyp=1, meshtype=2, meshedge=0.03)
You can now close the Mesh Generator by clicking on Quit. Then close the BoundaryTool by clicking on Quit.

Step 11: The bi-linear fem discretization module and the Jacobi CG linear solver module are the default specifications for the 2-D FEM session, so it is not necessary to modify them. However, you can click on the SolnMenu and look at the parameters for the Jacobi CG module which may be set by the user.

The Module Specification editor with the Jacobi CG module selected and the value of the zeta parameter displayed is shown to the left. Any parameter that may be set by the user is listed here. You can view and/or modify the values for any of the parameters. Dismiss the Module Specification editor by clicking on Cancel.

Step 11: We need to specify what output we want the system to produce for us. Click on the Output Specification Editor. Click on the Plot Function menu. Select u from the menu. Then select ux and uy. Now click on Output Filename. The Filename Selection dialog request the name of the file where the solution values for u will be stored. Enter a name for the output file. This file is a PDELab format file that can be loaded into the OutputTool enviroment where it can be used to visualize contour plots, solution plots, animations, etc. Now click on Continue. Two new segments are written to the Session: a save segment listing the name of the file, and an output segment listing the requested contents. Click on Continue. The new segments are shown below. The computed values for u, ux, and uy will be written to a file in PDELab format. 

SAVE.
   set output file (outputfile =  `reactor.out')
OUTPUT. 
   output  mesh
   output function(u)
   output function(ux)
   output function(uy)
Step 12: You can save the PDELab language program to a `.e' file by clicking on the SaveAs button in the Session toolkit. A Filename dialog requests a name for your .e file.
Here is the final .e file.

Execute the Problem
Step 13: We are now ready to enter the ExecuteTool environment, where the PDELab language program will be processed. A fortran main program will be generated from the .e file. It will be compiled and linked with the PDELab module libraries and the standard PDELab I/O libraries. Click on the ExecuteTool button in the Session toolkit.

Step 14: The PDELab .e file is loaded automatically. Click on the Run button The trace of the process for compiling and linking the PDELab executable appears in the trace window. When the compiling and linking process is complete, the executable program is run and the solution and trance files are generated and placed in your user directory. Click on Quit to exit the ExecuteTool.

Visualize the Solution
Step 15: The solution output has been generated during execution. You should find the output file containing the PDELab format data in the location specified by the save segment. If no path was specified, this file is located in the execution directory. The file contains the solution and its derivatives. Click on the OutputTool in the Session toolkit. Select PDELab Solution in the File menu. Select the output file that was generated by your program when the OutputTool requests the file name. The Select Solution Box is displayed after the file is loaded. The selection box contains the solution u and its derivatives ux and uy.


The temperature distribtuion.

Change of temperature in the x direction.
Change of temperature in the y direction.
Step 16: Select u from the Select Solution Box. Then choose the Visual2D visualization package and press Run Tool. The solution u is shown. Select Quit to exit Visual2D. You can return to the OutputTool and from the Solutions Selection Panel under Solutions in the OutputTool menu bar, de-select u and select ux. In the same manner, use Visual2D to visualize ux. Select Quit to exit the tool. Similary, view uy with Visual2D.