FEATool
v1.8 Finite Element Analysis Toolbox |

Tutorials

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

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

where is the density, the heat capacity, is the thermal conductivity, heat source term, and a vector valued convective velocity field.

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 *T _{inf} = 17 C* and a heat transfer coefficient of

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

- 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

In the command window type

featool

to start the graphical user interface (GUI).

- Click on the
**New Model**button in the upper horizontal toolbar to clear all data and start defining a new model. - 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. - 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). - 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. - 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.
- 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. - 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 or using the slider control, above the grid toolbar buttons. Then click on the**Generate**button to call the automatic grid generation function. - 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 to**7500**, heat capacity to**470**, and thermal conductivity to**44**in the corresponding edit fields. Also set the initial temperature here to -10. - Now set the material parameters for the Tungsten part by selecting the other subdomains (1 and 3). Here set the density to
**19000**, heat capacity to**134**, thermal conductivity to**163**, and the initial temperature to**84**. Press**OK**to finish. - 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 (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 . - 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. - 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. 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]')

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.

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

A well known benchmark, test, and validation problem suite for incompressible fluid flows are the DFG cylinder benchmark problems. 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 in a by rectangular channel. The fluid density is taken as *1* and the viscosity is *0.001*. A fully developed parabolic velocity profile is prescribed at the inlet, , with the maximum velocity . This results in a mean velocity and a Reynolds number and thus the flow field will be laminar.

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

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 , and the coefficients of drag and lift , defined as

The drag and lift forces, and , can be computed as

where is the velocity in the tangential direction .

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

- 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

In the command window type

featool

to start the graphical user interface (GUI).

- Click on the
**New Model**button in the upper horizontal toolbar to clear all data and start defining a new model. - 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. - To create the outer rectangle, first 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.
- 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 - 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**. - To create the inner circle, first left click on the
**Create circle/ellipse**button in the left hand side*Tools*toolbar frame. - 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.
- 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 - 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**. - 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. - Switch to
*Grid*mode by click on the corresponding mode toolbar button. - 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**button to call the automatic grid generation function. - 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 to**1**and viscosity 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. - Switch to
*boundary condition specification*mode by clicking on the mode button - 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. - 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 . - 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 with the**Outflow/pressure**condition). Finish by clicking the**OK**button. - 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. - 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. - 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**button . - Now go back to
**Solve**mode and solve the problem again with the refined grid. - Now that the solution looks better we can calculate the drag and lift coefficients. To do this select
**Boundary Integration...**from the*Post*menu. - 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. - 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.

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.

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].

