WANtaroHP (FEM: Axisymmetric Stress Analysis)

Contents

• Outline of a program 'f90_fem_asnt.f90'
• Format for input data file ('csv' format)
• Format for output data file ('csv' format)
• Comparison with theoretical solution
• Comparison with 2D FEM
• Programs and sample data

Outline of a program 'f90_fem_asnt.f90'

• This is a program for Axisymmetric Stress Analysis with the function of no-tension stress 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.
• 4 noded isoparametric element with 4 Gauss points is used. 1 node has 2 degrees of freedom in axial direction and radius direction.
• Thickness of circumference direction of 1 element is 1 radian.
• Nodal forces, nodal displacements and nodal temperature changes can be given as loads. Inertia force can not be treated. If you give zero displacements to the nodes, it means completely restricted boundary.
• Only an isotropic material can be treated. However, since the tensile strength of a element is inputted as the element characteristic, when the tensile principal stress of an element exceeds the tensile strength of it, the behavior of the element turns into a behavior as a No-tension element which does not share tensile stress. If all elements have sufficiently large tensile strength, the behavior of a structure becomes a behavior as an elastic body.
• In coordinate system, z-direction is defined as axial direction and r-direction is defined as radius direction.
• Simultaneous linear equations are solved using Cholesky method for banded matrix. Although this solving process has iterative process for no-tension analysis, stiffness matrix is not changed using stress transfer method. Therefore, to create the triangular matrix is carried out only one time, and forward elimination and backward substitution process are repeated depending on the values of unbalanced forces.
• Convergence criterion is that the case ratio of displacement increment to total displacement for all degrees of freedom becomes less than 1$\times$10$^{-6}$. Upper limit of iteration is set to 2000 times.
• Input/Output file name can be defined arbitrarily with the format of 'csv' and those are inputted from command line of a terminal.
• Used language for program is 'Fortran 90' and used compiler is 'GNU gfortran.'

Format for input data file ('csv' format)

comment                           # Comment
NODT,NELT,MATEL,KOZ,KOR,NF,IPR    # Basic values for analysis
Em,po,gamma,gkz,alpha,ts          # Material properties
      ....(1 to MATEL)....        #
node-1,node-2,node-3,node-4,matno # Element connectivity, Material set number
      ....(1 to NELT)....         # (counterclockwise order of node numbers)
z,r,deltaT                        # Node coordinates, Temperature change of node
      ....(1 to NODT)....         # (z: axial direction, r: radius direction)
nokz,rdis                         # Restricted node number and displacement in z-direction
      ....(1 to KOZ)....          # (Omit data input if KOZ=0)
nokr,rdis                         # Restricted node number and displacement in r-direction
      ....(1 to KOR)....          # (Omit data input if KOY=0)
node,fz,fr                        # Loaded node number, Load in z-direction, Load in r-direction
      ....(1 to NF)....           # (Omit data input if NF=0)

Notice

• Temperature changes are inputted as them of nodes ans they must be written in the same row written the node coordinates. For this item, temperature rising is positive.
• When tensile strength of element is not large and temperature decrease is large, equilibrium state may not be obtained because of indeterminate of structure.
• Restricted node means the node which has known (given) displacement. As a known (given) value of nodal displacement, any value can be given including zero for a restricted node.

Load input method

We consider the cylinder under the internal pressure. The cylinder has the radius of r=3000mm, the length in axial direction of z=200mm and internal pressure of p=1 N/mm$^2$ is acted. Total load acted to 1 element with thickness of 1 radian is shown below.

So, below loads must be acted for each node according to the thinking of equivalent nodal force,

Format for output data file ('csv' format)

Comment
NODT,NELT,MATEL,KOZ,KOR,NF,IPR
(Each value for above items)
*node characteristics
node,z,r,fz,fr,fix-z,fix-r,rdis-z,rdis-r,deltaT
     node   Node number
     z,r    : z & r-coordinates
     fz,fr  : z & r-direction nodal forces
     fix-z  : z-direction restricted condition (1: restricted, 0: not restricted)
     fix-r  : r-direction restricted condition (1: restricted, 0: not restricted)
     rdis-z : z-direction specified displacement
     rdis-r : r-direction specified displacement
     deltaT : Temperature change of node
  .....(1 to NODT).....
*element characteristics
element,node-1,node-2,node-3,node-4,E,po,gamma,gkz,alpha,ts,matno
     element : Element number
     node-1,node-2,node-3,node-4 : Element-nodes relationship
     E       : Elastic modulus of element
     po      : Poisson's ratio of element
     gamma   : Unit weight of element
     gkz     : Acceleration in z-direction (ratio to 'g')
     alpha   : Coefficient of thermal expansion of element
     ts      : Tensile strength of element
     matno   : material set number
  .....(1 to NELT).....
*displacement and force
node,coord-z,coord-r,dist-z,dist-r,reac-z,reac-r,ftvec-z,ftvec-r
     node            : Node number
     coord-z,coord-r : Coordinates of z & r-direction
     dist-z,dist-r   : Displacement in z & r-direction
     reac-z,reac-r   : Internal force in z & r direction
     ftvec-z,ftvec-r : Unbalanced force in z & r-direction
  .....(1 to NODT).....
