Outline of this page
In this page, typical FEM programs by Fortran90 are introduced.
- Simultaneous linear equations are solved using Cholesky method for banded matrix .
- 3 noded linear triangular element is used as triangular element.
- 4 noded isoparametric element with 4 Gauss points is used as quadrilateral element.
- Programs shown in this page can not solve the mixed model problem of 3 node and 4 node elements.
You need to choose either a triangular or a quadrilateral element, before solving a problem.
- Optimization programs of node number are prepared for 2D-Heat Conduction, 2D-Seepage Flow and 2D-Stress Analysis programs.
The purpose of these programs is to reduce the calculation time by reducing band width of stiffness matrix.
Another webpage (http://www.archi.hiro.kindai.ac.jp/laboratory/SAL/dfujii/Report/fem/fem_5.pdf) was referred for main subroutine of the optimization program of node number in this page.
2D-Nonsteady Heat Conduction Analysis
Outline of the program
- This is a program for 2D-nonsteady heat conduction analysis. (source code: f90_FEM_HEAT.f90)
- 4 noded isoparametric element with 4 Gauss points is used. 1 node has 1 degree of freedom.
- The thickness of element is constant. It cannot be changed.
- Following conditions are available.
Thermal insulation boundary | No input means thermal insulation boundary. |
Temperature given boundary | Specify the nodes and time histories of temperature at specified nodes. |
Heat transfer boundary | Specify the sides of element, heat transfer coefficients and time histories of temperature at outside. |
Heating material | Heating material can be used like a cement concrete. |
- Formula for heating material
T=K*[1-exp(-α*t)]
T: Temperature rise, K: Maximum Temperature rise, t: time, α: coefficient of heat release rate
- Crank-Nicolson method is used for the discretization of the time.
- If you use the optimization program of node number (source code: f90_NUMheat.f90), band width of matrix and calculation time can be reduced.
- Batch commands for command prompt of Windows are shown below.
001 gfortran -o f90_NUMheat.exe f90_NUMheat.f90
002 gfortran -o f90_FEM_HEAT.exe f90_FEM_HEAT.f90
003
004 f90_NUMheat inp_org_heat.csv inp_ren_heat.csv
005 f90_FEM_HEAT inp_ren_heat.csv out_ren_heat.csv
- 001
- To compile the program of 'f90_NUMheat' (optimization program of node number)
- 002
- To compile the program of 'f90_FEM_HEAT' (FEM program for heat conduction analysis)
- 004
- To execute the program of 'f90_NUMheat' (optimization program of node number), Input-file: inp_org_heat.csv, Output-file: inp_ren_heat.csv
- 005
- To execute the program of 'f90_FEM_HEAT' (FEM program for heat conduction analysis), Input-file: inp_ren_heat.csv, Output-file: out_ren_heat.csv
Source code and related documents
2D-Steady Seepage Flow Analysis (saturated-unsaturated analysis)
Outline of the program
- This is a program for 2D-Steady Seepage Flow Analysis. (source code: f90_FEM_SEEP.f90)
- Linear triangular element or 4 noded isoparametric element with 4 Gauss points can be used. 1 node has 1 degree of freedom.
This program can not solve the mixed model problem of triangular elements and quadrilateral elements. You need to choose either triangular model or quadrilateral model.
- This program can do the saturated-unsaturated seepage flow analysis in the vertical section, and saturated seepage flow analysis in the horizontal section.
- Von Genuchten model is used as a unsaturated ground model.
- Following conditions are available.
Impervious boundary | No-input means impervious boundary. |
Total head given boundary | Specify the nodes and total heads of the nodes. |
Discharge given boundary | Specify the nodes and discharges. Positive value of discharge means in-flow, and negative value of discharge means out-flow. |
Saturation point boundary | Specify the nodes which have the possibility of becoming saturation point. Pressure head at these nodes are zero or negative value and signs of discharge are negative (out-flow). |
Section for analysis | Choose either vertical or horizontal section. Pressure head in vertical section is given as the difference between total head and position head, and pressure head in horizontal section is same as total head in it. |
- If you use the optimization program of node number (source code: f90_NUMseep.f90), band width of matrix and calculation time can be reduced.
- Batch commands for command prompt of windows are shown below.
001 gfortran -o f90_NUMseep.exe f90_NUMseep.f90
002 gfortran -o f90_FEM_SEEP.exe f90_FEM_SEEP.f90
003
004 f90_NUMseep inp_org_seep.csv inp_ren_seep.csv
005 f90_FEM_HEAT inp_ren_seep.csv out_ren_seep.csv
- 001
- To compile the program of 'f90_NUMseep' (optimization program of node number)
- 002
- To compile the program of 'f90_FEM_SEEP' (FEM program for seepage flow analysis)
- 004
- To execute the program of 'f90_NUMseep' (optimization program of node number), Input-file: inp_org_seep.csv, Output-file: inp_ren_seep.csv
- 005
- To execute the program of 'f90_FEM_SEEP' (FEM program for seepage flow analysis), Input-file: inp_ren_seep.csv, Output-file: out_ren_seep.csv
Source code and related documents
2D-Stress Analysis (no-tension analysis)
Outline of the program
- This is a program for 2D-stress analysis. (source code: f90_FEM_PLNT.f90)
- Linear triangular element or 4 noded isoparametric element with 4 Gauss points can be used. 1 node has 2 degrees of freedom.
This program can not solve the mixed model problem of triangular elements and quadrilateral elements. You need to choose either triangular model or quadrilateral model.
- This program can do plane stress analysis or plane strain analysis.
- This program can also do no-tension analysis .
- As a method to simulate no-tension state, stress transfer method by Zienkiewicz is used.
In this method, equilibrium state can be obtained by iterative calculation method using linear stiffness matrix and coordinate transformation method of internal forces. It is assumed that Poisson ratio in no-tension state element is zero.
- Following conditions are available.
Load | Specify the loaded node and load value. |
Temperature | Specify the temperature change for all nodes. |
Inertia force | Specify the inertia force as a ratio to acceleration of gravity in horizontal and vertical direction for all elements. |
Displacement | Specify the nodes and displacements at the nodes. Any value can be inputted including zero. |
No-tension material | Specify the tensile strength for all elements.
If the value of maximum principal stress exceeds the tensile strength of an element, it is assumed that the element becomes no-tension material. |
- If you use the optimization program of node number (source code: f90_NUMseep.f90), band width of matrix and calculation time can be reduced.
- Since the program 'f90_NUMseep.exe' calls 'f90_num.exe' in it, it is necessary to exist program 'f90_num.exe' in the same directory as the program 'f90_NUMseep.exe.'
- Batch commands for command prompt of Windows are shown below.
001 gfortran -o f90_num.exe f90_num.f90
002 gfortran -o f90_NUMplnt.exe f90_NUMplnt.f90
003 gfortran -o f90_FEM_PLNT.exe f90_FEM_PLNT.f90
004
005 f90_NUMplnt inp_org_plnt.csv inp_ren_plnt.csv
006 f90_FEM_PLNT inp_ren_plnt.csv out_ren_plnt.csv
- 001
- To compile the program 'f90_num' (Main part of optimization program of node number)
- (If the program 'f90_num.exe' exists in the same directory as the program 'f90_NUMplnt.exe,'
it is called by 'f90_NUMplnt.exe' automatically.)
- 002
- To compile the program 'f90_NUMplnt' (Control part of optimization program of node number)
- 003
- To compile the program 'f90_FEM_PLNT' (FEM program for 2D stress analysis)
- 004
- To execute the program 'f90_NUMplnt' (optimization of node number), Input-file: inp_org_plnt.csv, Output-file: inp_ren_plnt.csv
- 006
- To execute the program 'f90_FEM_PLNT' (FEM program for 2D stress analysis), Input-file: inp_ren_plnt.csv, Output-file: out_ren_plnt.csv
Source code and related documents
Axisymmetric Stress Analysis (no-tension analysis)
Outline of the program
- This is a program for Axisymmetric Stress Analysis. (source code: f90_FEM_ASNT.f90)
- 4 noded isoparametric element with 4 Gauss points is used. 1 node has 2 degrees of freedom.
- This program can do no-tension analysis .
- As a method to simulate no-tension state, stress transfer method by Zienkiewicz is used.
- Following conditions are available. Inertia force is not available in this program.
Load | Specify the loaded node and load value. |
Temperature | Specify the temperature change for all nodes. |
Displacement | Specify the nodes and displacement at the node. Any value can be inputted including zero. |
No-tension material | Specify the tensile strength for all elements. If the value of maximum principal stress exceeds the tensile strength of an element, it is assumed that the element becomes no-tension material. |
- Batch commands for command prompt of windows are shown below.
001 gfortran -o f90_FEM_ASNT.exe f90_FEM_ASNT.f90
002
003 f90_FEM_ASNT inp_org_asnt.csv out_org_asnt.csv
- 001
- To compile the program 'f90_FEM_ASNT' (FEM program for axisymmetric stress analysis)
- 003
- To execute the program 'f90_FEM_ASNT' (FEM program for axisymmetric stress analysis), Input-file: inp_org_asnt.csv, Output-file: out_org_asnt.csv
Source code and related documents
2D-Frame Analysis
Outline of the program
- This is a program for 2D-Frame Analysis. (Source code: f90_FEM_FRAME.f90)
- This program can do the linear elastic analysis of frame using beam element. Beam element has 2 nodes for 1 element and 3 degrees of freedom for 1 node.
- Following conditions are available.
Load | Specify the loaded nodes and load values. |
Temperature | Specify the temperature change for all nodes. |
Inertia force | Specify the inertia force as a ratio to acceleration of gravity in horizontal and vertical direction for all elements. |
Displacement | Specify the nodes and displacements at the nodes. Any value can be inputted including zero. |
- Batch commands for command prompt of windows are shown below.
001 gfortran -o f90_FEM_FRAME.exe f90_FEM_FRAME.f90
002
003 f90_FEM_FRAME inp_org_frame.csv out_org_frame.csv
- 001
- To compile the program 'f90_FEM_FRAME' (FEM program for 2D frame analysis)
- 003
- To execute the program 'f90_FEM_FRAME' (FEM program for 2D frame analysis), Input-file: inp_org_frame.csv, Output-file: out_org_frame.csv
Source code and related documents
2D Mesh generator for FEM
Outline of the programs
- These are programs to generate triangular or quadrilateral elements for 2D FEM analysis.
- These programs have been arranged refering to the textbook by Dr. Takeo TANIGUCHI.
- Drawings in the document 'TeX_automesh.pdf' are drawn using GMT. So, support program for GMT drawing 'f90_GMT_MESH' is shown in this article.
- GMT drawing commands are included in the sample of execution bach file 'bat_a1T' and 'bat_a1Q'.
- The traiangular elements can be obtein to execute the program of 'f90_FEM_MESH3' only one time.
But sometimes you cannot get good result and it is depend on your input data.
- The quadrilateral elements can be obtain to execute the 3 programs of 'f90_FEM_MESH4A', 'f90_FEM_MESH4B' and 'f90_FEM_MESH4C'.
But sometimes you cannot get good result and it is depend on your input data.
It is difficult to obtain good shape quadrilateral elements comparison with to obtain good shape triangular elements.
- In addition, program 'f90_FEM_DOMAIN' is shown in this article. This one is to define to domain for each element.
Although domain number (material number) is defined in the program of 'f90_FEM_MESH3' for triangular elements, some users may want to re-try to define the domain for each triangular elements.
Therefore, this program was made with the purpose of definition of the domain for each triangular and quadrilateral elements.
Source code and related documents
For Trianglar elements
For Quadrilateral elements
Document and support program for GMT drawing
Filename | Description |
TeX_automesh.pdf | Document for useage of programs |
f90_GMT_MESH.txt | Support program for GMT drawing (only mesh drawing) |
f90_FEM_DOMAIN.txt | Support program to define domains for each element |
a_domain.txt | Batch file sample for execution of f90_FEM_DOMAIN including GMT commands |
out_test3T.txt | Output sample by f90_FEM_MESH3 for triangular elements |
out_test3Q.txt | Output sample by f90_FEM_MESH4C for quadrilateral elements |
inp_dom.txt | input sample of definition file of domain |
out_test3T.txt | Out put sample by f90_FEM_DOMAIN for triangular elements with domain number for each element |
out_test3Q.txt | Output sample by f90_FEM_DOMAIN for quadrilateral elements with domain number for each element |
The result drawings by 'a_domain.bat' and 'f90_FEM_DOMAIN' are shown below.