[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.

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.

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 . The size domain is chosen according to the data below. The inlet velocity function is given as where is the channel height, and maximum velocity 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.

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 external OpenFOAM CFD Solver.

- 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

In the command window type

featool

to start the graphical user interface (GUI).

- Click on the
**New Model**button in the upper horizontal toolbar to clear all data and start defining a new model. - 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. - To create the outer rectangle, first 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.
- 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 - 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**. - 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. - Switch to
*Grid*mode by click on the corresponding mode toolbar button. - Set a target grid size by using the
*Grid Size*slider control, or alternatively enter a specific grid size value manually in the corresponding edit field above. Then click on the**Generate**button to call the automatic grid generation function. - Although both FEATool and the
*OpenFOAM CFD*solvers support all grid cell types, it is sometimes desirable to work with structured quadrilateral and hexahedral cells. To convert to quadrilateral cells, select**Convert Grid Cells**from the**Grid**menu and apply grid smoothing by using the corresponding**Grid Smoothing**menu options. - 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 to**1**and viscosity 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. - Click on the
**Model Constrains and Expressions**and enter the constants and expressions shown below. - 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. When using the default built-in solver enter**4*u_max*(y-hstep)*(1-y)/h_inlet^2**, in the edit field for the x-velocity coefficient , to prescribe a parabolic laminar inlet profile, or when using OpenFOAM enter**3/2*u_max**to specify a constant mean velocity. - Select the right outflow boundary (here number 1) in the left hand side
*Boundaries*list box and choose the**Outflow/pressure**boundary condition from the drop-down list. Finish by clicking the**OK**button. - Now that the problem is specified, press the mode button to switch to
*solve*mode. To use the default built-in solver, press the button with an*equals*sign . Alternatively, press the button labeled**OpenFOAM**to open the OpenFOAM CFD solver settings and control panel. - In the
*OpenFOAM*solver settings dialog box, use the default settings, and press the**Solve**button to call the external*OpenFOAM*solver. When the solver is done the computed solution will automatically be imported, and FEATool will switch to Postprocessing mode. - To see the recirculation zone clearly, enter the expression for the normalized recirculation zone length __(u<-eps)*x/hstep*(y<hstep)*(y>0)__ in the
*Surface Plot*expression edit field. The**Arrow Plot**option can also be used to visualize the flow field. - The resulting solution shows a recirculation zone length of about 6 length units which can be improved with a finer and better grid.

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.

[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.

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.

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.

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

- 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

In the command window type

featool

to start the graphical user interface (GUI).

- Click on the
**New Model**button in the upper horizontal toolbar to clear all data and start defining a new model. - 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. - To create the first block for the bracket, first click on the
**Create cube/block**button in the left hand side*Tools*toolbar frame. - 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. - Click on the
**Create cube/block**button again to create the second section. - 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. - To create the cylindrical section click on the
**Create cylinder**button . - 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. - To create the final geometry a composite object must be generated. To do this select
**Combine Objects...**from the**Geometry**menu. - Enter the formula
**B1+B2-C1**in the opened dialog box and finish by clicking**OK**. - Switch to
*Grid*mode by clicking on the corresponding mode toolbar button. - 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**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. - 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 Poisson's ratio 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. - Switch to
*boundary condition specification*mode by clicking on the mode button - 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. - 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. - 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. - 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. - 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. - In the Postprocessing settings dialog box choose to plot the
**z-displacement**as a surface plot. - Three dimensional models can be rotated by using the
**Rotate**button. - 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**. - 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] )

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 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.

- 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. - 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 ];

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

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.

- 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;

- 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

- 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 )

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

Model example of an axisymmetric electrostatic spherical capacitor compared with analytical solutions.

- Start a new model and select
*Axi 2D*and*Electrostatics*. - Start in
*Geometry*mode and create three circles all centered at*0 0*but with radii*0.003*,*0.01*, and*0.012*. - Create three rectangles on the left side of the symmetry axis (
*r<0*) so that they cover the left side of the circles (extending a somewhat above, below and to the left). - One by one subtract the rectangles from the corresponding circles so that only half circles on right side half plane remain.
- Switch to
*Grid*mode by clicking on the corresponding mode toolbar button, and generate a grid with a the subdomain*Grid Size*equal to*0.001*. - 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, select all*Subdomains*(1-3) and enter the following coefficients in the*Equation Settings*dialog boxName Expression Permittivity, sigma+epsr*eps0/tscale Polarization, x-component, 0 Polarization, y-component, 0 Space charge density, rho/tscale - Select
*(P2/Q2) second order conforming*in the*FEM Discretizations*drop down box to icrease the spatial discretization accuracy. - Open the
*Model Constants and Expressions*dialog box and enter the following expressions (note that expressions for several subdomains are separated with one or more spaces)Name Expression r1 0.003 r2 0.01 r3 0.012 sigma 6e7 0 6e7 eps0 8.85e-12 epsr 1 3.9 1 tscale 1e-17 rho q0*3/4/pi/(r3^3-r2^3) 0 -q0*3/4/pi/r1^3 q0 6e-11 - Switch to
*boundary condition specification*mode by clicking on the mode button and select the following boundary conditions*Ground/antisymmetry V=0*for the outer boundaries (3 and 4).*Insulation/symmetry*for all symmetry boundaries*r=0*(1,2,5,6, and 7).

