Description of FEM theory

Contents

• TeX documants including PDF
• Interpolation Function and Gauss-Legendre Quadrature
• General equilibrium equation for static structural analysis using principle of virtual work
• 2D Frame Analysis
• 2D Truss Analysis
• Grid Girder Analysis
• 3D Frame Analysis
• 2D Stress Analysis
• Axisymmetric Stress Analysis
• 2D Saturated-Unsaturated Seepage Flow Analysis
• 2D Thermal Conductivity Analysis
• 2D Frame Vibration Analysis
• 2D Frame Geometrically Nonlinear Analysis
• 2D Frame Buckling Analysis
• 1D Thermal Conductivity Analysis

TeX documents including PDF

Document files

Command for execution

Make png images from TeX pdf using ImageMagick command 'convert'

platex tex_fig_rotate.tex
dvipdfm tex_fig_rotate.dvi
convert -trim -density 400 tex_fig_rotate.pdf -bordercolor 'transparent' -border 20x20 -quality 100 fig_rotate.png

platex tex_flow.tex
dvipdfm tex_flow.dvi
convert -trim -density 400 tex_flow.pdf -bordercolor 'transparent' -border 20x20 -quality 100 fig_flow.png

Make main TeX pdf

platex en_theory.tex
dvipdfm en_theory.dvi

Interpolation Function and Gauss-Legendre Quadrature

Introduction

In the finite element method, it is necessary to express the physical quantities of any point in the element with nodal physical quantities using interpolation function $[\boldsymbol{N}]$.

For exmple, in the thermal conductivity analysis, a temperature at any point $T$ can be expressed as follow.

where, $\phi_i$, $\phi_j$, $\phi_k$, $\phi_l$ mean the temperature at node $i$, $j$, $k$, $l$ for each.

If vectoe expression is applied,

where, $\{\boldsymbol{\phi}\}$ means nodal temperature vector.

In this section, it is considered to use the Gauss-Legendre quadrature for 4 nodes element and to express the physical quantities using the normalized parameter coordinate $(a,b)$ with range of $[−1,1]$.

Interpolation function and its derivative

In actural calculation, since a mathematical expression and derivative of interpolation function are required, they will be derived.

After this, the interpolation function $N(a, b)$ is used, where $(a,b)$ is normalized parameter coordinate with range of $a=[−1,1]$ and $b=[−1,1]$.

Mathematical expression of interpolation function

For a 4 nodes element, following interpolation function can be applied.

Derivative of interpolation function

The partial derivatives of $N$ about $a$ and $b$ can be obtained s follows.

Next, above equations can be expressed using Jacobi matrix $[J]$ as follows.

From above,

the each element of $[J]$ can be shown as follows, where $x_{i,j,k,l}$ and $y_{i,j,k,l}$ means nodal coordinates which form a finite element.

From above, the elements of $[J]$ were defined and derivative of $[\boldsymbol{N}]$ can be calculated.

At this point, it is necessary to note that $[J]$, $[\boldsymbol{N}]$ and its derivative are function of variables $a$ and $b$.

Gauss-Legendre quadrature

Gauss-Legendre quadrature is used for integration.

Surface integral

When 4 Gauss points ($n=2$) are considered in an finite element, the values of $a$, $b$ and weight $H$ can be shown in below table. In this case, the approximation of integral value can be calculated as a summation of 4 times calculation depending on the coordinate $(a, b)$. Furthermore, since the weight $H$ is equal to $1$ for all coordinate, the calculation can be more simplified.

Curvilinear integral

In the thermal conductivity analysis, curvilinear integral is requied on the element side which has heat transfer boundary condition. The method of curvilinear integral along the element side is shown below.

General equilibrium equation for static structural analysis using principle of virtual work

In this section, a derivation of general equilibrium equation for static structural analysis using principle of virtual work is descrived.