*stresses
element,kk,sig-z,sig-r,sig-t,tau-zr,ps1,ps2,ang,noten,matno
         element     : Element number
         kk          : Gauss point number (IPR=0: 1 to 4, IPR=1: Average stress of element)
         sig-z,sig-r,sig-t,tau-zr : Normal stresses in z & r & $\theta$-direction, shearing stress
         ps1,ps2,ang : 1st & 2nd principal stresses, principal direction (degree)
         noten       : flag for No-tension (0: Elastic, 1: ps1=no-tension, 2: ps1 and ps2=no-tension)
         matno       : Material set number
  .....(1 to NELT, depending on the value of IPR).....
NODT=(Number of nodes), nt=(nt), mm=(mm), ib=(ib)
nnn=(Number of itaration), icount=(Number of converged degrees of freedom)
Calculation time=(calculation time)
Date_time=(date of execution)
    nt : Total degrees of freedom of FE equation
    mm : Dimension of reduced FE equation
    ib : band width of reduced FE equation

Comparison with theoretical solution

Consider the circular cylindrical shell with thick wall under uniform internal pressure and uniform external pressure in plane strain state. Material has the properties of elastic modulus $E=25,000N/mm^2$, Poisson' ratio $\nu=0.2$. $\sigma_r$, $\sigma_{\theta}$ and $u$ mean normal stress in radius direction, normal stress in circumference direction and displacement in radius direction for each.

\begin{equation*} \begin{cases} u=\cfrac{a^2}{b^2-a^2}\left(\cfrac{(1+\nu)(1-2\nu)}{E}\cdot r+\cfrac{1+\nu}{E}\cdot\cfrac{b^2}{r}\right)\cdot P_a \\ ~~~~-\cfrac{b^2}{b^2-a^2}\left(\cfrac{(1+\nu)(1-2\nu)}{E}\cdot r+\cfrac{1+\nu}{E}\cdot\cfrac{a^2}{r}\right)\cdot P_b \\ \sigma_r=\cfrac{a^2}{b^2-a^2}\left(1-\cfrac{b^2}{r^2}\right)\cdot P_a-\cfrac{b^2}{b^2-a^2}\left(1-\cfrac{a^2}{r^2}\right)\cdot P_b \\ \sigma_{\theta}=\cfrac{a^2}{b^2-a^2}\left(1+\cfrac{b^2}{r^2}\right)\cdot P_a-\cfrac{b^2}{b^2-a^2}\left(1+\cfrac{a^2}{r^2}\right)\cdot P_b \end{cases} \end{equation*}

• In axisymmetric FEM model, number of division in radius direction is 5, and number of division is 1. Displacements in axial direction are restricted in order to satisfy the condition of plane strain state.
• In theoretical solution, normal stresses at inner surface and outer surface of wall are obtained under strict boundary conditions. The other hand, in axisymmetric FEM, normal stresses shown above tables are mean values of inner side element or outer side element. So, the difference is found between theoretical solution and results of axisymmetric FEM. But the difference is not so large.
• In displacement, the difference between theoretical solution and axisymmetric FEM result is small and it has the accuracy of 0.1%.

Comparison with 2D FEM

Model

• By axisymmetric FEM and 2D FEM, the comparison of results were performed.
• Characteristics of model for analysis is shown in following table.
• In axisymmetric FEM model, number of division in radius direction is the same value as 2D FEM. (number of division is 33)
• Stress state was plane strain.
• The comparable result was obtained by axisymmetric FEM and 2D FEM.

Section Double reinforced Internal diameter of tunnel 4,000 mm Thickness of concrete wall 800 mm Outer diameter of domain 50,000 mm Cover for Re-bar 100 mm Inner re-bar (equivalent steel thickness) D32@200x2 3.97 mm Outer re-bar (equivalent steel thickness) D32@200x2 3.97 mm Elastic modulus of concrete 25,000 N/mm$^2$ Poisson's ratio of concrete 0.2 ( zero at No-tension state) Thermal expansion coefficient of concrete 10$\times$10$^{-5}$ $^\circ$C$^{-1}$ Elastic modulus of re-bar 200,000 N/mm$^2$ Poisson's ratio of re-bar 0.3 Thermal expansion coefficient of re-bar 10$\times$10$^{-6}$ $^\circ$C$^{-1}$ Elastic modulus of bed rock 1 to 100,000 N/mm$^2$ Poisson's ratio of bed rock 0.25 Thermal expansion coefficient of bed rock 7$\times$10$^{-6}$ $^\circ$C$^{-1}$ Internal water pressure 1 MPa Temperature change in the domain (by heat conduction analysis) (see right figure) ($-$10 $^\circ$C at inner surface)

Results

Axisymmetric FEM model	2D plane strain model
Model for analysis of whole area (wall thickness: 800mm)

Axisymmetric FEM model	2D plane strain model}
Model for analysis around concrete wall (wall thickness :800mm, double reinforcement bars)

Programs and sample data

Shell script and f90 programs

Input data smple for 'f90_asnt' (thick cylindrical shell model)

Input data sample for 'f90_asnt' (reinforcing bar model)

Input data sample for 'f90_plnt' (reinforcing bar model: for comparison)