- Now that the problem is specified, press the mode button to switch to
*solve*mode. Then press the button with an*equals*sign to solve the problem. - After the problem has been solved, FEATool automatically changes to Postprocessing mode and plots the electric potential
*V*. - To verify the solution, choose the
*Subdomain Integration...*option in the*Post*menu, and integrate the expression*eeng*=*eps0*epsr*(Vr^2+Vz^2)*pi*r*over all subdomains. The capacitance can then be calculated as*q0^2/(2*eeng)*. - An alternative way to calculate the capacitance is
*q0/( Vmax - Vmin )*where*Vmax*and*Vmin*can be found by using the*Max Min Evaluation...*option from the*Post*menu. - Compare the computed capacitances with the analytical expression
*4*pi*epsr*eps0/( 1/r1 - 1/r2 )*which in this example should be equal to*1.8588e-12*.

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

Magnetostatics example for calculating the magnetic field around a horseshoe magnet.

- Start a new model and select
*2D*and*Magnetostatics*. - In
*Geometry*mode first create two circles, centered at*0 0*with radii*0.05*and*0.025*. - Create four rectangles with dimensions according to the table below
R1 R2 R3 R4 x_min -0.06 -0.05 0.025 -0.15 x_max 0.06 -0.025 0.05 0.15 y_min 0 0 0 -0.2 y_max 0.06 0.06 0.06 0.2 - Subtract rectangle
*R1*and the smaller circle*C2*from the larger one*C1*using the formula**C1-C2-R1**in the*Combine Objects...*option of the*Geometry*menu. - Switch to
*Grid*mode by clicking on the corresponding mode toolbar button, and generate a grid with a the subdomain*Grid Size*equal to*0.01*. - 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*Magnetization, M*coefficient to_{y}**1**and**-1**in the magnet ends, and**0**in the outer and magnet curved domains. - Select
*(P2/Q2) second order conforming*in the*FEM Discretizations*drop down box to icrease the spatial discretization accuracy. - In
*Boundary*mode, leave all the boundary conditions to the default*Magnetic insulation/antisymmetry*selection. - Now that the problem is specified, press the mode button to switch to
*solve*mode. Then press the button with an*equals*sign to solve the problem. - After the problem has been solved, FEATool automatically changes to
*Postprocessing*mode and plots the Magnetic potential*Vm*. - Contours and arrow plots of the magnetic field can be selected by using the the
*Postprocessing Settings*dialog box. - To verify the solution, choose the
*Subdomain Integration...*option in the*Post*menu, and integrate the*Magnetic potential, Az*and*Magnetic field*over all subdomains, the integrals should result in values around*-9.6e-11*and*0.0051*, respectively.

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

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.

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

- 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

In the command window type

featool

to start the graphical user interface (GUI).