Generally, in case that a elastic body subjected to the external force is under the equilibrium state, the virtual work due to the internal stresses in the body is equal to the virtual work due to the surface force or body force. This relationship can be expressed as follow.

At this point, it is assumed that a body deems an finite element and strain $\{\boldsymbol{\epsilon}\}$ and displacement $\{\boldsymbol{u}\}$ at any point in a fifite element can be expressed using nodal displacement $\{\boldsymbol{u}\}$ as follows.

where, $[\boldsymbol{B}]$ is a strain-displacement relationship matrix, $[\boldsymbol{N}]$ is an interpolation function matrix.

And a stress-strain relationship of elastic body is applied shown below.

where, $[\boldsymbol{D_e}]$ is stress-strain relationship matrix, $\{\boldsymbol{\epsilon_0}\}$ is initial strain due to temperature change and so on.

Next, it is assumed that $\{\boldsymbol{F}\}$ can be expressed using nodal inertia force $\{\boldsymbol{w_{nd}}\}$ as follow.

From above, left side of virtual work equation which means virtual work due to internal force becomes as follow.

And right side of virtual work equation which means virtual work due to external force becomes as follow.

As a result, general form of stiffness equation can be obtained shown below.

Above equation is to calculate the displacement of the structure. To calculate the stress of the element, following equation above mentioned shall be used.

where, $\{\boldsymbol{\epsilon}\}$ means the strain vector of an element calculated from the displacements which are obtained by solving a stiffness equation.

2D Frame Analysis

Element stiffnes equation

Element stiffness equation is shown below.

Element displacements and forces	local coordinate system (x-y) and global coordinate system (X-Y)

Coordinate transformation matrix

Coordinate transformation matrix from global coordinate system to local coordinate system is shown below.

Stiffness equation in global coordinate system

Stiffness equation in global coordinate system is shown below. This shall be assemblied for all elements.

Since element nodal thermal load vector $\{\boldsymbol{f_t}\}$ is a correction term of axial force of member, the coordinate transformation is requied when assemblied. Whereas, since element nodal inertia force vector $\{\boldsymbol{f_b}\}$ is defined using accelerations in global coordinate system, the coordinate transformation is not required when assemblied.

Furthermore, when section force is calculated from the solution of total stiffness equation, axial force shall be corrected as follows.

2D Truss Analysis

Element stiffness equation

The element stiffness matrix of 2D truss consists of the matrix which doesn't include the item of rotation or moment of 2D frame element stiffness matrix. The coordinate transformation matrix also doesn't include the item of rotation of that of 2D frame.

The element stiffness equation of 2D truss is shown below.

The coordinate transformation matrix of 2D truss is shown below.

Grid Girder Analysis

Element stiffness equation

The relationships of section forces between grid girder structure and 2D frame structure are shown below.

The element stiffness equation of grid girder structure is shown below.

Since the coordinate transformation is carried out on X-Y plane, the coordinate transformation matrix is the same as that for 2D frame analysis.

Global coordinate system	Local coordinate system

(Reference) Torsional constant of rectangular cross section $J$

The torsional constant $J$ such as for concrete rectangular cross section can be calculated as follow.

When the coefficient $\eta$ is introduced, the relationship between coefficient $\eta$ and the member dimension $h/b$ becomes shown below.

Above calculation was done by Python program shown below.

import numpy as np

### Result ###
# h/b  1.0    2.0    3.0    5.0    10.0   20.0
# eta  7.114  4.373  3.798  3.433  3.202  3.098
# J=h*b**3/eta

n=100
b=1.0
#hh=np.arange(1.0,21.0,1.0)
hh=np.array([1.0,2.0,3.0,5.0,10.0,20.0])
for idx,elem in enumerate(hh):
    h=elem
    s=0.0
    for i in range(1,n+1):
        s=s+1/(2*i-1)**5*np.tanh((2*i-1)*np.pi*h/2/b)
    x=1/3*(1-192/np.pi**5*b/h*s)
    print('{0:10.3f} {1:10.3f}'.format(h/b,1/x))

