FEATool  v1.7
Finite Element Analysis Toolbox
 All Files Functions Pages
Tutorials

A selection of tutorial models and examples are presented in this section.


Heat Transfer - Shrink Fitting of an Assembly

FEATool supports modeling heat transfer through both conduction, that is heat transported by a diffusion process, and also convection, which is heat transported through a fluid through convection by a velocity field. The heat transfer physics mode supports both these processes, and defines the following equation

\[ \rho C_p\frac{\partial T}{\partial t} + \nabla\cdot(-k\nabla c) = Q - \rho C_p\mathbf{u}\cdot\nabla T \]

where $\rho$ is the density, $C_p$ the heat capacity, $k$ is the thermal conductivity, $Q$ heat source term, and $\mathbf{u}$ a vector valued convective velocity field.

ex_ht5_geom_50.svg

This example models heat conduction in the form of transient cooling for shrink fitting of a two part assembly. A tungsten rod heated to 84 C is inserted into a -10 C chilled steel frame part. The time when the maximum temperature has cooled to 70 C should be determined. The assembly is cooled due to convection through a surrounding medium kept at Tinf = 17 C and a heat transfer coefficient of h = 750 W/m2K. The surrounding cooling medium is not modeled directly, thus the convective term is omitted, but the effects are incorporated into the model by the use of natural convection boundary conditions [1].


Heat Transfer - Shrink Fitting of an Assembly using the GUI