- 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. - 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. - 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. - 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 - 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**. - 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. - Join both rectangles by selecting both
**R1**and**R2**and clicking on the**Add geometry objects**button . - 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**button to call the automatic grid generation function. - 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 and viscosity 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). - 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. - In the
*Equation Settings**ht*tab, set the density , specific heat , and heat conductivity 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. - 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). - Select the left over boundary (number 4) and set an
**Outflow/pressure**with set to**0**here. This prescribes a zero fixed pressure at this particular boundary and makes the problem well posed numerically. - 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). - Select a
**Temperature**boundary condition for the left boundary (here boundaries 3 and 6) with a fixed temperature set to**1**. - Similarly, also select the
**Temperature**condition for the right boundary (boundary 1) but with a zero temperature, set to**0**. Finish by clicking the**OK**button. - 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. - 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. - Press
**OK**or**Apply**to show the resulting plot for the velocity field. - Similarly plot the temperature by selecting
**Temperature, T**in the drop down boxes for*Surface*and*Contour*plots. - Press
**OK**or**Apply**again to show the resulting plot for the temperature field. - 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). - 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).
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 liney = [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.

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.

[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.

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.

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

where is a stress tensor, represents body forces, electric displacement, and distributed charges.

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

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

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' )

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

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

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

with boundary conditions and on the left and right sides of the domain, respectively. Initial condition is taken as For this problem an exact analytical solution exists, namely . Consult the Custom Equation section for details on the equation syntax.

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

- 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

In the command window type

featool

to start the graphical user interface (GUI).

- Click on the
**New Model**button in the upper horizontal toolbar to clear all data and start defining a new model. - 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. - To create a line grid, first click on the
**Create new line grid**button in the left hand side*Tools*toolbar frame. - 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. - Press the mode button in the
*Mode*toolbar to switch from*grid*mode to*physics and equation/subdomain*specification mode. - Press the
**edit eqn**button in the*Equation Settings*dialog box that automatically opens. - 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. - Enter
**x^5**in the*Initial Conditions*edit field. Press**OK**to finish with the equation specifications and close the dialog box. - Switch to
*boundary condition specification*mode by clicking on the mode button - 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.

- 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. - Now that the problem is specified, press the mode button to switch to
*solve*mode. Then press the**Solver settings**toolbar button. - Select the
**Time-Dependent - Monolithic (Fully coupled)**solver option, then press**Apply**and**Solve**to start the solution process. - 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. - In the Postprocessing settings dialog box enter
**u-(x-t)^5**to show the difference between the computed and exact analytical solution. - 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.

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.

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.

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

for a scalar field on a circle with radius in two dimensions. Both the diffusion coefficient and right hand side source term 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

Moreover, homogeneous Dirichlet boundary conditions are prescribed on all boundaries of the domain, that is . The exact solution for this problem is which can be used to measure the accuracy of the computed solution.

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.

- 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

In the command window type

featool

to start the graphical user interface (GUI).

- 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. - 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. - To create a circle, first left click on the
**Create circle/ellipse**button in the left hand side*Tools*toolbar frame. - 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.
- 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 - 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**. - Press the mode button in the
*Mode*toolbar to switch from*geometry*mode to*grid*generation mode. - Click on the
**Generate**button to call the grid generation function which automatically generates a grid of triangles for the circle. - 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 both the diffusion coefficient and source term coefficient 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. - Switch to
*boundary condition specification*mode by clicking on the mode button - 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 coefficientequal to*r***0**in the*Boundary Coefficients*frame and finish by clicking on**OK**. - 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. - 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. - Activate
*Contour*and*Arrow*plots by marking the corresponding check boxes and press**OK**or**Apply**to show the new plot. - The solution is now visualized as surface, contour, and arrow plots in the main plot axes window.
- 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. - 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.
- 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.

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.

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.

- 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 };

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

- Assign names to the space dimension coordinates (these names can then be used in equation coefficients and postprocessing expressions)
fea.sdim = { 'x' 'y' };

- Use addphys to add the predefined Poisson physics mode poisson
fea = addphys( fea, @poisson );

- 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

- 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.

- 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 zerofea.phys.poi.bdr.coef

- 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};

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

- Call the stationary solver solvestat with the default settings
fea.sol.u = solvestat(fea);

- Plot with postplot to display the scalar solution as surface and contour plots, and the solution gradient (
*ux*,*uy*) as arrow plotspostplot( fea, 'surfexpr', 'u', ... 'isoexpr', 'u', ... 'arrowexpr', {'ux','uy'} )

- Plot the difference between the computed and exact solution
u_diff = 'u-(1-x^2-y^2)/4'; clf postplot( fea, 'surfexpr', u_diff )

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).

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.

- 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 };

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

- Assign names to the space dimension coordinates (these names can then be used in equation coefficients and postprocessing expressions)
fea.sdim = { 'x' 'y' };

- Assign a name to the dependent variable (in this case
*u*)fea.dvar = { 'u' };

- 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' };

- 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 };

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

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

- Call the stationary solver solvestat with the default settings
fea.sol.u = solvestat( fea );

- Plot the results with postplot to display the scalar solution as and contour plots, and the solution gradient (
*ux*,*uy*) as arrow plotspostplot( fea, 'surfexpr', 'u', ... 'isoexpr', 'u', ... 'arrowexpr', {'ux','uy'} )