3D Frame Analysis

Element stiffness equation

Nodal displacement vector

Nodal force vector

Element stiffness matrix

Transformation matrix

Direction cosine of element axis (x-axis)

• An beam element has 2 nodes which are node-$i$ and node-$j$.
• The coordinate of node-$i$ is $(x_i,y_i,z_i)$, and the coordinate of node-$j$ is $(x_j,y_j,z_j)$.
• A direction cosine of an element axis which is the same as x-axis of an element can be expressed as following equations.

Submatrix of transformation matrix

$\boldsymbol{[t]}$ is a submatrix (3 x 3) of transformation matrix (12 x 12) which transforms the coordinate from global coordinate system to local coordinate system.

where, $\boldsymbol{\{x\}}$ means a vector in local coordinate system, and $\boldsymbol{\{X\}}$ means a vector in global coordinate system.

(1) Case that x-axis in local coordinate system is not parallel to Z-axis in global coordinate system

(2) Case that x-axis in local coordinate system is parallel to Z-axis in global coordinate system

In this case, since $x_j-x_i=0$，$y_j-y_i=0$ and $l=0$，$m=0$, it is necessary to re-write $\boldsymbol{[t]}$.

Full matrix of transformation matrix

Local coordinate system (x-y-z) with Zero chord angle in global coordinate system (X-Y-Z)

2D Stress Analysis

Finite element equation formulation

Calculation formula for nodal displacements

Calculation formula for element stresses

In the calculation formula for nodal displacements, the thermal load vector $\{\boldsymbol{f_t}\}$ is included in the items of load vector. However, in the element stresses calculation, the initial strain due to temperature change shall be subtracted from the strain which are calculated from the nodal displacements shown in above formula.

Equation of relationship between nodal displacement and strain at any points

Equation of nodal displacement and displacement at any point

Stress-strain relationship on 2D elastic problem

AS well known as Hooke's law, stress-strain relationship for 3D isotropic elastic body can be expressed as follow, where $E$ is elastic modulus, $\nu$ is Poisson's ratio, $\alpha$ is thermal expansion coefficient, $T$ is temperature change.

When considering $x-y$ plane, in the plane stress state,

Therefore,

In the plane strain state,

Therefore,

The rearranged result of above can be shown as follows.

Stress and strain components in 2D elastic problem

where, tensile stress and tensile strain have positive sign, and temperature increase has positive sign.

Stress-strain relationship in 2D elastic problem

Formulation as isoparametric element with 4 nodes 4 Gouss points

Introduction of strain-nodal displacement relationship matrix $[\boldsymbol{B}]$

The displacements $u$,$v$ at any point in a quadrilateral element are assumed as follows, where coordinate ($a$,$b$) is normalized parameter coordinate with range of [$-1$,$1$] for each of $a$ or $b$, $u_{i,j,k,l}$ and $v_{i,j,k,l}$ are nodal displacements which forms an element.

When the displacements at any point in an element is defined as $\{\boldsymbol{u}\}$, and the nodal displacements is defined as $\{\boldsymbol{u_{nd}}\}$, following expression can be obtained using matrix expression.

The strainat any point in an element $\{\boldsymbol{\epsilon}\}$ can be expressed using nodal displacement $\{\boldsymbol{u_{nd}}\}$.

Element stiffness matrix

A stiffness matrix of 4 nodes isoparametric element $[\boldsymbol{k}]$ can be expressed as follow using the constant element thickness $t$ and stress-strain relationship matrix $[\boldsymbol{D_e}]$.

Nodal force vector due to temperature change

Nodal force vector {\boldsymbol{f_t}\}$ due to thermal strain$\{\boldsymbol{\epsilon_0}\}$ can be expressed as follow.

The component of strain due to temperature change is shown below.