This section describes how to set up and solve the thermal shrink fitting example with the FEATool graphical user interface (GUI).

  1. Start Octave or Matlab, and if you have not run the installation script (which automatically adds the FEATool directory paths at startup) then change your working directory to where your FEATool installation is, for example
    cd C:\featool
    
  2. In the command window type

    featool
    

    to start the graphical user interface (GUI).

  3. Click on the New Model button in the upper horizontal toolbar to clear all data and start defining a new model.

    ex_poi1_01_50.png


  4. In the opened New Model dialog box, click on the 2D radio button in the Select Space Dimensions frame, and select Heat Transfer from the Select Physics drop-down list. Leave the space dimension and dependent variable names to their default values. Finish and close the dialog box by clicking on the OK button.

    ex_ht5_04_50.png


  5. To create the first steel part, use the Create square/rectangle , Create circle/ellipse tools to create one rectangle and three circles. The Inspect/edit selected geometry object toolbar button allows access to edit geometry object dimensions and properties. Give the rectangle dimensions so that it spans it spans between 0 and 0.11 in the x-direction, and 0 and 0.12 in the y-direction. Similarly, set the circles radii to 0.015, 0.035, 0.025, and centers to (0.065, 0), (0.11, 0.12), and (0, 0.06).

    ex_ht5_05_50.png


  6. Generate a composite object for the part by using Combine Objects... from the Geometry menu. By using the formula R1-E1-E2-E3 all circles will be subtracted from the rectangle.

    ex_ht5_06_50.png


  7. Create another rectangle with dimensions (0.065, 0.16, 0.05, 0.07) and circle with radius 0.01 and center (0.065, 0.06) for the second part.

    ex_ht5_07_50.png


  8. In this case the objects should be joined. Click on them so that they are highlighted in red. Alternatively, they can be selected by holding the Ctrl key while clicking on the labels R2 and E4 in the Selection list box (or simply press Ctrl+a to select all objects). When the geometry objects are selected, press the Add geometry objects button to generate the combined shape.

    ex_ht5_08_50.png


  9. Switch to Grid mode by click on the corresponding mode toolbar button. Set the target grid size by entering 0.003 in the Grid Size edit field above the grid toolbar buttons. Then click on the Generate unstructured grid button to call the automatic grid generation function.

    ex_ht5_09_50.png


  10. Press the mode button in the Mode toolbar to switch from grid mode to physics and equation/subdomain specification mode. In the Equation Settings dialog box, first select the subdomain that corresponds to the steel part (here 2), then set the density $\rho$ to 7500, heat capacity $C_p$ to 470, and thermal conductivity $k$ to 44 in the corresponding edit fields. Also set the initial temperature $T_0$ here to -10.

    ex_ht5_10_50.png


  11. Now set the material parameters for the Tungsten part by selecting the other subdomains (1 and 3). Here set the density $\rho$ to 19000, heat capacity $C_p$ to 134, thermal conductivity $k$ to 163, and the initial temperature $T_0$ to 84. Press OK to finish.

    ex_ht5_11_50.png


  12. Switch to boundary condition specification mode by clicking on the mode button. In the Boundary Settings dialog box, select the Heat flux boundary condition for all the boundaries. Enter k_ht for the convective transfer coefficient $h$ (this will use the actual value for the subdomain parameter k so that we don't have to enter the values again), and 17 for the surrounding reference temperature $T_{inf}$.

    ex_ht5_12_50.png


  13. Now that the problem is specified, press the mode button to switch to solve mode. Since this is a time dependent study, open the solver settings and select the Time-Dependent - Monolithic (Fully Coupled) solver. Set the Time step to 0.25, Simulation time to 16, Time stopping criteria to 0, also Maximum non-linear iterations to 2. Then press Apply and Solve to start the solution process.

    ex_ht5_13_50.png


  14. After the problem has been solved the temperature at the final time will be shown. To find the time where the maximum temperature is 70 one can study the solution at different times by using the corresponding option in the Postprocessing Settings dialog box. Note that both the colorbar and Limits field will show the minimum and maximum surface plot value.

    ex_ht5_14_50.png


  15. To see this more clearly we can process the simulation data on the command line by using the Export > FEA Problem Struct To Main Workspace from the File menu. Then switch to the main command window and use the following script commands to see that T = 70 is reached at about t = 11.

    tlist = fea.sol.t;
    for i=1:numel(tlist)
      T_min(i) = min(fea.sol.u(:,i));
      T_max(i) = max(fea.sol.u(:,i));
    end
    ix = find(T_max<70);
    i1 = ix(1);
    i2 = i1 - 1;
    s  = ( T_max(i2) - 70 )/( T_max(i2) - T_max(i1) );
    t_70 = tlist(i2) + s*( tlist(i1) - tlist(i2) )
    u_70 = fea.sol.u(:,i2) + ...
           s*( fea.sol.u(:,i1) - fea.sol.u(:,i2) )
    
    figure
    plot( tlist, T_min, 'b-' )
    hold on
    plot( tlist, T_max, 'r-' )
    grid on
    title('Maximum and minimum temperatures')
    ylabel('Temperature [C]')
    xlabel('Time [s]')
    



    ex_ht5_15_50.png



Heat Transfer - Shrink Fitting of an Assembly using the CLI

How to set up and solve the thermal shrink fitting problem on the command line interface is illustrated in the ex_heattransfer5 script file which can be found in the FEATool examples directory.


References

[1] Krysl P. A Pragmatic Introduction to the Finite Element Method for Thermal and Stress Analysis. Pressure Cooker Press, USA, 2005.


Fluid Dynamics - Flow Around a Cylinder

A well known benchmark, test, and validation problem suite for incompressible fluid flows are the DFG cylinder benchmark problems [2]. Although it is not possible to derive analytical solutions to these test cases, accurate numerical solutions to benchmark reference quantities have been established for a number of configurations [3], [4].

The test configuration used in the following places a solid cylinder centered at (0.2, 0.2) with diameter $D=0.1$ in a $L=2.2$ by $H=0.41$ rectangular channel. The fluid density $\rho$ is taken as 1 and the viscosity $\mu$ is 0.001. A fully developed parabolic velocity profile is prescribed at the inlet, $\mathbf{u}_{inflow} = \mathbf{u}(0, y) = 4U_{max}H^{-2}(y(H-y), 0)$, with the maximum velocity $U_{max}=0.3$. This results in a mean velocity $U_{mean}=0.2$ and a Reynolds number $Re=\rho U_{mean}D/\mu=20$ and thus the flow field will be laminar.

ex_ns1_geom.svg

As the fluid is considered incompressible the problem is governed by the Navier-Stokes equations. That is,

\[ \left\{\begin{aligned} \rho\left( \frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\right) - \nabla\cdot(\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^T)) + p &= \mathbf{F} \\ \nabla\cdot\mathbf{u} &= 0 \end{aligned}\right. \]

where in this case the time dependent term can be neglected. The benchmark quantities that should be computed include the pressure difference between the front and rear of the cylinder $\Delta p = p(0.15,0.2)−-p(0.25,0.2)$, and the coefficients of drag $c_d$ and lift $c_l$, defined as

\[ c_d = \frac{2F_d}{\rho U_{mean}^2 D},\qquad c_l = \frac{2F_l}{\rho U_{mean}^2 D} \]

The drag and lift forces, $F_d$ and $F_l$, can be computed as

\[ F_d = \phantom{-}\int_S \left( \mu\frac{\partial\mathbf{u}_{\tau}(t)}{\partial n}n_y-pn_x \right) dS,\qquad F_l = -\int_S \left( \mu\frac{\partial\mathbf{u}_{\tau}(t)}{\partial n}n_x+pn_y \right) dS \]

where $\mathbf{u}_{\tau}$ is the velocity in the tangential direction $\tau = (n_y, -n_x, 0)^T$.


Flow Around a Cylinder using the GUI

This section describes how to set up and solve the Flow Around a Cylinder example with the FEATool graphical user interface (GUI).

  1. Start Octave or Matlab, and if you have not run the installation script (which automatically adds the FEATool directory paths at startup) then change your working directory to where your FEATool installation is, for example
    cd C:\featool
    
  2. In the command window type

    featool
    

    to start the graphical user interface (GUI).

  3. Click on the New Model button in the upper horizontal toolbar to clear all data and start defining a new model.

    ex_poi1_01_50.png


  4. In the opened New Model dialog box, click on the 2D radio button in the Select Space Dimensions frame, and select Navier-Stokes Equations from the Select Physics drop-down list. Leave the space dimension and dependent variable names to their default values. Finish and close the dialog box by clicking on the OK button.

    ex_ns1_02_50.png


  5. To create the outer rectangle, first click on the Create square/rectangle button in the left hand side Tools toolbar frame.

    ex_ns1_03_50.png


  6. Then left click in the main plot axes window, hold the mouse button, and move the mouse pointer to show a rectangle with red outlines.

    ex_ns1_04_50.png


  7. Release the mouse button to finalize and create a solid geometry object. The object properties must now be edited to set the correct size and position of the rectangle. To do this, click on the rectangle R1 to select it and highlight it in red. Then click on the Inspect/edit selected geometry object toolbar button

    ex_ns1_05_50.png


  8. In the opened Edit Geometry Object dialog box, edit the min and max point coordinates to specify a rectangle with length 2.2 and height 0.41. Finish editing the geometry object and close the dialog box by clicking OK.

    ex_ns1_06_50.png


  9. To create the inner circle, first left click on the Create circle/ellipse button in the left hand side Tools toolbar frame.

    ex_ns1_07_50.png


  10. Then click and hold the left mouse button anywhere in the main plot axes window, move the mouse pointer to show red outlines of a circle or ellipse. Release the button to finalize and create a solid geometry object.

    ex_ns1_08_50.png


  11. The object properties of the E1 must also be changed to make a circle with radius 0.05 centered at (0.2, 0.2). To do this, click on the E1 to select it which also highlights it in red. Then again click on the Inspect/edit selected geometry object toolbar button

    ex_ns1_09_50.png


  12. In the opened Edit Geometry Object dialog box change the center coordinates to 0.2 0.2, and the x and y radius 0.05 in the corresponding edit fields. Finish editing E1 and close the dialog box by clicking OK.

    ex_ns1_10_50.png


  13. To subtract the circle from the rectangle first select both geometry objects by clicking on them so both are highlighted in red. Alternatively, if the circle is obscured by the rectangle they can be selected by holding the Ctrl key while clicking on the labels R1 and E1 in the Selection list box (or simply press Ctrl+a to select all objects). When the geometry objects are selected, press the Subtract geometry objects button to generate the final geometry shape.

    ex_ns1_10b_50.png


  14. Switch to Grid mode by click on the corresponding mode toolbar button.

    ex_ns1_10c_50.png


  15. Set a target grid size by entering 0.025 in the Grid Size edit field above the grid toolbar buttons. Then click on the Generate unstructured grid button to call the automatic grid generation function.

    ex_ns1_11_50.png


  16. Press the mode button in the Mode toolbar to switch from grid mode to physics and equation/subdomain specification mode.

    ex_ns1_12_50.png


  17. In the Equation Settings dialog box that automatically opens, set the density $\rho$ to 1 and viscosity $\mu$ to 0.001 in the corresponding edit fields. The other coefficients can be left to their zero default values. Press OK to finish and close the dialog box.

    ex_ns1_13_50.png


  18. Switch to boundary condition specification mode by clicking on the mode button

    ex_ns1_14_50.png


  19. In the Boundary Settings dialog box, select all boundaries except the right and left ones (here boundaries 1, 3, and 5-8) in the left hand side Boundaries list box and choose the Wall/no-slip boundary conditions from the drop-down list.

    ex_ns1_15_50.png


  20. Now select the left inflow boundary (here number 4) in the left hand side Boundaries list box and choose the Inlet/velocity boundary condition from the drop-down list. To specify a parabolic velocity profile enter the expression 4*0.3/0.41^2*y*(0.41-y) in the edit field for the x-velocity coefficient $u_o$.

    ex_ns1_16_50.png


  21. Select the right outflow boundary (here number 2) in the left hand side Boundaries list box and choose the Neutral outflow/stress boundary condition from the drop-down list (alternatively one can prescribe a pressure $p_o$ with the Outflow/pressure condition). Finish by clicking the OK button.

    ex_ns1_17_50.png


  22. Now that the problem is specified, press the mode button to switch to solve mode. Then press the button with an equals sign in the Tools toolbar frame to call the solver with the default solver settings.

    ex_ns1_18_50.png


  23. After the problem has been solved FEATool will automatically switch to postprocessing mode and display the computed solution. The default postprocessing visualization shows the magnitude of the velocity field. Since the solution does not look quite smooth and resolved it might be worthwhile to refine the grid. Press the mode button in the Mode toolbar to switch back to grid mode again.

    ex_ns1_19_50.png


  24. Although we could use the uniform grid refinement option, it is usually better to use the automatic grid generation routine with a smaller grid size for better grid quality and approximation of curved boundaries. Set a new target grid size by entering 0.015 in the Grid Size edit field above the grid toolbar buttons. Then click on the Generate unstructured grid button .

    ex_ns1_20_50.png


  25. Now go back to Solve mode and solve the problem again with the refined grid.

    ex_ns1_21_50.png


  26. Now that the solution looks better we can calculate the drag and lift coefficients. To do this select Boundary Integration... from the Post menu.

    ex_ns1_22_50.png


  27. In the Boundary Integration dialog box, select boundaries which make up the circle (5-8) in the left hand side Boundaries list box. Then enter the expression for the drag coefficient 2*(nx*p+miu_ns*(-2*nx*ux-ny*(uy+vx)))/(1*0.2^2*0.1) in the Integration Expression edit field. Press the Apply button to show the result in the lower Integration Result frame and in the FEATool terminal window.

    ex_ns1_23_50.png


  28. The lift coefficient can similarly be computed by evaluating the expression 2*(T_y)/(1*0.2^2*0.1) where Ty is the total force, y-component expression (which can be selected in the drop down box), defined by ny*p+miu_ns*(-nx*(vx+uy)-2*ny*vy). The magnitudes of the computed drag and lift coefficients are here, 5.4392 and 0.0137, which can be compared with the benchmark reference solutions, that is 5.5795 and 0.010619. The deviation from the reference values suggest that it might be desirable to recompute a solution with an even finer grid, possibly refined around the circle. See for example ex_navierstokes3 where a very adapted benchmark grid is constructed manually.


Flow Around a Cylinder using the CLI

The process to set up and solve the flow around cylinder test problem on the command line interface is illustrated in the ex_navierstokes3 script file which can be found in the examples directory.


Instationary Flow Around a Cylinder

In addition to the stationary test case described above. An instationary benchmark test case is also available (ex_navierstokes6). This test case uses the same geometry but instead applies an inflow condition that varies with time (u_inflow = 6*sin(pi*t/8)*(y*(0.41-y))/0.41^2) so that 0<Re(t)<100 [5]. Computations with FEATool show that the drag and lift coefficients, and pressure difference between front and rear of the cylinder agrees very well with the reference values [6].


References

[2] Schäfer M, Turek S. Benchmark computations of laminar flow around cylinder. Flow Simulation with High–Performance Computers II (Notes on Numerical Fluid Mechanics), 52, 547–566, Vieweg, 1996.

[3] Nabh G. On higher order methods for the stationary incompressible Navier-Stokes equations. PhD Thesis 1998; Universität Heidelberg. Preprint 42/98, 1998.

[4] John V, Matthies G. Higher-order finite element discretizations in a benchmark problem for incompressible flows. International Journal for Numerical Methods in Fluids 2001; 37(8):885–903, DOI: 10.1002/fld.195.

[5] John V. Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder. International Journal for Numerical Methods in Fluids 2004; 44:777-788.

[6] John V, Rang J. Adaptive time step control for the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Engrg. 199 (2010) 514–524.


Fluid Dynamics - Axisymmetric Fluid Flow

The Equation Editing Example - Axisymmetric Fluid Flow example in the quickstart guide shows how to edit predefined equations in the physics modes to model rotationally symmetric or axisymmetric flow in a narrowing pipe.


Fluid Dynamics - Flow Over a Backwards Facing Step

Flow over a backwards facing step is another classic computational fluid dynamics test problem which has in various forms been used quite extensively for validation of simulation codes. The test problem essentially consists of studying how a fully developed flow profile reacts to a sudden expansion or step in a channel. The expansion will cause a break in the flow and a recirculation or separation zone will form. To measure and compare results the resulting length of the recirculation or separation zone is most often used.

The stationary incompressible Navier-Stokes equations are to be applied corresponding to a Reynold number of $Re = 389$. The size domain is chosen according to the data below. The inlet velocity function is given as $u_{in}(0, y) = 4u_{max}h^{-2}(y(h-y), 0)$ where $h$ is the channel height, and maximum velocity $ u_{max} = 1$ No-slip zero velocity conditions are applied to all walls and a suitable outflow condition must also be applied. The reference recirculation zone length found in reference [7] and [8] is estimated to be 7.93 length units.


Flow over a Backwards Facing Step with the external FeatFlow CFD Solver

This section describes how to set up and solve the Flow over a Backwards Facing Step example with the FEATool graphical user interface (GUI) and using the FeatFlow external CFD Solver.

  1. Start Octave or Matlab, and if you have not run the installation script (which automatically adds the FEATool directory paths at startup) then change your working directory to where your FEATool installation is, for example
    cd C:\featool
    
  2. In the command window type

    featool
    

    to start the graphical user interface (GUI).

  3. Click on the New Model button in the upper horizontal toolbar to clear all data and start defining a new model.

    ex_poi1_01_50.png


  4. In the opened New Model dialog box, click on the 2D radio button in the Select Space Dimensions frame, and select Navier-Stokes Equations from the Select Physics drop-down list. Leave the space dimension and dependent variable names to their default values. Finish and close the dialog box by clicking on the OK button.

    ex_ns1_02_50.png


  5. To create the outer rectangle, first click on the Create square/rectangle button in the left hand side Tools toolbar frame.

    ex_ns1_03_50.png


  6. Then left click in the main plot axes window, hold the mouse button, and move the mouse pointer to show a rectangle with red outlines.

    ex_ns1_04_50.png


  7. Release the mouse button to finalize and create a solid geometry object. The object properties must now be edited to set the correct size and position of the rectangle. To do this, click on the rectangle R1 to select it and highlight it in red. Then click on the Inspect/edit selected geometry object toolbar button

    ex_ns1_05_50.png


  8. In the opened Edit Geometry Object dialog box, edit the min and max point coordinates to specify a rectangle with xmin -1.9802, xmax 7.9802, ymin 0, and ymax 1. Finish editing the geometry object and close the dialog box by clicking OK.

  9. Repeat the steps to create another smaller rectangle with xmin -1.9802, xmax 0, ymin 0, and ymax 0.4851. Subtract the rectangles from each other by selecting both rectangles and clicking on the button.

    ex_bfs_01_50.png


  10. Switch to Grid mode by click on the corresponding mode toolbar button.

  11. Set a target grid size by entering 0.5 in the Grid Size edit field above the grid toolbar buttons. Then click on the Generate unstructured grid button to call the automatic grid generation function.

    ex_bfs_02_50.png


  12. The FeatFlow CFD solver requires quadrilateral in 2D or hexahedral grids in 3D. Select the Convert Grid Cells from the Grid menu and apply grid smoothing by using the corresponding Grid Smoothing menu options.

    ex_bfs_03_50.png


  13. Press the mode button in the Mode toolbar to switch from grid mode to physics and equation/subdomain specification mode. In the Equation Settings dialog box that automatically opens, set the density $\rho$ to 1 and viscosity $\mu$ to 2/3/389 in the corresponding edit fields. The other coefficients can be left to their zero default values. Press OK to finish and close the dialog box.

    ex_bfs_04_50.png


  14. Click on the Model Constrains and Expressions and enter the constants and expressions shown below.

    ex_bfs_05_50.png


  15. Switch to boundary condition specification mode by clicking on the mode button. Select the left inflow boundary (here number 3) in the left hand side Boundaries list box and choose the Inlet/velocity boundary condition from the drop-down list. To specify a parabolic velocity profile enter the expression 4*u_max*(y-hstep)*(1-y)/h_inlet^2 in the edit field for the x-velocity coefficient $u_o$.

    ex_bfs_06_50.png


  16. Select the right outflow boundary (here number 1) in the left hand side Boundaries list box and choose the Neutral outflow/stress boundary condition from the drop-down list. Finish by clicking the OK button.

    ex_bfs_07_50.png


  17. Now that the problem is specified, press the mode button to switch to solve mode. Then press the button labeled FeatFlow to open the FeatFlow CFD solver settings and control panel.

  18. In the FeatFlow solver settings dialog box set the Grid refinement level to 4. This controls the number of uniform refinements used by the multigrid solver and also the density of the resulting output grid. First click on the Export button to export the FeatFlow grid and data files corresponding to the defined problem. Then press the Solve button to call the external FeatFlow-cc2d solver. Finally, when the solver is done press the Import button to import the computed solution.
    ex_bfs_08_50.png


  19. To clearer see the recirculation zone length use the Surface Plot expression __(u<-eps)*x/hstep*(y<hstep)*(y>0)__. The Arrow Plot option can also be used to more clearly visualize the flow field.

    ex_bfs_09_50.png


  20. The resulting solution shows a recirculation zone length of about 6 length units which can be improved with a finer and better grid.

    ex_bfs_10_50.png


Flow over a Backwards Facing Step using the CLI

The process to set up and solve the flow over a backwards facing step test problem on the command line interface is illustrated in the ex_navierstokes4 script file which can be found in the examples directory.


References

[7] P.M. Gresho and R.L. Sani, Incompressible Flow and the Finite Element Method, Volume 1 & 2, John Wiley & Sons, New York, 2000.

[8] A. Rose and B. Simpson: Laminar, Constant-Temperature Flow Over a Backward Facing Step, 1st NAFEMS Workbook of CFD Examples, Glasgow, UK, 2000.


Structural Mechanics - Thin Plate with Hole

An example of setting up and solving the Thin Plate with Hole example with the FEATool graphical user interface (GUI) is described in the FEATool quickstart guide Structural Mechanics Example - Thin Plate with Hole.


Structural Mechanics - 3D Deflection of a Bracket

In this section a three dimensional model of a bracket with a hole will be described. The vertical side of the bracket is fixed to a wall (thus zero displacement boundary conditions are suitable here), while the outer end of the bracket is loaded with a negative force of -10000 in the z-direction. The resulting displacements will be studied together with creating a displacement plot visualizing the deformation.


3D Deflection of a Bracket using the GUI

This section describes how to set up and solve the bracket stress problem with the FEATool graphical user interface (GUI).

  1. Start Octave/Matlab, and if you have not run the installation script (which automatically adds the FEATool directory paths at startup) then change your working directory to where your FEATool installation is, for example
    cd C:\featool
    
  2. In the command window type

    featool
    

    to start the graphical user interface (GUI).

  3. Click on the New Model button in the upper horizontal toolbar to clear all data and start defining a new model.

    ex_poi1_01_50.png


  4. In the opened New Model dialog box, click on the 3D radio button in the Select Space Dimensions frame, and select Linear Elasticity from the Select Physics drop-down list. Leave the space dimension and dependent variable names to their default values. Finish and close the dialog box by clicking on the OK button.

    ex_le2_02_50.png


  5. To create the first block for the bracket, first click on the Create cube/block button in the left hand side Tools toolbar frame.

    ex_le2_03_50.png


  6. In the open Create block dialog box, set the extents of the block to [0 0.02], [0 0.2], [0 0.2], for the x, y, and z min and max coordinates. Press OK to finish and generate the first block.

    ex_le2_04_50.png


  7. Click on the Create cube/block button again to create the second section.

    ex_le2_05_50.png


  8. In the open Create block dialog box, now set the extents of the block to [0 0.2], [0 0.2], [0.09 0.11], for the x, y, and z min and max coordinates. Press OK to finish and generate the second block.

    ex_le2_06_50.png


  9. To create the cylindrical section click on the Create cylinder button .

    ex_le2_07_50.png


  10. In the open Create block dialog box, now set the cylinder center to [0.1 0.1 0.08], the radius to 0.08, the length to 0.04, and alignment axis to 3 (pointing in the z-direction). Press OK to finish and generate the cylinder.

    ex_le2_08_50.png


  11. To create the final geometry a composite object must be generated. To do this select Combine Objects... from the Geometry menu.

    ex_le2_09_50.png


  12. Enter the formula B1+B2-C1 in the opened dialog box and finish by clicking OK.

    ex_le2_10_50.png


  13. Switch to Grid mode by clicking on the corresponding mode toolbar button.

    ex_le2_11_50.png


  14. In grid mode, first set the Grid Size to 0.01 in the left toolbar edit field (this is half of the thickness so that a minimum of two cells are generated in each direction). Then click on the Generate unstructured grid button to call the grid generation function which automatically generates a grid of tetrahedra. Note that this can take some time depending on your system configuration.

    ex_le2_12_50.png


  15. Press the mode button in the Mode toolbar to switch from grid mode to physics and equation/subdomain specification mode.

    ex_le2_13_50.png


  16. In the Equation Settings dialog box that automatically opens, set Poisson's ratio $\nu$ to 0.3 and the modulus of elasticity E to 200e9 in the corresponding edit fields. The other coefficients can be left to their default values. Press OK to finish and close the dialog box.

    ex_le2_14_50.png


  17. Switch to boundary condition specification mode by clicking on the mode button

    ex_le2_15_50.png


  18. In the Boundary Settings dialog box, first select all boundaries (this can be done by pressing Ctrl+a after clicking in the Boundaries list box) and set all conditions to Edge loads with a value of zero 0. Press Apply to save these changes.

    ex_le2_16_50.png


  19. Now select the left vertical boundary (in this case boundary number 1) and select Fixed displacement in all directions. Press Apply again to save the changes.

    ex_le2_18_50.png


  20. Finally select the far right boundary (in this case boundary number 62) and select set an edge load in the z-direction with value -1e4. Press OK to finish prescribing boundary conditions.

    ex_le2_20_50.png


  21. Now that the problem is specified, press the mode button to switch to solve mode. Then press the button with an equals sign in the Tools toolbar frame to call the solver with the default solver settings.

    ex_le2_22_50.png


  22. After the problem has been solved FEATool will automatically switch to postprocessing mode and display the computed von Mieses stress. To change the plot, open the postprocessing settings dialog box by clicking on the Postprocessing settings button in the Tools toolbar frame.

    ex_le2_23_50.png


  23. In the Postprocessing settings dialog box choose to plot the z-displacement as a surface plot.

    ex_le2_24_50.png


  24. Three dimensional models can be rotated by using the Rotate button.

    ex_le2_25_50.png


  25. To create a displacement plot the model must be processed on the command line. To export the model choose Export, FEA Problem Struct To Main Workspace from the File Menu.

    ex_le2_26_50.png


  26. The following code scales and evaluates the displacements in the grid points and creates a distorted grid for visualization of the deformation.
    scaling = 5e3;
    dp = zeros(size(fea.grid.p));
    for i=1:3
      dp(i,:) = scaling*evalexpr( fea.dvar{i}, ...
                  fea.grid.p, fea );
    end
    fea_disp = fea;
    fea_disp.grid.p = fea_disp.grid.p + dp;
    
    f = figure;
    left_cells = selcells(fea,['x>',num2str(1.5*0.02)]);
    plotgrid( fea, 'facecolor', [.95 .95 .95], ...
                   'edgecolor', [.8 .8 1], ...
                   'selcells', left_cells, ...
                   'parent', gca )
    hold on
    plotgrid( fea_disp, 'parent', gca )
    
    title( 'Displacement plot' )
    view( [30 20] )
    


    ex_le2_27_50.png



3D Deflection of a Bracket using the CLI

The process to set up and solve the bracket stress test problem on the command line interface is illustrated in the ex_linearelasticity2 script file which can be found in the examples directory.


Parametric Study of Deflection of a Bracket

Parametric studies can be very useful tools to clearly see how a model will behave under varying conditions. The following section presents a parametric study of the bracket model varying both the geometry, in this case plate thickness, and also the applied load force.

  1. The first step in creating a parametric model is to save the model as an m-script model file. This can be done by selecting Save as M-Script Model from the File menu.

    ex_le2_28_50.png


  2. Open the file in your favorite editor. It is now possible to see exactly which steps and functions were used in creating the model. The first change to make is to add vectors of parameters at the top of the file
     thickness = [ 0.02 0.015 0.01 0.008 ];
     loads     = [ 1e4 2e4 1e5 1e6 ];
    
  3. The second change is to replace the geometry and grid operations with the following code so that the thickness of the second block will be varied with the iteration counter i
     for i=1:length(thickness)
    
     % Geometry operations.
     fea.geom = struct;
     gobj = gobj_block( 0, 0.02, 0, 0.2, 0, 0.2, 'B1' );
     fea  = geom_add_gobj( fea, gobj );
     gobj = gobj_block( 0, 0.2, 0, 0.2, ...
                   0.1-thickness(i)/2, 0.1+thickness(i)/2, 'B1' );
     fea  = geom_add_gobj( fea, gobj );
     gobj = gobj_cylinder( [ 0.1, 0.1, 0.08 ], 0.08, 0.04,3,'C1');
     fea  = geom_add_gobj( fea, gobj );
     fea  = geom_apply_formula( fea, 'B1+B2-C1' );
    
     % Grid generation.
     fea.grid = gridgen( fea, 'hmax', thickness(i)/2 );
    
  4. As second loop can now be introduced where the varying load is extracted and prescribed by using a constant variable. The equation settings can be kept from the saved model.

     for j=1:length(loads)
    
     % Constants and expressions.
     fea.expr = { 'load', loads(j) };
    
     % Equation settings.
     fea.phys.el.dvar = { 'u', 'v', 'w' };
     fea.phys.el.sfun = { 'sflag1', 'sflag1', 'sflag1' };
     fea.phys.el.eqn.coef = { 'nu_el', ', '', { '0.3' };
     ...
    

    Two loops have now been set up around the model, the first creating a geometry with varying thickness, and the second will prescribe the varying load. Note that setting up the loops in this way allows reusing the grid in the inner loop and saving computational time.

  5. When the grid changes the boundary numbering may also change. The static boundary condition prescription must therefore be replaced with one that identifies the correct boundaries. This can be achieved in the following way
     % Boundary conditions.
     n_bdr = max(fea.grid.b(3,:));
     ind_fixed = findbdr( fea, 'x<0.005' );
     ind_load  = findbdr( fea, 'x>0.19'  );
     fea.phys.el.bdr.sel = ones(1,n_bdr);
    
     % Define Homogeneous Neumann BC type.
     n_bdr  = max(fea.grid.b(3,:));
     bctype = num2cell( zeros(3,n_bdr) );
    
     % Set Dirichlet BCs for the fixed boundary.
     [bctype{:,ind_fixed}] = deal( 1 );
     fea.phys.el.bdr.coef{1,end-2} = bctype;
    
     % Define all zero BC coefficients.
     bccoef = num2cell( zeros(3,n_bdr) );
    
     % Apply negative z-load to outer edge.
     [bccoef{3,ind_load}] = deal( '-load' );
    
     % Store BC coefficient array.
     fea.phys.el.bdr.coef{1,end} = bccoef;
    
  6. When the problem is fully specified we can call the solver, apply postprocessing to calculate the maximum stress, and close the loops
    % Solver call.
    fea = parsephys( fea );
    fea = parseprob( fea );
    fea.sol.u = solvestat( fea );
    
    % Postprocessing.
    s_vm = fea.phys.el.eqn.vars{1,2};
    vm_stress = evalexpr( s_vm, fea.grid.p, fea );
    max_vm_stress(i,j) = max(vm_stress);
    
    end,end
    
  7. Finally visualize the all computed maximum stresses at once as
    % Visualization.
    [t,l] = meshgrid(thickness,loads);
    surf( t, l, max_vm_stress, log(max_vm_stress) )
    xlabel( 'Thickness' )
    ylabel( 'Load' )
    zlabel( 'Maximum stress' )
    view( 45, 30 )
    


    ex_le2_29_50.png


The complete parameter study model can be found as the ex_linearelasticity3 script file in the examples directory.


Multiphysics - Natural Convection in a Square Cavity

This section contains an example how to easily define and couple multiphysics models in FEATool. The chosen example is a benchmark model by De Vahl [9], [10] which simulates natural convection in a unit square through the Boussinesq approximation. The model consists of a Navier-Stokes equations physics mode representing the fluid flow with solid wall or no-slip boundary conditions everywhere. In addition a heat transfer physics mode is added for the temperature. The top and bottom boundaries are perfectly insulated while the left boundary has a temperature of 1 and the right zero.

The physics modes are two way coupled through the vertical source term in the Navier-Stokes equations, Pr*Ra*T (here the Prandtl and Rayleigh numbers are, Pr=0.71 and Ra=1e3 respectively), and the velocities transporting the temperature coming directly from the fluid flow. The references contain benchmark reference and comparison results for a number of quantities such as maximum velocities and the Nusselt number.


Natural Convection in a Square Cavity using the GUI

In the following it is described how to set up and solve a multiphysics natural convection benchmark problem the FEATool graphical user interface (GUI).

  1. Start Octave or Matlab, and if you have not run the installation script (which automatically adds the FEATool directory paths at startup) then change your working directory to where your FEATool installation is, for example
    cd C:\featool
    
  2. In the command window type

    featool
    

    to start the graphical user interface (GUI).

  3. Either select New... from the File menu, or click on the New Model button in the upper horizontal toolbar, to clear all data and start defining a new model.

    ex_poi1_01_50.png


  4. In the opened New Model dialog box, click on the 2D radio button in the Select Space Dimensions frame, and select Navier-Stokes Equations from the Select Physics drop-down list. Leave the space dimension and dependent variable names to their default values. Finish and close the dialog box by clicking on the OK button.

    ex_ns1_02_50.png


  5. To create the square geometry, click on the Create square/rectangle button in the left hand side Tools toolbar frame. Then left click in the main plot axes window, hold the mouse button, and move the mouse pointer to show a rectangle with red outlines.

    ex_ns1_04_50.png


  6. Release the button to finalize and create a solid geometry object. The object properties must now be edited to set the correct size and position of the rectangle. To do this, click on the rectangle R1 to select it and highlight it in red. (Alternatively you select it by clicking on R1 in the selection list box under the left side toolbar buttons.) Then click on the Inspect/edit selected geometry object toolbar button

    ex_ns1_05_50.png


  7. In the opened Edit Geometry Object dialog box, edit the minimum coordinates to 0 and maximum 1 to specify a unit square. Finish editing the geometry object and close the dialog box by clicking OK.

    ex_natconv_02_50.png


  8. An artificial point must now be created to prescribe a fixed pressure since there are no natural outflow conditions in this model (this is a requirement of the Navier-Stokes equations with enclosed domains). In the following we will create small boundary segments to approximate a point. To do this add another small rectangle in the lower left corner (in this case with side length 0.0025). Alternatively, it is also possible to use a point constraint for pressure which can be set by using Add Point Constraints... in the Boundary menu.

    ex_natconv_03_50.png


  9. Join both rectangles by selecting both R1 and R2 and clicking on the Add geometry objects button .

    ex_natconv_04_50.png


  10. Switch to Grid mode by clicking on the corresponding Mode toolbar button. Then enter 0.03 in the Grid Size edit field and click on the Generate unstructured grid button to call the automatic grid generation function.

    ex_natconv_04b_50.png


  11. Press the mode button to switch from grid mode to physics and equation/subdomain specification mode. In the Equation Settings dialog box that automatically opens, set the fluid density $\rho$ and viscosity $\mu$ both to 1, and the body/volume force in the y-direction to 0.71*1e3*T (this corresponds to the Pa*Ra*T, where T is the dependent variable for temperature which will be added next).

    ex_natconv_05_50.png


  12. To access the multiphysics selection and add another physics mode press the plus + tab (only available in the Multiphysics and Professional editions of FEATool). Now select Heat Transfer from the Select Physics drop down list. Finish by pressing the Add Physics >>> button.

    ex_natconv_06_50.png


  13. In the Equation Settings ht tab, set the density $\rho$, specific heat $Cp$, and heat conductivity $k$ to 1, while the source term Q is set to 0. The convective velocities should be coupled from the Navier-Stokes equations physics mode, to do this enter u and v in the corresponding edit fields (as these are the default names of the dependent variables for the velocities). Press OK to finish with the equation specifications.

    ex_natconv_07_50.png


  14. Switch to boundary condition specification mode by clicking on the mode button. Switch to the ns tab, which corresponds to the boundary conditions prescribed to the Navier-Stokes equations physics mode. Then select the Wall/no-slip for all boundaries except one of the small boundaries in the lower left corner (here all except boundary number 4).

    ex_natconv_08_50.png


  15. Select the left over boundary (number 4) and set an Outflow/pressure with $p_o$ set to 0 here. This prescribes a zero fixed pressure at this particular boundary and makes the problem well posed numerically.

    ex_natconv_08b_50.png


  16. Click on the ht tab to change to specifying boundary conditions for the heat transfer physics mode. Select Thermal insulation/symmetry for the top and bottom boundaries (here boundary numbers 2, 4, and 5).

    ex_natconv_09_50.png


  17. Select a Temperature boundary condition for the left boundary (here boundaries 3 and 6) with a fixed temperature $T_o$ set to 1.

    ex_natconv_10_50.png


  18. Similarly, also select the Temperature condition for the right boundary (boundary 1) but with a zero temperature, $T_o$ set to 0. Finish by clicking the OK button.

    ex_natconv_11_50.png


  19. Now that the problem is specified, press the mode button to switch to solve mode. Then press the button with an equals sign in the Tools toolbar frame to call the solver with the default solver settings.

  20. FEATool will automatically switch to postprocessing mode after a solution has been found. To change the plot, open the postprocessing settings dialog box by clicking on the Postprocessing settings button . To see the velocity select Velocity field in the predefined Surface Plot and Arrow Plot listbox options.

    ex_natconv_12_50.png


  21. Press OK or Apply to show the resulting plot for the velocity field.

    ex_natconv_13_50.png


  22. Similarly plot the temperature by selecting Temperature, T in the drop down boxes for Surface and Contour plots.

    ex_natconv_14_50.png


  23. Press OK or Apply again to show the resulting plot for the temperature field.

    ex_natconv_15_50.png


  24. It is possible to enter user defined postprocessing expressions, here showing how to enter the temperature gradient in the x-direction by simply entering Tx in the edit fields (T for the dependent variable name for the temperature and x for the derivative in the x-direction).

    ex_natconv_16_50.png


    ex_natconv_17_50.png


  25. The professional edition of FEATool also features a boundary integration option which can be used to calculate the average Nusselt number for the vertical boundaries. Here it is showing the resulting value of 1.0651 which can be compared with the reference value of 1.118).

    ex_natconv_18_50.png


  26. The lower command prompt can be used to input Octave and Matlab commands for advanced processing (or alternatively the whole FEATool fea struct can be exported to the main workspace for further processing). The following code shows how to calculate the maximum velocity and position along the vertical center line

    y = [0:0.01:1];
    x = 0.5*ones(size(y));
    
    u_eval = evalexpr( 'u', [x;y], fea );
    [u_max,ix] = max( u_eval )
    y_max = y( ix )
    

    The resulting values of u_max and y_max are here 2.6819 and 0.81, respectively, which is a little low compared to the reference solution (3.649 at a height of 0.813) which can indicate the the grid resolution is somewhat too coarse.

    ex_natconv_19_50.png



Natural Convection in a Square Cavity using the CLI

The process to set up and solve the natural convection problem on the command line interface is illustrated in the ex_natural_convection script file which can be found in the examples directory.


References

[9] D. de Vahl Davis, Natural Convection of Air in a Square Cavity: A Bench Mark Solution, Int. J. Numer. Meth. Fluids, vol. 3, pp. 249-264, 1983.

[10] D. de Vahl Davis and I. P. Jones, Natural Convection of Air in a Square Cavity: A Comparison Exercise, Int. J. Numer. Meth. Fluids, vol. 3, pp. 227-248, 1983.


Multiphysics - Heat Exchanger

An example of setting up and modeling a heat exchanger with coupled fluid flow and temperature with both free and forced convection is described in the FEATool quickstart guide Multiphysics Example - Heat Exchanger.


Multiphysics - Piezoelectric Bending of a Beam

This example describes how to model piezoelectric bending of a beam on the command line as an m-script model in FEATool. The test case comes from a research paper by Hwang and Park [11]. For this problem, the following coupled system of equations should be solved

\[ -\nabla\cdot \left[\begin{array}{c} \sigma \\ D \end{array}\right] = \left[\begin{array}{c} f \\ -\rho \end{array}\right] \]

where $\sigma$ is a stress tensor, $f$ represents body forces, $D$ electric displacement, and $\rho$ distributed charges.

By using the constitutive relations for a piezoelectric material [11], the stresses can in two dimensions be converted to the following form

\[ \nabla\cdot\left[ C \left[\begin{array}{c} \nabla u \\ \nabla v\ \\ \nabla V \end{array}\right]\right] = \left[\begin{array}{c} 0 \\ 0 \\ -\rho \end{array}\right] \]

where $C$ is an array of constitutive relations and material parameters. This leads to a system of three equations for the unknown displacements $u$ and $v$, and the potential $V$.

There are no readily available physics modes in FEATool yet so to model this we are defining the problem from scratch. In this case we will be looking at composition of two 1.2 cm long and 2 mm thick beams with opposite polarization.

First we assign names to the space coordinates (x and y), and create a grid. Note, that the shape is so simple here that we don't need to explicitly define a geometry, we can use the 'rectgrid' command directly

fea.sdim = { 'x' 'y' };
fea.grid = rectgrid( 20, 20, [ 0 12e-3; -1e-3 1e-3] );

Next is to define the equation coefficients and constitutive relations. See the attached file for their definitions, one thing to note that to define the polarization the coefficient is multiplied with a switch expression 2*(y<0)-1 which evaluates to 1 in the top half and -1 in the bottom (here we have also used the name y for the space coordinate in the 2nd direction as we defined earlier).

Now we define the dependent variables and assign 2nd order Lagrange finite element shape functions for them all

fea.dvar  = { 'u' 'v' 'V' };
fea.sfun  = { 'sflag2' 'sflag2' 'sflag2' };

The material parameters and constitutive relations are defined as

% Equation coefficients.
Emod =  2e9;       % Modulus of elasticity
nu   =  0.29;      % Poissons ratio
Gmod =  0.775e9;   % Shear modulus
d31  =  0.22e-10;  % Piezoelectric sn coefficient
d33  = -0.3e-10;   % Piezoelectric sn coefficient
prel = 12;         % Relative electrical permittivity
pvac = 0.885418781762-11; % Elec. perm. of vacuum

% Constitutive relations.
constrel   = [ Emod/(1-nu^2)    nu*Emod/(1-nu^2) 0    ;
               nu*Emod/(1-nu^2) Emod/(1-nu^2)    0    ;
               0                0                Gmod ];
piezoel_st = [ 0 d31 ;
               0 d33 ;
               0 0   ];
piezoel_sn = constrel*piezoel_st;
dielmat_st = [ prel 0    ;
               0    prel ]*pvac ;
dielmat_sn = [ dielmat_st - piezoel_st'*piezoel_sn ];

% Populate coefficient matrices (negative sign
% due to fem partial integration).
c{1,1} = { constrel(1,1)   constrel(1,3) ;
           constrel(1,3)   constrel(3,3) };
c{1,2} = { constrel(1,3)   constrel(1,2) ;
           constrel(3,3)   constrel(2,3) };
c{2,1} = c{1,2}';
c{2,2} = { constrel(3,3)   constrel(1,3) ;
           constrel(1,3)   constrel(2,2) };
c{1,3} = { piezoel_sn(1,1) piezoel_sn(1,2) ;
           piezoel_sn(3,1) piezoel_sn(3,2) };
for i=1:4
  c{1,3}{i} = [num2str(c{1,3}{i}),'*(2*(y<0)-1)'];
end
c{3,1} = c{1,3}';
c{2,3} = { piezoel_sn(3,1) piezoel_sn(3,2);
           piezoel_sn(2,1) piezoel_sn(2,2) };
for i=1:4
  c{2,3}{i} = [num2str(c{2,3}{i}),'*(2*(y<0)-1)'];
end
c{3,2} = c{2,3}';
c{3,3} = { dielmat_sn(1,1) dielmat_sn(2,1) ;
           dielmat_sn(2,1) dielmat_sn(2,2) };

To define the finite element bilinear forms we use the following code

bilinear_form = [ 2 2 3 3 ;
                  2 3 2 3 ];
for i=1:length(fea.dvar)
  for j=1:length(fea.dvar)
    fea.eqn.a.form{i,j} = bilinear_form;
    fea.eqn.a.coef{i,j} = c{i,j}(:)';
  end
end

where each bilinear form has four terms, the top line in the definition defines the finite element shape for the dependent variable, and second line for the test function space (a 2 indicates x-derivative, and 3 y-derivative).

In this case all three source terms are zero, thus

fea.eqn.f.form = { 1 1 1 };
fea.eqn.f.coef = { 0 0 0 };

The boundary conditions are defined as follows

n_bdr                   = max(fea.grid.b(3,:));
fea.bdr.d               = cell(3,n_bdr);
fea.bdr.n               = cell(3,n_bdr);
i_top                   = 3;
i_bottom                = 1;
i_left                  = 4;
[fea.bdr.d{3,i_top}]    = deal( 200 );
[fea.bdr.d{3,i_bottom}] = deal( 0 );
[fea.bdr.d{1:2,i_left}] = deal( 0 );

In this way we have set a potential difference of 200, and fixed the left hand side boundary. Now we parse the problem struct and solve the system with

fea       = parseprob( fea );
fea.sol.u = solvestat( fea );

After the problem has been solved we can for example visualize the displacement in the y-direction with the following commands

plotsubd( fea, 'labels', 'off', 'setaxes', 'off' )
u = evalexpr( 'u', fea.grid.p, fea );
v = evalexpr( 'v', fea.grid.p, fea );

DSCALE = 5000;
fea_disp = fea;
fea_disp.grid.p = fea.grid.p + [ u v ]'*DSCALE;
postplot( fea_disp, 'surfexpr', 'v', 'isoexpr', 'v' )



ex_piezo1_50.png


An example m-script file for this model can be found in the ex_piezoelectric1 m-script file.


References

[11] Hwang, W. S., Park, H. C. Finite Element Modeling of Piezoelectric Sensors and Actuators. AIAA Journal, 31(5), pp 930-937, 1993.


Custom Equation - Black-Scholes model equation

This section contains an example how to easily define a custom equation or model in FEATool. In this case a Black-Scholes model equation is studied which is often used in financial analytics. The equation to be modeled reads

\[ \frac{\partial u}{\partial t} - \frac{\partial}{\partial x}(\frac{1}{2}\frac{\partial u}{\partial x}) - \frac{\partial u}{\partial x} = -u + ( (x-t)^5 - 10(x-t)^4 - 10(x-t)^3) \]

with boundary conditions $u(0,t) = -t^5$ and $u(x,t) = (x-t)^5$ on the left and right sides of the domain, respectively. Initial condition is taken as $u(x,t=0) = x^5$ For this problem an exact analytical solution exists, namely $u(x,t) = (x-t)^5$. Consult the Custom Equation section for details on the equation syntax.


Black-Scholes model equation using the GUI

This section describes how to set up and solve the custom equation example with the FEATool graphical user interface (GUI).

  1. Start Octave or Matlab, and if you have not run the installation script (which automatically adds the FEATool directory paths at startup) then change your working directory to where your FEATool installation is, for example
    cd C:\featool
    
  2. In the command window type

    featool
    

    to start the graphical user interface (GUI).

  3. Click on the New Model button in the upper horizontal toolbar to clear all data and start defining a new model.

    ex_poi1_01_50.png


  4. In the opened New Model dialog box, click on the 1D radio button in the Select Space Dimensions frame, and select Custom Equation from the Select Physics drop-down list. Leave the space dimension and dependent variable names to their default values. Finish and close the dialog box by clicking on the OK button.

    ex_ce1_01_50.png


  5. To create a line grid, first click on the Create new line grid button in the left hand side Tools toolbar frame.

    ex_ce1_02_50.png


  6. Enter 0:0.05:1 in the Line Coordinates edit field of the Create Line Grid dialog box. This will specify a grid between 0 and 1 with 20 (1/0.05) cells. Press OK to finish and close the dialog box.

    ex_ce1_03_50.png


  7. Press the mode button in the Mode toolbar to switch from grid mode to physics and equation/subdomain specification mode.

    ex_ce1_04_50.png


  8. Press the edit eqn button in the Equation Settings dialog box that automatically opens.

    ex_ce1_05_50.png

  9. Enter u' - 1/2*ux_x - ux_t = -u + ( (x-t)^5 - 10*(x-t)^4 - 10*(x-t)^3 ) in the Edit Equations dialog box. Press OK to finish.

    ex_ce1_06_50.png


  10. Enter x^5 in the Initial Conditions edit field. Press OK to finish with the equation specifications and close the dialog box.

    ex_ce1_07_50.png


  11. Switch to boundary condition specification mode by clicking on the mode button

    ex_ce1_08_50.png


  12. In the Boundary Settings dialog box, first select the left boundary, boundary 1 in the Boundaries: list box. Choose the Dirichlet, r_u condition, and enter -t^5 in the corresponding edit field.

    ex_ce1_09_50.png


  1. Similarly for the right boundary, select boundary 2 in the Boundaries: list box and choose the Dirichlet, r_u condition, and enter (1-t)^5 in the corresponding edit field. Press OK to finish the boundary condition specifications.

    ex_ce1_10_50.png


  2. Now that the problem is specified, press the mode button to switch to solve mode. Then press the Solver settings toolbar button.

    ex_ce1_11_50.png


  3. Select the Time-Dependent - Monolithic (Fully coupled) solver option, then press Apply and Solve to start the solution process.

    ex_ce1_12_50.png


  4. After the problem has been solved FEATool will automatically switch to postprocessing mode and display the computed solution at the last time step. To change the plot open the postprocessing settings dialog box by clicking on the Postprocessing settings button in the Tools toolbar frame.

    ex_ce1_13_50.png


  5. In the Postprocessing settings dialog box enter u-(x-t)^5 to show the difference between the computed and exact analytical solution.

    ex_ce1_14_50.png


  6. To zoom in on the error open the Axis/Grid Settings... in the Options menu. Switch off Axis equal and manually set the Axis limits to 0 1 -2e-2 2e-2. We can now see that the maximum error has a magnitude of 1e-2 which can be further improved by using a finer grid and smaller time steps.

    ex_ce1_15_50.png



Black-Scholes model equation using the CLI

The process to set up and solve the custom equation problem on the command line interface is illustrated in the ex_custom_equation1 script file which can be found in the examples directory.


Custom Equation - Wave Equation

In the quickstart guide there is described how to set up a model of the wave equation with the custom equations physics mode (Custom Equation Example - Wave Equation on a Circle). In this case the custom equation uses two equations/dependent variables do describe the system.


Classic Equation - Poisson Equation on a Circle

The classic Poisson equation is one of the most fundamental partial differential equations (PDEs). Although one of the simplest equations, it is a very good model for the process of diffusion and comes up again and again in many applications (for example fluid flow, heat transfer, and chemical transport). It is thus very fundamental to many simulation codes to be able to solve it correctly and efficiently.

This example shows how to up and solve the Poisson equation

\[ d_{ts}\frac{\partial u}{\partial t} + \nabla\cdot(-D\nabla u) = f \eqno(1) \]

for a scalar field $u=u(\bf{x})$ on a circle $\Omega$ with radius $r=1$ in two dimensions. Both the diffusion coefficient $D$ and right hand side source term $f$ are assumed constant and equal to 1. The Poisson problem is also considered stationary meaning the time dependent term can be neglected. With these assumptions equation (1) simplifies to

\[ - \Delta u = 1. \]

Moreover, homogeneous Dirichlet boundary conditions are prescribed on all boundaries of the domain, that is $u=0\ on\ \partial\Omega$. The exact solution for this problem is $u(x,y)=(1-x^2-y^2)/4$ which can be used to measure the accuracy of the computed solution.


Poisson Equation using the GUI

This section describes how to set up and solve the Poisson equation (1) with the FEATool graphical user interface (GUI) which is available when using FEATool together with Octave version 4.0 or later and Matlab.

  1. Start Octave/Matlab, and if you have not run the installation script (which automatically adds the FEATool directory paths at startup) then change your working directory to where your FEATool installation is, for example
    cd C:\featool
    
  2. In the command window type

    featool
    

    to start the graphical user interface (GUI).

  3. Either select New... from the File menu, or click on the New Model button in the upper horizontal toolbar, to clear all data and start defining a new model.

    ex_poi1_01_50.png


  4. In the opened New Model dialog box, click on the 2D radio button in the Select Space Dimensions frame, and select Poisson Equation from the Select Physics drop-down list. Leave the space dimension and dependent variable names to their default values. Finish and close the dialog box by clicking on the OK button.

    ex_poi1_02_50.png


  5. To create a circle, first left click on the Create circle/ellipse button in the left hand side Tools toolbar frame.

    ex_poi1_03_50.png


  6. Then click and hold the left mouse button anywhere in the main plot axes window, and move the mouse pointer to show red outlines of a circle or ellipse. Release the button to finalize and create a solid geometry object.

    ex_poi1_04_50.png


  7. The object properties must be changed to make a circle with radius 1 centered at the origin. To do this, click on the ellipse E1 to select it which also highlights it in red (alternatively you select it by clicking on R1 in the selection list box under the left side toolbar buttons). Then click on the Inspect/edit selected geometry object toolbar button

    ex_poi1_05_50.png


  8. In the opened Geometry Object dialog box change the center coordinates edit field 0 0, and the x radius and y radius to 1 in the corresponding fields. Finish editing the geometry object and close the dialog box by clicking OK.

    ex_poi1_05b_50.png


  9. Press the mode button in the Mode toolbar to switch from geometry mode to grid generation mode.

    ex_poi1_06_50.png


  10. Click on the Generate unstructured grid button to call the grid generation function which automatically generates a grid of triangles for the circle.

    ex_poi1_06b_50.png


  11. Press the mode button in the Mode toolbar to switch from grid mode to physics and equation/subdomain specification mode.

    ex_poi1_08_50.png


  12. In the Equation Settings dialog box that automatically opens, set both the diffusion coefficient $D$ and source term coefficient $f$ to 1 in the corresponding edit fields. All other coefficients can be left to their default values. Press OK to finish and close the dialog box.

    ex_poi1_09_50.png


  13. Switch to boundary condition specification mode by clicking on the mode button

    ex_poi1_10_50.png


  14. In the Boundary Settings dialog box, select all boundaries in the left hand side Boundaries list box and choose Dirichlet boundary condition in the drop-down list. Set the Dirichlet boundary coefficient r equal to 0 in the Boundary Coefficients frame and finish by clicking on OK.

    ex_poi1_11_50.png


  15. Now that the problem is specified, press the mode button to switch to solve mode. Then press the button with an equals sign in the Tools toolbar frame to call the solver with the default solver settings.

    ex_poi1_13_50.png


  16. After the problem has been solved FEATool will automatically switch to postprocessing mode and display the computed solution. To change the plot, open the postprocessing settings dialog box by clicking on the Postprocessing settings button in the Tools toolbar frame.

    ex_poi1_14_50.png


  17. Activate Contour and Arrow plots by marking the corresponding check boxes and press OK or Apply to show the new plot.

    ex_poi1_15_50.png


  18. The solution is now visualized as surface, contour, and arrow plots in the main plot axes window.

    ex_poi1_16_50.png


  19. By going back to the Postprocessing settings dialog box and entering the expression u-(1-x^2-y^2)/4 in the Surface Plot expression edit field it is possible to plot and visualize the difference between the computed and exact reference solution.

    ex_poi1_17_50.png


  20. The error in the computed solution is now plotted, which clearly shows that the errors are localized to the top and bottom boundaries. From the colorbar we can also see that the maximum error has a magnitude of 1.5e-3.

    ex_poi1_18_50.png


  21. To further reduce the error one can recompute the problem with a finer grid, and also increase the solution space approximation order by choosing a higher finite element shape function in the subdomain settings.


Poisson Equation using the CLI and Physics Modes

This section describes how to set up and solve the Poisson equation (1) with the Octave or Matlab command line interface (CLI). To make it easier and faster for users to set up model problems several physics modes have been predefined with equations and boundary conditions for physics such as convection, diffusion, and reaction of chemical species, heat transport through convection and conduction, and incompressible fluid flow with the Navier-Stokes equations. In this example the Poisson physics mode will be used.

  1. Start Octave/Matlab and make sure you have either run the installation script to set up the FEATool directory paths or alternatively run

    addpaths
    

    on the command line from the FEATool installation directory before you start.

  2. Define a circle geometry object
    r    = 1;     % Circle radius.
    xc   = 0;     % x-center coordinate.
    yc   = 0;     % y-center coordinate.
    
    gobj = gobj_circle( [xc yc], r );
    fea.geom.objects  = { gobj };
    
  3. Call gridgen to create an unstructured grid with maximum grid cell size hmax
    hmax     = 0.1;   % Maximum grid cell size.
    fea.grid = gridgen( fea, 'hmax', hmax );
    
  4. Assign names to the space dimension coordinates (these names can then be used in equation coefficients and postprocessing expressions)
    fea.sdim = { 'x' 'y' };
    
  5. Use addphys to add the predefined Poisson physics mode poisson
    fea = addphys( fea, @poisson );
    
  6. To find and set the diffusion coefficient and source term to 1 first Inspect the 'eqn.coef' field and verify that the d_poi coefficient and f_poi source term are assigned rows 2 and 3, respectively. The coefficients values are specified in the fourth column of eqn.coef as a cell array of values per subdomain (in this case only one subdomain)
    fea.phys.poi.eqn
    
  7. Now set the correct coefficients
    fea.phys.poi.eqn.coef{2,4} = {1};   % Diffusion coefficient.
    fea.phys.poi.eqn.coef{3,4} = {1};   % Source term.
    
  8. To find and set the Dirichlet boundary coefficients to 0 first inspect the bdr.coef field and verify that the Dirichlet conditions are used with the coefficient set to zero
    fea.phys.poi.bdr.coef
    
  9. Then set the correct boundary coefficients
    % Use first boundary condition specification.
    fea.phys.poi.bdr.sel       = 1;
    
    % Set Dirichlet coefficient to zero.
    fea.phys.poi.bdr.coef{1,7} = {0};
    
    % Set flag to use Dirichlet BC for first selection.
    fea.phys.poi.bdr.coef{1,5} = {1};
    
    % Set flag not to use Neumann BC for first selection.
    fea.phys.poi.bdr.coef{2,5} = {0};
    
  10. Check and parse the physics and problem structs by calling parsephys and parseprob (this verifies that everything in the problem definition struct is correct and outputs the necessary lower hierarchy fields)
    fea = parsephys( fea );
    fea = parseprob( fea );
    
  11. Call the stationary solver solvestat with the default settings
    fea.sol.u = solvestat(fea);
    
  12. Plot with postplot to display the scalar solution as surface and contour plots, and the solution gradient (ux, uy) as arrow plots
    postplot( fea, 'surfexpr', 'u', ...
                   'isoexpr',  'u', ...
                   'arrowexpr', {'ux','uy'} )
    
  13. Plot the difference between the computed and exact solution
    u_diff = 'u-(1-x^2-y^2)/4';
    clf
    postplot( fea, 'surfexpr', u_diff )
    


Poisson Equation using the CLI without Physics Modes

This section describes how to set up and solve the Poisson equation (1) with the Octave or Matlab command line interface (CLI) by directly setting up the necessary equation and boundary definition fields (that is without using the predefined Poisson physics mode).

  1. Start Octave/Matlab and make sure you have either run the installation script to set up the FEATool directory paths or alternatively run

    addpaths
    

    on the command line from the FEATool installation directory before you start.

  2. Define a circle geometry object
    r    = 1;     % Circle radius.
    xc   = 0;     % x-center coordinate.
    yc   = 0;     % y-center coordinate.
    
    gobj = gobj_circle( [xc yc], r );
    fea.geom.objects  = { gobj };
    
  3. Call gridgen to create an unstructured grid with maximum grid cell size hmax
    hmax     = 0.1;   % Maximum grid cell size.
    fea.grid = gridgen( fea, 'hmax', hmax );
    
  4. Assign names to the space dimension coordinates (these names can then be used in equation coefficients and postprocessing expressions)
    fea.sdim = { 'x' 'y' };
    
  5. Assign a name to the dependent variable (in this case u)
    fea.dvar = { 'u' };
    
  6. Define a finite element shape function to use for the discretization (here sflag1 corresponds to linear conforming P1 finite element shape functions, while sflag2 would correspond to second order P2 shape functions.)
    fea.sfun = { 'sflag1' };
    
  7. Define the equation system, bilinear form in fea.eqn.a and right hand side linear form in fea.eqn.f
    % Define bilinear form. The first row indicates test function
    % space, and second row trial function space. Each column
    % defines a term in the bilinear form and 1 corresponds to
    % a basis function evaluation, 2 = x-derivative,
    % 3 = y-derivative, (d 4 = z-derivative in 3D).
    fea.eqn.a.form = { [2 3; 2 3] };
    
    % Define coefficients used in assembling the bilinear forms.
    fea.eqn.a.coef = { [1 1] };
    
    % Test function space to evaluate in linear form.
    fea.eqn.f.form = { 1 };
    
    % Coefficient value used in assembling the linear form.
    fea.eqn.f.coef = { 1 };
    
  8. Define the boundary conditions, Dirichlet conditions in fea.bdr.d and Neumann (flux conditions) in fea.bdr.n
    % Number of boundary segments.
    n_bdr          = max(fea.grid.b(3,:));
    
    % Allocate space for n_bdr boundaries.
    fea.bdr.d      = cell(1,n_bdr);
    
    % Assign u=0 to all Dirichlet boundaries.
    [fea.bdr.d{:}] = deal(0);
    
    % No Neumann boundaries (fea.bdr.n{i} empty).
    fea.bdr.n      = cell(1,n_bdr);
    
  9. Check and parse the problem struct by calling parseprob (this verifies that everything in the problem definition struct is correct and outputs the necessary lower hierarchy fields)
    fea = parseprob( fea );
    
  10. Call the stationary solver solvestat with the default settings
    fea.sol.u = solvestat( fea );
    
  11. Plot the results with postplot to display the scalar solution as and contour plots, and the solution gradient (ux, uy) as arrow plots
    postplot( fea, 'surfexpr', 'u', ...
                   'isoexpr',  'u', ...
                   'arrowexpr', {'ux','uy'} )
    
  12. Plot the difference between the computed and exact solution
    u_diff = 'u-(1-x^2-y^2)/4';
    clf
    postplot( fea, 'surfexpr', u_diff )
    


Poisson Equation using the CLI and core FEM library functions

This final section shows how to define and solve the Poisson equation (1) by directly calling the core finite element library functions. The low level functions are used to assemble the system matrix and right hand side vector to which boundary conditions are applied, after which the system can be solved.

  1. Start Octave/Matlab and make sure you have either run the installation script to set up the FEATool directory paths or alternatively run

    addpaths
    

    on the command line from the FEATool installation directory before you start.

  2. Define a circle geometry object
    r    = 1;     % Circle radius.
    xc   = 0;     % x-center coordinate.
    yc   = 0;     % y-center coordinate.
    
    gobj = gobj_circle( [xc yc], r );
    fea.geom.objects  = { gobj };
    
  3. Call gridgen to create an unstructured grid with maximum grid cell size hmax
    hmax     = 0.1;   % Maximum grid cell size.
    fea.grid = gridgen( fea, 'hmax', hmax );
    
  4. Use the function assemblea to assemble the finite element stiffness matrix
    % Bilinear form to assemble (u_x*v_x+u_y*v_y).
    form  = [2 3;2 3];
    
    % FEM shape functions used in test and trial function spaces
    % (here first order conforming P1 Lagrange functions).
    sfun  = {'sflag1';'sflag1'};
    
    % Coefficients to use for each term in the bilinear form.
    coef  = [1 1];
    i_cub = 3;   % Numerical quadrature rule to use.
    
    [vRowInds,vColInds,vAvals,n_rows,n_cols] = ...
         assemblea( form, sfun, coef, i_cub,  ...
           fea.grid.p, fea.grid.c, fea.grid.a );
    
    % Construct sparse matrix.
    A = sparse( vRowInds, vColInds, vAvals, n_rows, n_cols );
    
  5. Assemble right hand side load vector by calling assemblef
    % Linear form to assemble (1 = evaluate function values).
    form  = [1];
    
    % Finite element shape function.
    sfun  = {'sflag1'};
    
    % Coefficient to use in the linear form.
    coef  = [1];
    
    % Numerical quadrature rule to use.
    i_cub = 3;
    
    f = assemblef( form, sfun, coef, i_cub, ...
          fea.grid.p, fea.grid.c, fea.grid.a );
    
  6. Find indexes to boundary nodes
    bind = [];
    for i_b=1:size(fea.grid.b,2)    % Loop over boundary edges.
      i_c    = fea.grid.b(1,i_b);   % Cell number.
      i_edge = fea.grid.b(2,i_b);
      j_edge = mod( i_edge, size(fea.grid.c,1) ) + 1;
      bind   = [bind fea.grid.c([i_edge j_edge],i_c)'];
    end
    bind = unique(bind);
    
  7. Set homogeneous Dirichlet boundary conditions
    A = A';   % Transpose matrix for faster row elimination.
    A(:,bind) = 0;          % Zero/remove entire rows.
    for i=1:length(bind)    % Set diagonal entries of Dirichlet
                            % boundary condition rows to 1.
      i_a        = bind(i);
      A(i_a,i_a) = 1;
    end
    A = A';   %'
    
    f(bind) = 0;            % Set corresponding RHS entries
                            % to the prescribed boundary values.
    
  1. Solve (invert) problem by calling the built-in linear solver
    u = A\f;
    
  2. Create mock fea fields for postprocessing
    fea.sdim  = { 'x' 'y' };
    fea.dvar  = { 'u' };
    fea.sfun  = { 'sflag1' };
    fea       = parseprob( fea );
    
  3. Use postplot to plot the solution
    fea.sol.u = u;
    postplot( fea, 'surfexpr', 'u' )
    
  4. Plot the difference between the computed and exact solution
    u_diff = u - ( 1 - fea.grid.p(1,:)'.^2 - ...
                       fea.grid.p(2,:)'.^2 )/4;
    clf
    fea.sol.u = u_diff;
    postplot( fea, 'surfexpr', 'u' )
    


Classic Equation - Poisson Equation with a Point Source

The Classic Equation Example - Poisson Equation with a Point Source example in the quickstart guide shows how to model a circle with a point source at its center.


Classic Equation - Poisson Equation for a Sphere

The following example illustrates how to set up and solve the Poisson equation on a unit sphere $-\Delta u=1$ ( $r=0..1$) with homogeneous boundary conditions, that is $u=0|_{r=1}$ with the FEATool FEM assembly primitives.

To derive the finite element discretization of the general Poisson equation $ -\Delta u = f$ first multiply the equation with an arbitrary function $v$ (called test function), and integrate over the whole domain $\Omega$, that is $\int_{\Omega} -\Delta u v\ d\Omega=\int_{\Omega} fv\ d\Omega$. By applying the Gauss theorem or partial integration to the left side we get

\[ \int_{\Omega} \nabla u\cdot\nabla v\ d\Omega + \int_{\partial\Omega} -\hat{\bf{n}}\cdot\nabla u\ dS=\int_{\Omega} fv\ d\Omega \]

The boundary term can be neglected assuming that we only have prescribed or fixed value (Dirichlet) boundary conditions. In 3D the equation will look like

\[ \int_{\Omega} \frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y}+\frac{\partial u}{\partial z}\frac{\partial v}{\partial z}\ dxdydz = \int_{\Omega} fv\ dxdydz \]

By discretizing the solution (trial function space) variable as $u=\sum_{i=1}^N \phi_i u_i$ and similarly test function $v=\sum_{j=1}^N \phi_j$, where $\phi_{i/j}$ are the finite element shape or basis functions (usually taken the same for $i$ and $j$ in the Galerkin approximation), and $u_i$ is the value of the solution at node $i$. Inserting this we get

\[ \int\sum_{i=1}^N\sum_{j=1}^N \frac{\partial \phi_i}{\partial x}\frac{\partial \phi_j}{\partial x}+\frac{\partial \phi_i}{\partial y}\frac{\partial \phi_j}{\partial y}+\frac{\partial \phi_i}{\partial z}\frac{\partial \phi_j}{\partial z}\ u_i = \int\sum_{j=1}^N f\phi_j \]

which after discretization and integration with numerical quadrature gives us a matrix system to solve

\[ Au=b \]

To start implementing the model pragmatically we first have to create a grid. Instead of creating a grid by using a geometry object we use a geometric sphere grid primitive directly

grid = spheregrid();
p = grid.p;   % Extract grid points.
c = grid.c;   % Extract grid cells.
a = grid.a;   % Extract adjacency information.

To assemble the matrix A we can identify three bi-linear terms (with a coefficient of one) $1\cdot\phi_{i,x}\phi_{i,y}+1\cdot\phi_{i,y}\phi_{j,y}+1\cdot\phi_{i,z}\phi_{j,z}$. To specify this in FEATool first create a cell array with coefficient values (or expressions) for each term, in this case coef = {1 1 1};. Function values and/or derivatives for the trial and test functions are specified in a 2 by n-terms vector, here form = [2 3 4;2 3 4]; where the first row specifies trial functions to use for the three terms (1=function value, 2=x-derivative, 3=y-derivative, 4=z-derivative), and the bottom row the test functions. Lastly, the actual shape or basis functions are specified in a cell array as sfun = {'sflag1';'sflag1'};, which prescribes first order linear Lagrange shape functions for both the trial and test function spaces. Given a grid with grid points p and cell connectivity c the call to the fem assembly routine looks like

coef   = { 1 1 1 };
form   = [ 2 3 4 ; ...
           2 3 4 ];
sfun   = { 'sflag1' ; ...
           'sflag1' };
i_cub  = 3;              % Numerical quadrature rule.
cind   = 1:size(c,2);    % Index to cells to assemble.
[i,j,s,m,n] = assemblea( form, sfun, coef, i_cub, p, c, a );
A      = sparse( i, j, s, m, n );

Similarly, the right hand side linear form or load vector b can be assembled as

coef   = { 1 };
form   = [ 1 ];
sfun   = { 'sflag1' };
i_cub  = 3;              % Numerical quadrature rule.
cind   = 1:size(c,2);    % Index to cells to assemble.
b      = assemblef( form, sfun, coef, i_cub, p, c, a );

Prescribing homogeneous (zero value) Dirichlet boundary conditions on the boundary of the domain can be done by first finding the indexes to the nodes on the boundary (in this case the unit sphere), and then setting the corresponding values in the load vector to zero. Moreover, the corresponding rows in the matrix A must also be modified so that they are all zero except for the diagonal entries which are set to one (this preserves the boundary values specified in the b vector).

radi = sqrt(sum(p.^2,1));        % Radius of grid points.
bind = find(radi>1-sqrt(eps));   % Find grid points
                                 % on the boundary.

b(bind)   = 0;            % Set 'bind' rhs entries to zero.

A(bind,:) = 0;            % Set 'bind' rows to zero.
for i=1:length(bind)      % Set 'bind' diagonal entries to 1
  i_a        = bind(i);
  A(i_a,i_a) = 1;
end

The problem is now fully specified and can be solved as

u = A\b;

where u now contains the values of the solution in the grid points/nodes. To visualize the solution we put everything in a FEATool fea struct

fea.grid  = grid;
fea.dvar  = { 'u' };
fea.sfun  = sfun;
fea.sol.u = u;
fea = parseprob( fea );

postplot( fea, 'surfexpr', 'u', 'selcells', 'y>-0.1' )



poisson3d_50.png


Other Examples

The following selection of command line examples and test cases can be found in the examples directory.

FieldDescription
ex_brinkman1Coupled porous and incompressible laminar flow.
ex_convdiff1Convection and diffusion on a unit square.
ex_convdiff2Periodic convection and diffusion on a line with exact solution.
ex_convdiff3Infinite time dependent convection and diffusion on a line.
ex_convdiff4One dimensional Burgers equation with steady solution.
ex_convdiff5Dominating convection example requiring artificial stabilization.
ex_convdiff61D stationary convection-diffusion-reaction with exact solution.
ex_convreact1Time dependent one dimensional convection and reaction model.
ex_custom_equation1Black-Scholes model equation implemented as a custom equation.
ex_euler_beam11D Euler-Bernoulli beam example.
ex_euler_beam21D Euler-Bernoulli beam vibration mode example.
ex_nonlinear_pde11D nonlinear PDE with analytical solution.
ex_nonlinear_pde21D nonlinear PDE with analytical solution extended to 2D.
ex_diffusion1Diffusion equation on a unit square with different solutions.
ex_diffusion2Diffusion equation on a line with exact solutions.
ex_heattransfer12D heat conduction with natural convection and radiation.
ex_heattransfer2One dimensional stationary heat transfer with radiation.
ex_heattransfer3One dimensional transient heat conduction.
ex_heattransfer4Two dimensional heat transfer with convective cooling.
ex_heattransfer5Two dimensional transient cooling shrink fitting example.
ex_heattransfer6Axisymmetric steady state heat conduction of a cylinder.
ex_heattransfer7One dimensional transient heat conduction with analytic solution.
ex_laplace1Laplace equation on a unit square.
ex_laplace2Laplace equation on a unit circle.
ex_linearelasticity1Linear elasticity test case, solid stress-strain for a cube.
ex_linearelasticity2Three dimensional example of stress on a bracket.
ex_linearelasticity3Parametric study of the bracket deformation model.
ex_linearelasticity4Stress of loaded I-beam supported by two brackets.
ex_multiphase1Multiphase flow example of a static bubble.
ex_multiphase2Multiphase flow example of a rising bubble.
ex_multiphase3Multiphase flow example of a breaking dam.
ex_heat_exchanger1Free and forced convection in a heat exchanger.
ex_natural_convectionNatural convection in a square cavity.
ex_darcy1Porous media flow in a packed bed reactor using Darcy's law.
ex_navierstokes1Incompressible fluid flow in a channel.
ex_navierstokes2Incompressible driven cavity flow.
ex_navierstokes3Incompressible fluid flow around a cylinder in a channel (Re=20).
ex_navierstokes4Incompressible fluid flow over a backwards facing step.
ex_navierstokes5Two dimensional decay of a standing vortex.
ex_navierstokes6Time dependent flow around a cylinder in a channel (0<=Re<=100).
ex_navierstokes7Laminar flow in a curved three dimensional pipe.
ex_navierstokes8Axisymmetric flow in a circular pipe.
ex_navierstokes9Axisymmetric impinging jet.
ex_navierstokes103D flow in a pipe using the FeatFlow solver.
ex_navierstokes113D cavity flow in a cube using the FeatFlow solver.
ex_navierstokes123D flow over a backwards facing step using the FeatFlow solver.
ex_navierstokes133D flow past a cylinder using the FeatFlow solver.
ex_piezoelectric1Bending of a beam due to piezoelectric effects.
ex_electrostatics1Electrostatics test problem.
ex_electrostatics2Axisymmetric model of a spherical capacitor.
ex_resistive_heating1Resistive heating of a spiral tungsten filament.
ex_magnetostatics1Magnetostatics test problem.
ex_magnetostatics2Magnet field around a horseshoe magnet.
ex_magnetostatics3Cylindrical magnet without electrical currents.
ex_magnetostatics4Axisymmetric model of cylindrical magnet.
ex_magnetostatics52D Magnetostatics test model.
ex_planestress1Stress computation for a hole in a thin plate.
ex_planestress2NAFEMS plane stress benchmark example.
ex_planestress3NAFEMS plane stress benchmark of a tapered plate.
ex_planestress4NAFEMS benchmark of plate deformation due to thermal stress.
ex_planestress5NAFEMS benchmark for plane stress in an elliptic membrane.
ex_planestress62D Plane stress analysis of a pressure vessel.
ex_planestrain12D Plane strain analysis of a pressure vessel.
ex_axistressstrain1Axisymmetric stress-strain of a hollow cylinder.
ex_axistressstrain2Axisymmetric stress-strain of a hollow sphere.
ex_axistressstrain3Axisymmetric disc with fixed edge and point load.
ex_axistressstrain4Stress formed in brake disk due to heat build up.
ex_poisson1Poisson equation on a line.
ex_poisson2Poisson equation on a circle.
ex_poisson3Poisson equation on a unit square.
ex_poisson4Poisson equation on a rectangle with complex solution.
ex_poisson5Poisson equation on a sphere.
ex_poisson6Poisson equation on a unit cube.
ex_poisson7Poisson equation on a unit circle with a point constraint.
ex_waveequation1Wave equation on a unit circle.