- Plot the difference between the computed and exact solution
u_diff = 'u-(1-x^2-y^2)/4'; clf postplot( fea, 'surfexpr', u_diff )

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.

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.

- 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 };

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

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

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

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

- 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.

- Solve (invert) problem by calling the built-in linear solver
u = A\f;

- Create mock fea fields for postprocessing
fea.sdim = { 'x' 'y' }; fea.dvar = { 'u' }; fea.sfun = { 'sflag1' }; fea = parseprob( fea );

- Use postplot to plot the solution
fea.sol.u = u; postplot( fea, 'surfexpr', 'u' )

- 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' )

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.

The following example illustrates how to set up and solve the Poisson equation on a unit sphere ( ) with homogeneous boundary conditions, that is with the FEATool FEM assembly primitives.

To derive the finite element discretization of the general Poisson equation first multiply the equation with an arbitrary function (called test function), and integrate over the whole domain , that is . By applying the Gauss theorem or partial integration to the left side we get

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

By discretizing the solution (trial function space) variable as and similarly test function , where are the finite element shape or basis functions (usually taken the same for and in the Galerkin approximation), and is the value of the solution at node . Inserting this we get

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

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) . 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' )

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

Field | Description |
---|---|

ex_brinkman1 | Coupled porous and incompressible laminar flow. |

ex_convdiff1 | Convection and diffusion on a unit square. |

ex_convdiff2 | Periodic convection and diffusion on a line with exact solution. |

ex_convdiff3 | Infinite time dependent convection and diffusion on a line. |

ex_convdiff4 | One dimensional Burgers equation with steady solution. |

ex_convdiff5 | Dominating convection example requiring artificial stabilization. |

ex_convdiff6 | 1D stationary convection-diffusion-reaction with exact solution. |

ex_convreact1 | Time dependent one dimensional convection and reaction model. |

ex_custom_equation1 | Black-Scholes model equation implemented as a custom equation. |

ex_euler_beam1 | 1D Euler-Bernoulli beam example. |

ex_euler_beam2 | 1D Euler-Bernoulli beam vibration mode example. |

ex_nonlinear_pde1 | 1D nonlinear PDE with analytical solution. |

ex_nonlinear_pde2 | 1D nonlinear PDE with analytical solution extended to 2D. |

ex_diffusion1 | Diffusion equation on a unit square with different solutions. |

ex_diffusion2 | Diffusion equation on a line with exact solutions. |

ex_heattransfer1 | 2D heat conduction with natural convection and radiation. |

ex_heattransfer2 | One dimensional stationary heat transfer with radiation. |

ex_heattransfer3 | One dimensional transient heat conduction. |

ex_heattransfer4 | Two dimensional heat transfer with convective cooling. |

ex_heattransfer5 | Two dimensional transient cooling shrink fitting example. |

ex_heattransfer6 | Axisymmetric steady state heat conduction of a cylinder. |

ex_heattransfer7 | One dimensional transient heat conduction with analytic solution. |

ex_laplace1 | Laplace equation on a unit square. |

ex_laplace2 | Laplace equation on a unit circle. |

ex_linearelasticity1 | Linear elasticity test case, solid stress-strain for a cube. |

ex_linearelasticity2 | Three dimensional example of stress on a bracket. |

ex_linearelasticity3 | Parametric study of the bracket deformation model. |

ex_linearelasticity4 | Stress of loaded I-beam supported by two brackets. |

ex_spanner | Gmsh import of a mesh and stress calculation in a spanner. |

ex_multiphase1 | Multiphase flow example of a static bubble. |

ex_multiphase2 | Multiphase flow example of a rising bubble. |

ex_multiphase3 | Multiphase flow example of a breaking dam. |

ex_heat_exchanger1 | Free and forced convection in a heat exchanger. |

ex_natural_convection | Natural convection in a square cavity. |

ex_darcy1 | Porous media flow in a packed bed reactor using Darcy's law. |