where, $\alpha$ is thermal expansion coefficient, $\nu$ is Poisson's ratio, $T_p$ is temperature change at any point in an element, $\{\phi_i~~\phi_j~~\phi_k~~\phi_l\}^T$ is nodal temperature change vector, $[N_1~~N_2~~N_3~~N_4]$ is elements of interpolation function $[\boldsymbol{N}]$. And temperature increase has positive sign.

Nodal force vector due to body force

Body force (inertia force) at any point in a element $\{\boldsymbol{F}\}=\{F_{Bx},F_{By}\}^T$ can be expressed using nodal force vector $\{w\}$ and interpolation function$[\boldsymbol{N}]$.

where, $\gamma$ is element unit weight (mass $\times$ gravity acceleration), $k_h$ and $k_v$ are horiaontal and vertical acceleration (ratio to 'g').

From above, nodal force vector due to body force (inertia force) $\{\boldsymbol{f_b}\}$ becomes can be calculated as follow.

The vector $\{\boldsymbol{w}\}$ is defined as follow by pulling out $\gamma$.

Axisymmetric Stress Analysis

Stress-strain relationship and strain-displacement relationship in axisymmetric problem

Regarding the coordinate, it is defined that the direction of rotation axis is set to z-direction, the radius direction is set to r-direction and circumferential direction is set to $\theta$-direction. Using this coordinate system, stress-strain relationshipin axisymmetric problem is shown below.

where $E$ is elastic modulus, $\nu$ is Poisson's ratio, $\alpha$ is thermal expansion coefficient, $T$ is temperature change.

By transforming above equations, followings can be obtained.

Using matrix and vector format, above can be expressed as follows.

Strain-displacement relationship is shown below, where $w$ is the displacement in rotation axis direction and $u$ is the displacement in radius direction.

Formulation as isoparametric element with 4 nodes 4 Gouss points

Introduction of strain-nodal displacement relationship matrix $[\boldsymbol{B}]$

The displacements in rotation axis direction $w$ and the displacement in radial direction $u$ at any point in a axisymmetric quadrilateral element are assumed as follows, where coordinate ($a$,$b$) is normalized parameter coordinate with range of [$-1$,$1$] for each of $a$ or $b$, $w_{i,j,k,l}$ and $u_{i,j,k,l}$ are nodal displacements which forms an element.

When the displacements at any point in an element is defined as $\{\boldsymbol{u}\}$, and the nodal displacements is defined as $\{\boldsymbol{u_{nd}}\}$, following expression can be obtained using matrix expression.

The strain at any point in an element $\{\boldsymbol{\epsilon}\}$ can be expressed as follow using nodal displacement $\{\boldsymbol{u_{nd}}\}$.

The value of $r$ coordinate at any point in an element is evaluated using following equation.

where, $r_i$, $r_j$, $r_k$, $r_l$ are the value of $r$ coordinates which forms an element.

Element stiffness matrix

When an axisymmetric finite element is considered, the volume of small element which is located at $r$ in radial direction from the rotation center becomes $r\times d\theta\times dr\times dz$. When integrals and load vector is evaluated with angle of 1 radian, the result of an integral variable conversion for Gauss-Legendre quadrature can be shown as follow.

As a result, element stiffness matrix of axisymmetric element $[\boldsymbol{k}]$ can be expressed as follow using stress-strain relationship matrix $[\boldsymbol{D_e}]$.

2D Saturated-Unsaturated Seepage Flow Analysis

Finite element equation of 2D steady saturated seepage flow problem

Governing equation of steady saturated seepage flow problem with 2D orthotropic material followed Darcy's law is shown below in $x-y$ rectangular coordinate system, where the depth is not considered.

To define the unknown $\phi$ which satisfies following weak formulation of governing equation for any $\delta \phi$ is considered using Galerkin method.

When the total head $\phi$ at any point can be expressed as follow using interpolation function matrix $[\boldsymbol{N}]$ and nodal total head vector $\{\boldsymbol{h}\}$,

the weak formulation of governing equation can be expressed as follow.

At this time, the surface integral value shall be zero, because $\{\boldsymbol{\delta h}\}$ is arbitrary　any nodal total headvector.

Next, using a partial integration formula

and Green's theorem

following equation can be derived.

And when the velocity in $x$ and $y$-directions are defined as $v_x$ and $v_y$, following equations can be obtained using Darcy's law.

As a result, element finite element equation can be obtained as follow.

where, $\{\boldsymbol{q}\}$ is integral value of velocity vector along the side, and it is easy to understand that it means discharge.

Formulation as isoparametric element with 4 nodes 4 Gouss points

Permeability matrix

When an isoparametric element has 4 Gauss integral points, an element permeability matrix can be expressed as follow.

Regarding nodal discharge vector $\{\boldsymbol{q}\}$, it can be inputted as a equivalent nodal discharge vector directly.

Average velocity in an element

The total head at any point in an element $\phi$ can be expressed as follow using interpolation function $[\boldsymbol{N}]$ and nodal total head vector $\{\boldsymbol{h}\}$.

Therefore, the mean velocity at any point in an element $v_x$ and $v_y$ can be expressed as follows using Darcy's law.

In actural calculation, following formula can be used to obtain the average velocity in an element.

(Notice) Above discussion is for orthotropic material. However, a program introduced can treat only isotropic material.

Discussion of extension to unsaturated seepage flow analysis

Saturated seepage flow analysis

A saturated seepage flow anaalysis is a problem to solve following simultaneous linear equations.

where, $\{q\}$ is nodal discharge vector, $\{h\}$ is total head vector, $[k]$ is permeability matrix, and they have following characteristics.

• Permeability matrix $[k]$ is a constant.
• components of nodal total head vector $\{h\}$ which are corresponding to known components of nodal discharge vector $\{q\}$ are unknown variavles.
• components of nodal total head vector $\{h\}$ which are corresponding to unknown components of nodal discharge vector $\{q\}$ are known variavles.

Therefore, after treatment of the relationship between known and unknown variables, the solution can be obtained by only one time calculation.

Saturated-unsaturated seepage flow analysis

simultaneous equations in saturated-unsaturated seepage analysis have following characteristics.

• Three types boundary conditions shall be considered such as given discharge, given total head and seepage face which have the posibilities of occurance of free surface.
• Either discharge or total head becomes known variable except seepage face boundary.
• Both of discharge and total head become unknown variables on the seepage face boundary.
• However, seepage face boundary has only outflow condition, because it touchs to the air. Inaddition, pressure head on that boundary should be less than or equal to zero.
• The permeabolity matrix in the saturated area is constant. However, the permeability matrix in the unsaturated area becomes a function of suction head (negative pressure head). Therefore, simultaneous equations shall be solved by iterative calculation.

As a iterative calculation method, successive substitution method is applied. In this method, initial total head are given for all nodes, and ibtained solutions are used as input values for next iterative calculation.

Regarding the values of initial total heads, it is better to give the minimum value of z-coordinate of the model (minimum altitude in the model).

Unsaturated permeability characteristics (van Genuchten model)

As unsaturated permeability characteristics, van Genuchten model shown below is applied, because it is convenient for programing to treat contineous function as material characteristics.

2D Thermal Conductivity Analysis

Finite element equation of 2D unsteady thermal conductivity problem

Governing equation of 2D unsteady thermal conductivity problem with 2D isotropic material is shown below in $x-y$ rectangular coordinate system, where the depth is not considered.

Heat flux at boundary $q$ can be expressed as follow using Fourier's law. where $n$ is an outward normal on the boundary.

To define the unknown $T$ which satisfies following weak formulation of governing equation for any $\delta T$ is considered using Galerkin method.

When the temperature at any point $T$ can be expressed as follow using interpolation function matrix $[\boldsymbol{N}]$ and nodal temperature vector $\{\boldsymbol{\phi}\}$,

the weak formulation of governing equation can be expressed as follow.