ex_navierstokes1 | Incompressible fluid flow in a channel. |

ex_navierstokes2 | Incompressible driven cavity flow. |

ex_navierstokes3 | Incompressible fluid flow around a cylinder in a channel (Re=20). |

ex_navierstokes3b | Flow around a cylinder solver benchmark. |

ex_navierstokes4 | Incompressible fluid flow over a backwards facing step. |

ex_navierstokes5 | Two dimensional decay of a standing vortex. |

ex_navierstokes6 | Time dependent flow around a cylinder in a channel (0<=Re<=100). |

ex_navierstokes7 | Laminar flow in a curved three dimensional pipe. |

ex_navierstokes8 | Axisymmetric flow in a circular pipe. |

ex_navierstokes8b | Axisymmetric pipe flow solver benchmark. |

ex_navierstokes9 | Axisymmetric impinging jet. |

ex_navierstokes10 | 3D flow in a pipe using the OpenFOAM solver. |

ex_navierstokes11 | 3D cavity flow in a cube using the OpenFOAM solver. |

ex_navierstokes12 | 3D flow over a backwards facing step using the OpenFOAM solver. |

ex_navierstokes13 | 3D flow past a cylinder using the OpenFOAM solver. |

ex_navierstokes14 | 2D Pousille flow due to pressure gradient. |

ex_swirl_flow1 | 2D axisymmetric swirl flow in a tube with analytic solution. |

ex_swirl_flow2 | 2D axisymmetric swirl flow in a stepped narrowing tubular region. |

ex_swirl_flow3 | 2D axisymmetric time dependent Taylor-Couette swirl flow. |

ex_swirl_flow4 | 2D axisymmetric swirl flow around a rotating disk. |

ex_piezoelectric1 | Bending of a beam due to piezoelectric effects. |

ex_electrostatics1 | Electrostatics test problem. |

ex_electrostatics2 | Axisymmetric model of a spherical capacitor. |

ex_resistive_heating1 | Resistive heating of a spiral tungsten filament. |

ex_magnetostatics1 | Magnetostatics test problem. |

ex_magnetostatics2 | Magnet field around a horseshoe magnet. |

ex_magnetostatics3 | Cylindrical magnet without electrical currents. |

ex_magnetostatics4 | Axisymmetric model of cylindrical magnet. |

ex_magnetostatics5 | 2D Magnetostatics test model. |

ex_planestress1 | Stress computation for a hole in a thin plate. |

ex_planestress2 | NAFEMS plane stress benchmark example. |

ex_planestress3 | NAFEMS plane stress benchmark of a tapered plate. |

ex_planestress4 | NAFEMS benchmark of plate deformation due to thermal stress. |

ex_planestress5 | NAFEMS benchmark for plane stress in an elliptic membrane. |

ex_planestress6 | 2D Plane stress analysis of a pressure vessel. |

ex_planestrain1 | 2D Plane strain analysis of a pressure vessel. |

ex_axistressstrain1 | Axisymmetric stress-strain of a hollow cylinder. |

ex_axistressstrain2 | Axisymmetric stress-strain of a hollow sphere. |

ex_axistressstrain3 | Axisymmetric disc with fixed edge and point load. |

ex_axistressstrain4 | Stress formed in brake disk due to heat build up. |

ex_poisson1 | Poisson equation on a line. |

ex_poisson2 | Poisson equation on a circle. |

ex_poisson3 | Poisson equation on a unit square. |

ex_poisson4 | Poisson equation on a rectangle with complex solution. |

ex_poisson5 | Poisson equation on a sphere. |

ex_poisson6 | Poisson equation on a unit cube. |

ex_poisson7 | Poisson equation on a unit circle with a point constraint. |

ex_periodic1 | Moving 1D pulse in a periodic domain. |

ex_periodic2 | 2D Periodic Poisson equation example. |

ex_robinbc1 | 1D example with a Robin boundary condition. |

ex_waveequation1 | Wave equation on a unit circle. |