At this time, the surface integral value shall be zero, because $\{\boldsymbol{\delta\phi}\}$ is arbitrary nodal total head vector.

Next, using a partial integration formula

and Green's theorem

and using Fourier's law, the item related to $\kappa$ becomes shown below.

And following equation can be obtained.

As a result, element finite element equation can be obtained as follow.

Explanation about boundary condition

In case that heat flux $q_0$ is given on the boundary (side of element)

for adiabatic boundary, it becomes $q_0=0$.

In case that heat transfer boundary (side of element), using heat transfer rate $\alpha_c$ and outside temperature $T_c$,

Above items shall be added to heat conduction matrix or heat flux vector.

Finite element equation considered boundary condition effects

Explanation of heating material

Treatment of Heating material such as cement concrete is considered. The adiabatic temperature rise is assumed as follow.

$K$ (Tk) and $\alpha$ (Al) in above equation will be inputted as material characteristics.

Using above, heat value $Q$ and heat rate $\dot{Q}$ can be expressed as follow.

Solution of unsteady finite element equation

Unsteady finite element equation can be solved using Crank-Nicolson method.

The finite element equation of analysis model is exoressed as follow.

where, nodal temperature vector $\{\boldsymbol{\Phi}\}$, deviation of nodal temperature vector $\{\partial\boldsymbol{\Phi}/\partial t\}$ and heat flux vector $\{\boldsymbol{F}\}$ at the time $t+\Delta t/2$ are defined as follows.

By ubstituting above relationship to the finite element equation, following equation can be obtained. As a result, unknown nodal temperature at time $t+\Delta t$ can be calculated using known temperature at time $t$.

If material considered does not have the time dependence, it is possible to calculate the temperature at next step by only one time calculatio of inverse matrix. For the calculation of inverse matrix, numpy.linalg.inv(A) in Python can be used.

Formulation as isoparametric element with 4 nodes 4 Gouss points

The matrices and vectors can be expressed as follows for an isoparametric element with 4 nodes 4 Gouss points, where $\ell$ is the length of the side with specified boundary condition of an element.

Element heat conduction matrix

Element heat capacity matrix

Heat flux vector

2D Frame Vibration Analysis

Equation of motion and Nermark's βmethod

For simplification, SDOF (single degree of freedom) system is considered.

The equation of motion of SDOF system at the time $t+\Delta t$ is shown below, where dot means a differentiation by time $t$.

The displacement $u$ and the velocity $\dot{u}$ at the time $t+\Delta t$ can be expressed as follows using Taylor series.

When setting $1/3!=\beta$ for 4th term of the right side and $1/2!=\gamma$ for 3rd term of right side, following equation van be obtained omitting the higher terms.

Next, the linearization of $\dddot{u}$ is considered as follow.

Using above relationship, following equations can be obtained.

When the application of the equations to finite element method, it is better to solve the equation with unknown displacement. So, the equations can be rearranged as follows.

Substituting above to the equation of motion, followings can be obtained.

Usually, $\gamma=1/2$ is applied, and the case of $\beta=1/4$ means the average acceleration method, the case of $\beta=1/6$ means the linear acceleration method.

The equation to be solved with matrix expression

The equations for SDOF can be changed with matrix expressed as follows by fixing $\gamma=1/2$.

Calculation formula of displacement

Calculation formulas of velocity and acceleration

After calculation of the displacement $\{\boldsymbol{u}(t+\Delta t)\}$, the velocity $\{\boldsymbol{\dot{u}}(t+\Delta t)\}$ ane the acceleration $\{\boldsymbol{\ddot{u}}(t+\Delta t)\}$ can be calculated using below equations.

where, $\beta=1/4$ means average acceleration method, $\beta=1/6$ means linear acceleration method.

Stiffness matrix and mass matrix

The coordinate transformation matrix $[\boldsymbol{T}]$ transforms the displacements or forces in global coordinate system to those in local coordinate system. The relationship between the stiffness matrix in global coordinate system $[\boldsymbol{K}]$ and that in local coordinate system $[\boldsymbol{k}]$ is shown below.

The mass matrix in global coordinate system $[\boldsymbol{M}]$ and that in local coordinate system $[\boldsymbol{m}]$ have the same relationship as the stiffness matrices shown below.

Damping matrix

In the actual calculation, a damping matrix $[\boldsymbol{C}]$ should be defined to consider the effect of damping of the structure. In tis section, definition of the damping matrix is discussed.

For simplification, the equation of motion of SDOF system is considered.

Hear, Rayleigh damping is considered. In the theorem of Rayleigh damping, damping coefficient is poportional to the mass and stiffnes coefficient. This relationship can be expressed as follow.

Considering above, dampong ratio $h$ can be rearranged as follow, where $\omega$ is fundamental circular frequency, $f$ is natural frequency.

Accordingly, if damping rates $h_i$, $h_j$ and natural frequency $f_i$, $f_j$ of $i$th mode and $j$the mode, proportional constants $\zeta_m$ and $\zeta_k$ can be obtained by satisfying above equations simultaneously as shown below.

The natrix expression of damping coefficient is shown below using tow scalar constants $\zeta_m$ and $\zeta_k$.

As a method to define the constants of Rayleigh damping, following method can be used.

• To set damping rate as a popular value such as $h_1=h_2=0.05$.
• To obtain the natural frequencies $f_1$ and $f_2$ for 1st and 2nd mode using eigenvalue analysis.
• To define $\zeta_m$ and $\zeta_k$ using two values of natural frequency ($f_1$, $f_2$) and two values of damping rate ($h_1=h_2=0.05$)

Eigenvalue analysis

The problem to obtain the fundamental circular frequency is a problem to solve the characteristic equation shown below as a generalized eigenvalue problem.

After getting the fundamental circular frequency $\omega_n$, natural frequency $f_n$ can be obtained using following equation.

2D Frame Geometrically Nonlinear Analysis

Incremental stiffness equation

As well known, an incremental stiffness equation including geometrical nonlinear terms can be expressed shown below:

Calculation formulas for Internal forces

Calculation of internal force

The increments of internal force $\{\Delta f^*\}$ can be obtained using $[K_L]$ and $\{\Delta u^*\}$. Where, $[K_L]$ is a linear stiffness matrix of an element, and $\{\Delta u^*\}$ is an incremental displacement vector excluding rigid rotation in local coordinate system.

Where, $\Delta\ell$ means the difference of the previous element length and current element length. Regarding the rotation component, it is taken as the difference of the previous rotation angle (subscript $k-1$) and current rotation angle (subscript $k$).

Elimination of rigid rotation

The method of elimination of rigid rotation is shown below:

Using the addition formula for tangent,

Elimination of rigid rotation

Using above equations, $(\tan\theta_i^*)_k$, $(\tan\theta_j^*)_k$ for currenct step and $(\tan\theta_i^*)_{k-1}$, $(\tan\theta_j^*)_{k-1}$ for previous step can be obtained. From these, roration angles excluding rigid rotation in an element can be obtained as follows.

In these steps, transformation matrix which is calculated from original (initial) coordinates of nodes should be used.

Arc-Length Method

When simplified scalar load-displacement curve is considered, following equations can be obtained reffering below figure.

Concept of Arc-Length method

Initial value of $\Delta\lambda$

If $\Delta U_R=0$ is assumed, initial value of $\Delta\lambda$ can be obtained as following equation.

Above equation can take two values, and it should be noted that a sign of $\Delta\lambda_0$ is very important in Arc-Length method.

A sign of $\Delta\lambda_0$ is defined as shown below:

• Define a displacement increment vector ${\Delta U_{-1}}$ from previous equilibrium point to current equilibrium point.
• Calculate $\{\Delta U_0\}$ as $\{\Delta U_0\}=[K_T]^{-1}\{\Delta F\}$.
• Calculate an inner product $\{\Delta U_{-1}\}^T\{\Delta U_0\}=|\Delta U_{-1}|\cdot |\Delta U_0|\cdot \cos\theta$, where $\theta$ is an angle between 2 displacement increment vectors.
• If an inner product $\{\Delta U_{-1}\}^T\{\Delta U_0\} \geqq 0$, the angle $\theta$ is less than or equal to 90 degree. In this case, $\Delta\lambda_0$ has positive sign.
• If an inner product $\{\Delta U_{-1}\}^T\{\Delta U_0\} \lt 0$, the angle $\theta$ is greater than 90 degree. In this case, $\Delta\lambda_0$ has negative sign.

Regarding the scaling parameter $\phi$, it can be obtained as following equation. In this program, recommended value of $\alpha$ is one ($\alpha=1.0$). \begin{equation*} \phi=\sqrt{\frac{\alpha}{\{\Delta F\}^T\{\Delta F\}}} \end{equation*}

Regarding the arc length $\Delta s$, it can be obtained assuming $\Delta\lambda=1.0$.

Correction factor $\Delta\lambda$ for iterative calculation

Referring above conceptial figure and replacing $\Delta s$ to $\Delta L$, following equation can be obtained.

From the condition of minimization of $(\Delta L)^2$, $\Delta\lambda$ can be calculate as following equation.

Flowchart for analysis

Flowchart of Calculation

A treatment of $\Delta s$ and $\Delta\lambda$ in actual programing code is shown below:

• At the initial loading (nnn=1), $\Delta s$ is calculated and $\Delta\lambda_0$ is set to 1.0.
• If calculated $\Delta\lambda$ is much less than $\Delta\lambda_0$, the value of $\Delta s$ should be increased, because $\Delta s$ has been calculated with the assumption of $\Delta\lambda=1.0$.
• If calculated $\Delta\lambda$ is greater than initial value of $\Delta\lambda_0 (=1.0)$, the increase of the value of $\Delta s$ should be stoped, because huge value of $\Delta s$ causes unexpected behavior of the structure.

# Initial parameter setting for Arc-length method
if nnn==1:
    ds0=np.sqrt(np.sum(dis0*dis0)+spara*spara*np.sum(df*df))
    ds=ds0
    dlam0=1.0
dlam=np.sqrt(ds*ds/(np.sum(dis0*dis0)+spara*spara*np.sum(df*df)))
if np.abs(dlam) < 0.1*np.abs(dlam0): ds=ds*1.2
if np.abs(dlam0) < np.abs(dlam): ds=ds
dlam=np.sqrt(ds*ds/(np.sum(dis0*dis0)+spara*spara*np.sum(df*df)))
if np.sum(dis_ref*dis0)<0.0: dlam=-dlam

2D Frame Buckling Analysis

Characteristic equation for buckling analysis

Characteristic equation to obtain the buckling load and displacement mode of 2D frame structure is shown below.

Although the shapes of the matrices $[K_L]$ and $[K_G]$ are same as those of 2D frame geometrically nonlinear analysis, the compressive axial force $P$ included in $[K_G]$ shall have the positive sign.

The coefficient $\lambda$ can be obtained by solving the following characteristic equation. However, since $[K_G]$ includes the axial force $P$, it is necessary to carry out the one time linear analysis to obtain the axial force $P$ before eigenvalue analysis.

1D Thermal Conductivity Analysis

Solution of unsteady finite element equation

Equation to be solved by Crank-Nicolson method is shown below.

Each matrix and vector for an element can be expressed as follows.

Explanation of heating material

Treatment of Heating material such as cement concrete is considered. The adiabatic temperature rise is assumed as follow.

$K$ (Tk) and $\alpha$ (Al) in above equation will be inputted as material characteristics.

Using above, heat value $Q$ and heat rate $\dot{Q}$ can be expressed as follow.