line

| Curvature | Bending Moment | Slope EI | | --------- | -------------- | -------- | | dεf | 0 | 0 | | dεf | M0 | (dεf)e | | dεf | (dεf)p | (dεf)p |

Fig. 5.4 Moment curvature relationship for a Timoshenko beam. $$ H ^ {\prime} = \frac {d M}{(d \epsilon_ {f}) _ {p}}. $$ This can be interpreted as the slope of the strain-hardening portion of the moment-curvature curve after the removal of the elastic curvature component. Thus $$ H ^ {\prime} = \frac {d M}{d \epsilon_ {f} - (d \epsilon_ {f}) _ {e}} = \frac {(E I) _ {T}}{1 - [ (E I) _ {T} / E I ]}. \tag {5.32} $$ It is therefore possible to rewrite (5.31) as $$ d \epsilon_ {f} = \frac {d M}{E I} + \frac {d M}{H ^ {\prime}} = \frac {d M (H ^ {\prime} + E I)}{E I H ^ {\prime}} \tag {5.33} $$ and then the incremental moment–curvature relationship can be written in the form $$ d M = \frac {E I H ^ {\prime}}{(E I + H ^ {\prime})} d \epsilon_ {f}. \tag {5.34} $$ Thus during yielding the incremental stress–strain resultant relationship is $$ d M = E I \left(1 - \frac {E I}{E I \pm H ^ {\prime}}\right) d \epsilon_ {f} $$ $$ d Q = G \hat {A} d \epsilon_ {s}. \tag {5.35} $$ The shear force/shear strain relationship is always elastic whereas the moment–curvature relationship is elasto-plastic. After yielding the flexural rigidity EI is replaced by $$ E I \left(1 - \frac {E I}{E I + H ^ {\prime}}\right). $$ If the hardening parameter $H'$ is equal to zero then the material behaviour is elasto-perfectly plastic and as mentioned in Section 3.5 for elasto-plastic axial bar elements this may lead to tangential stiffness matrices which are singular. This difficulty can also be avoided by use of the initial stiffness method in which the elastic element stiffnesses are employed at every stage of the computation thereby guaranteeing a positive definite assembled stiffness matrix. # 5.4.3 Solution of nonlinear equations Let us now generate the nonlinear equilibrium equations using the virtual expression (5.11). In order to do this we require the global rather than the element expressions for the lateral displacements, rotation, curvature and shear strain. At any point in the finite element mesh the lateral displacement and rotation can be obtained from the expression $$ \left[ \begin{array}{l} w \\ \theta \end{array} \right] = N \varphi \tag {5.36} $$ where the shape function matrix is $$ N = \left[ \begin{array}{l l l l l l} N _ {1}, & 0, & N _ {2}, & 0, & \dots , & N _ {n}, & 0 \\ 0, & N _ {1}, & 0, & N _ {2}, & \dots , & 0, & N _ {n} \end{array} \right] \tag {5.37} $$ and the vector of nodal displacements is $$ \varphi = [ w _ {1}, \theta_ {1}, w _ {2}, \theta_ {2}, \dots , w _ {n}, \theta_ {n} ] ^ {T} \tag {5.38} $$ where $w_{i}$ , $\theta_{i}$ and $N_{i}$ are the lateral displacement, rotation and global shape functions associated with node i. The curvature and shear strain at any point within the entire finite element mesh is given as $$ - \frac {d \theta}{d x} = B _ {f} \varphi \quad \text { and } \quad \frac {d w}{d x} - \theta = B _ {s} \varphi \tag {5.39} $$ where $$ \boldsymbol {B} _ {f} = \left[ 0, - \frac {d N _ {1}}{d x}, 0, - \frac {d N _ {2}}{d x}, \dots , 0, - \frac {d N _ {n}}{d x} \right] \tag {5.40} $$ and $$ \boldsymbol {B} _ {s} = \left[ \frac {d N _ {1}}{d x}, - N _ {1}, \frac {d N _ {2}}{d x}, - N _ {2}, \dots , \frac {d N _ {n}}{d x}, - N _ {n} \right] \tag {5.41} $$ Virtual curvatures and shear strains are given as $$ - \frac {d (\delta \theta)}{d x} = \boldsymbol {B} _ {f} \delta \varphi \quad \text { and } \quad \frac {d (\delta w)}{d x} - \delta \theta = \boldsymbol {B} _ {s} \delta \varphi \tag {5.42} $$ respectively, where the vector of virtual nodal displacements is written as $$ \delta \varphi = [ \delta w _ {1}, \delta \theta_ {1}, \delta w _ {2}, \delta \theta_ {2}, \dots , \delta w _ {n}, \delta \theta_ {n} ] ^ {T}. \tag {5.43} $$ Thus the virtual work expression (5.11) can now be written as $$ \begin{array}{l} \int_ {0} ^ {l} [ \delta \varphi ] ^ {T} [ \boldsymbol {B} _ {f} ] ^ {T} M d x + \int_ {0} ^ {l} [ \delta \varphi ] ^ {T} [ \boldsymbol {B} _ {s} ] ^ {T} Q d x \\ - \int_ {0} ^ {l} [ \delta \varphi ] ^ {T} [ \bar {\mathbf {N}} ] ^ {T} q d x = 0 \tag {5.44} \\ \end{array} $$ where $\bar{\mathbf{N}} = [N_1, 0, N_2, 0, \dots, N_n, 0]$ . (5.45) Since (5.44) must be true for any set of virtual displacements $\delta \varphi$ then we have $$ \left\{\int_ {0} ^ {l} \left[ \boldsymbol {B} _ {f} \right] ^ {T} M d x + \int_ {0} ^ {l} \left[ \boldsymbol {B} _ {s} \right] ^ {T} Q d x \right\} - \int_ {0} ^ {l} \left[ \bar {\mathbf {N}} \right] ^ {T} q d x = 0 \tag {5.46} $$ or $$ p - f = 0. $$ In fact this equation is identical to (5.22) when there is no plasticity. Unfortunately in elasto-plastic problems M is a nonlinear function and in general we can only predict the vector p approximately. Thus (5.46) is nonlinear and since p is only approximately known than p-f will equal a residual value $\psi(\varphi)$ which we attempt to reduce to zero in our solution procedure. We evaluate contributions to p element by element and assemble in the usual manner. The contribution from element e has the form $$ \begin{array}{l} \boldsymbol {p} ^ {(e)} = \int_ {x _ {1} ^ {(e)}} ^ {x _ {2} ^ {(e)}} \left[ \begin{array}{c} 0 \\ \frac {1}{l ^ {(e)}} \\ 0 \\ - \frac {1}{l ^ {(e)}} \end{array} \right] M ^ {(e)} d x + \int_ {x _ {1} ^ {(e)}} ^ {x _ {2} ^ {(e)}} \left[ \begin{array}{c} - \frac {1}{l ^ {(e)}} \\ \frac {x ^ {(e)} - x _ {2} ^ {(e)}}{l ^ {(e)}} \\ \frac {1}{l ^ {(e)}} \\ \frac {x _ {1} ^ {(e)} - x ^ {(e)}}{l ^ {(e)}} \end{array} \right] Q ^ {(e)} d x \\ = \left[ - Q ^ {(e)}, M ^ {(e)} - \frac {(Q l) ^ {(e)}}{2}, Q ^ {(e)}, - M ^ {(e)} - \frac {(Q l) ^ {(e)}}{2} \right] ^ {T}. \tag {5.47} \\ \end{array} $$ \*The second integral evaluation is equivalent to using a 1-point Gauss rule. 

flowchart

```mermaid graph TD A["DATA"] --> B["Input data defining geometry, loading and boundary conditions, material properties, etc."] B --> C["INITAL"] C --> D["Input data for current increment."] D --> E["INCREM"] E --> F["Initialize accumulative arrays to zero. Update load vector."] F --> G["NONAL"] G --> H["Set indicator to identify type of solution algorithm."] H --> I{Is new element stiffness matrix required?} I -->|No| J["No"] I -->|Yes| K["STIFFB"] K --> L["Calculate the element stiffness matrices and store on disc."] L --> M["ASSEMB and GREDUC or RESOLV and BAKSUB"] M --> N["Assemble global stiffness matrix (or take previous one) and global load vector and solve the resulting equations for unknowns."] N --> O["REFORB"] O --> P["Calculate the residual force vector."] P --> Q["CONUND"] Q --> R["Has solution converged?"] R -->|No| S["No"] R -->|Yes| T["RESULT"] T --> U["Output the results."] U --> V["END"] V --> W["LOAD INCREMENT LOOP"] W --> X["ITERATION LOOP"] ```

Fig. 5.5 Overall structure of program TIMOSH. Note that the appropriate value of bending moment $M^{(e)}$ is inserted in (5.47). Table 5.1 shows the complete sequence of nonlinear equation solving which is very similar to the one adopted for the axially-loaded bars in Chapter 3. 1. Begin load increment. Set $f = f + \Delta f$ , iteration counter i = 0 and $\psi^{i} = \Delta f + \psi$ (that is, include equilibrium correction from previous increment). 2. Evaluate the new tangential stiffness matrix $K_{T}$ if necessary. 3. Solve $\psi^i = \mathbf{K}_T\Delta \varphi^i$ 4. Evaluate $\varphi = \varphi + \Delta \varphi^i$ . 5. For each element evaluate $M^{(e)}$ and $Q^{(e)}$ . Check $M^{(e)}$ and adjust its value accordingly to account for any plastic behaviour. Evaluate the element residual force vector $[\psi^{(e)}]^{i+1} = \mathfrak{p}^{(e)} - \mathbf{f}^{(e)}$ and assemble into the global residual force vector $\psi^{i+1}$ . 6. Check $\Delta\varphi^{i}$ for convergence. 7. If solution has converged set $\psi = \psi^{i+1}$ and go to step 1, otherwise set $i = i+1$ and go to step 2. Table 5.1 Solution procedure for elasto-plastic nonlayered Timoshenko beam analysis. # 5.4.4 Overall program structure of TIMOSH A modular approach is adopted for program TIMOSH. In fact the overall structure is identical to the structure in the programs of Chapter 3. Figure 5.5 shows the overall structure of TIMOSH. Routines DATA, INITIAL, INCREM, NONAL, ASSEMB, GREDUC, BAKSUB, CONUND, RESOLV and RESULT have already been described in Chapter 3. The only new routines are STIFFB, REFORB and, of course, the MASTER routine BEAM. 5.4.5 New routines for nonlayered elasto-plastic Timoshenko beam analysis Master BEAM The master calling routine BEAM simply organises the calling of the main routines as described in Fig. 5.5.

	MASTER BEAM	EPBM	1
C	**********	EPBM	2
C		EPBM	3
C ***	ELSTO-PLASTIC NONLAYERED TIMOSHENKO BEAM PROGRAM	EPBM	4
C		EPBM	5
C	**********	EPBM	6
.	COMMON/UNIM1/NPOIN,NELEM,NBOUN,NLOAD,NPROP,NNODE,IINCS,IITER,	EPBM	7
	KRESL,NCHEK,TOLER,NALGO,NSVAB,NDOFN,NINCS,NEVAB,	EPBM	8
	NITER,NOUTP,FACTO	EPBM	9
	COMMON/UNIM2/PROPS(5,4),COORD(26),LNODS(25,2),IFPRE(52),	EPBM	10
	FIXED(52),TLOAD(25,4),RLOAD(25,4),ELOAD(25,4),	EPBM	11
	MATNO(25),STRES(25,2),PLAST(25),XDISP(52),	EPBM	12

```csv TDISP(26,2),TREAC(26,2),ASTIF(52,52),ASLOD(52), REACT(52),FRESV(1352),PEFIX(52),ESTIF(4,4) CALL DATA CALL INITIAL DO 30 IINCS=1,NINCS CALL INCLOD DO 10 IITER=1,NITER CALL NONAL IF(KRESL.EQ.1) CALL STIFFB CALL ASSEMB IF(KRESL.EQ.1) CALL GREDUC IF(KRESL.EQ.2) CALL RESOLV CALL BAKSUB CALL REFORB CALL CONUND IF(NCHEK.EQ.0) GO TO 20 IF(IITER.EQ.1.AND.NOUTP.EQ.1) CALL RESULT IF(NOUTP.EQ.2) CALL RESULT 10 CONTINUE WRITE(6,900) 900 FORMAT(1H0,5X,'SOLUTION NOT CONVERGED') STOP 20 CALL RESULT 30 CONTINUE STOP END EPBM 13 EPBM 14 EPBM 15 EPBM 16 EPBM 17 EPBM 18 EPBM 19 EPBM 20 EPBM 21 EPBM 22 EPBM 23 EPBM 24 EPBM 25 EPBM 26 EPBM 27 EPBM 28 EPBM 29 EPBM 30 EPBM 31 EPBM 32 EPBM 33 EPBM 34 EPBM 35 EPBM 36 EPBM 37 EPBM 38 . ``` Subroutine STIFFB The purpose of this routine is to evaluate the element stiffness matrices and store them on disc prior to their use in the assembly and equation solving routines. ```csv SUBROUTINE STIFFB STFB 1 C**************************STFB 2 C STFB 3 C *** CALCULATES ELEMENT STIFFNESS MATRICES STFB 4 C STFB 5 C**************************STFB 6 COMMON/UNIM1/NPOIN,NELEM,NBOUN,NLOAD,NPROP,NNODE,IINCS,IITER, STFB 7 KRESL,NCHEK,TOLER,NALGO,NSVAB,NDOFN,NINCS,NEVAB, STFB 8 NITER,NOUTP,FACTO STFB 9 COMMON/UNIM2/PROPS(5,4),COORD(26),LNODS(25,2),IFPRE(52), STFB 10 FIXED(52),TLOAD(25,4),RLOAD(25,4),ELOAD(25,4), STFB 11 MATNO(25),STRES(25,2),PLAST(25),XDISP(52), STFB 12 TDISP(26,2),TREAC(26,2),ASTIF(52,52),ASLOD(52), STFB 13 REACT(52),FRESV(1352),PEFIX(52),ESTIF(4,4) STFB 14 REWIND 1 STFB 15 DO 20 IELEM=1,NELEM STFB 16 LPROP=MATNO(IELEM) STFB 17 EIVAL=PROPS(LPROP,1) STFB 18 SVALU=PROPS(LPROP,2) STFB 19 HARDS=PROPS(LPROP,4) STFB 20 NODE1=LNODS(IELEM,1) STFB 21 NODE2=LNODS(IELEM,2) STFB 22 ELENG=ABS(COORD(NODE2)-COORD(NODE1)) STFB 23 IF(PLAST(IELEM).NE.0.0) EIVAL=EIVAL*(1.0-EIVAL/(EIVAL+HARDS)) STFB 24 VALU1=0.5*SVALU STFB 25 VALU2=SVALU/ELENG STFB 26 VALU3=EIVAL/ELENG STFB 27 VALU4=0.25*SVALU*ELENG STFB 28 ESTIF(1,1)= VALU2 STFB 29 ESTIF(1,2)= VALU1 STFB 30 ```

ESTIF(1,3)=-VALU2	STFB	31
ESTIF(1,4)=VALU1	STFB	32
ESTIF(2,2)=VALU3+VALU4	STFB	33
ESTIF(2,3)=-VALU1	STFB	34
ESTIF(2,4)=-VALU3+VALU4	STFB	35
ESTIF(3,3)=VALU2	STFB	36
ESTIF(3,4)=-VALU1	STFB	37
ESTIF(4,4)=VALU3+VALU4	STFB	38
DO 10 ISTIF=1,4	STFB	39
DO 10 JSTIF=ISTIF,4	STFB	40
10 ESTIF(JSTIF,ISTIF)=ESTIF(ISTIF,JSTIF)	STFB	41
WRITE(1) ESTIF	STFB	42
20 CONTINUE	STFB	43
RETURN	STFB	44
END	STFB	45

STFB 15 Rewind disc ready for writing element stiffnesses. STFB 16–38 For each element evaluate the upper triangular portion of the element stiffness matrix $K^{(e)}$ . Note that if plasticity has taken place the elastic EI is replaced by the elasto-plastic $(EI)_{T}$ . STFB 39-41 Obtain using symmetry the lower triangular portion of $K^{(e)}$ . STFB 42 Write all element stiffness matrices on to disc. Subroutine REFORB This routine evaluates the equivalent nodal forces.

SUBROUTINE REFORB	RFRB	1
C**********	RFRB	2
C	RFRB	3
C *** CALCULATES INTERNAL EQUIVALENT NODAL FORCES	RFRB	4
C	RFRB	5
C**********	RFRB	6
COMMON/UNIM1/NPOIN.NELEM,NBOUN,NLOAD,NPROP,NNODE,IINCS,IITER,	RFRB	7
KRESL,NCHEK,TOLER,NALGO,NSVAB,NDOFN,NINCS,NEVAB,	RFRB	8
NITER,NOUTP,FACTO	RFRB	9
COMMON/UNIM2/PROPS(5,4),COORD(26),LNODS(25,2),IFPRE(52),	RFRB	10
FIXED(52),TLOAD(25,4),RLOAD(25,4),ELOAD(25,4),	RFRB	11
MATNO(25),STRES(25,2),PLAST(25),XDISP(52),	RFRB	12
TDISP(26,2),TREAC(26,2),ASTIF(52,52),ASLOD(52),	RFRB	13
REACT(52),FRESV(1352),PEFIX(52),ESTIF(4,4)	RFRB	14
DO 10 IELEM=1,NELEM	RFRB	15
DO 10 IEVAB=1,NEVAB	RFRB	16
10 ELOAD(IELEM,IEVAB)=0.0	RFRB	17
DO 70 IELEM=1,NELEM	RFRB	18
LPROP=MATNO(IELEM)	RFRB	19
EIVAL=PROPS(LPROP,1)	RFRB	20
SVALU=PROPS(LPROP,2)	RFRB	21
YIELD=PROPS(LPROP,3)	RFRB	22
HARDS=PROPS(LPROP,4)	RFRB	23
NODE1=LNODS(IELEM,1)	RFRB	24
NODE2=LNODS(IELEM,2)	RFRB	25
ELENG=ABS(COORD(NODE2)-COORD(NODE1))	RFRB	26
WNOD1=XDISP(NODE1*NDOFN-1)	RFRB	27
WNOD2=XDISP(NODE2*NDOFN-1)	RFRB	28
THTA1=XDISP(NODE1*NDOFN)	RFRB	29
THTA2=XDISP(NODE2*NDOFN)	RFRB	30
STRAN=(THTA1-THTA2)/ELENG	RFRB	31
STLIN=STRAN*EIVAL	RFRB	32
STCUR=STRES(IELEM,1)+STLIN	RFRB	33
PREYS=YIELD+HARDS*ABS(PLAST(IELEM))	RFRB	34
IF(ABS(STRES(IELEM,1)).GE.PREYS) GO TO 20	RFRB	35

```txt ESCUR=ABS(STCUR)-PREYS RFRB 36 IF(ESCUR.LE.0.0) GO TO 40 RFRB 37 RFACT=ESCUR/ABS(STLIN) RFRB 38 GO TO 30 RFRB 39 20 IF(STRES(IELEM,1).GT.0.0.AND.STLIN.LE.0.0) GO TO 40 RFRB 40 IF(STRES(IELEM,1).LT.0.0.AND.STLIN.GE.0.0) GO TO 40 RFRB 41 RFACT=1.0 RFRB 42 30 REDUC=1.0-RFACT RFRB 43 STRES(IELEM,1)=STRES(IELEM,1)+REDUC*STLIN+ RFRB 44 RFACT*EIVAL*(1.0-EIVAL/(EIVAL+HARDS))*STRAN RFRB 45 PLAST(IELEM)=PLAST(IELEM)+RFACT*STRAN*EIVAL/(EIVAL+HARDS) RFRB 46 GO TO 50 RFRB 47 40 STRES(IELEM,1)=STRES(IELEM,1)+STLIN RFRB 48 50 STRES(IELEM,2)=STRES(IELEM,2)+(SVALU/ELENG)*(WNOD2-WNOD1) RFRB 49 -0.5*SVALU*(THTA1+THTA2) RFRB 50 ELOAD(IELEM,1)=ELOAD(IELEM,1)-STRES(IELEM,2) RFRB 51 ELOAD(IELEM,2)=ELOAD(IELEM,2)+STRES(IELEM,1) RFRB 52 -0.5*ELENG*STRES(IELEM,2) RFRB 53 ELOAD(IELEM,3)=ELOAD(IELEM,3)+STRES(IELEM,2) RFRB 54 ELOAD(IELEM,4)=ELOAD(IELEM,4)-STRES(IELEM,1) RFRB 55 -0.5*ELENG*STRES(IELEM,2) RFRB 56 70 CONTINUE RFRB 57 RETURN RFRB 58 END RFRB 59 ``` RFRB 15-17 Zero space for storing $\pmb{p}$ . RFRB 18-57 For each element evaluate $p^{(e)}$ and assemble into $p$ . # 5.4.6 Examples of nonlayered elasto-plastic Timoshenko beam analysis Two numerical examples are considered. The first example, Example 5.1, involves the yielding of a rectangular simple beam under uniformly distributed load. The beam material has the following properties: $$ E = 2 1 0 \cdot 0 \mathrm{kN} / \mathrm{mm} ^ {2} $$ $$ \nu = 0 \cdot 3 $$ $$ \sigma_ {0} = 0 \cdot 2 5 \mathrm{kN} / \mathrm{mm} ^ {2} $$ $$ H ^ {\prime} = 0 \cdot 0 $$ and the beam proportions are: $$ b = 1 5 0 \mathrm{mm} $$ $$ t = 3 0 0 \mathrm{mm} $$ $$ l = 3 0 0 0 \mathrm{mm} $$ Typical input data is provided in Appendix IV. The problem, finite element idealisation employed and the results are illustrated in Fig. 5.6, which shows that the beam fails as soon as a plastic hinge forms at the centre of the beam. Note that the beam material is assumed to have no strain hardening. The second example considered, Example 5.2, is the clamped I beam shown in Fig. 5.7. The beam has the same material properties as those of Example 5.1. The dimensions and finite element discretisation of the beam are given in Fig. 5.7; the load-displacement relationship at the beam centre is also provided. It is seen that there is an initial loss of stiffness corresponding to the 

line

| Central deflection (mm) | Applied load (KN) | | ----------------------- | ----------------- | | 0 | 0 | | 10 | 1260 | | 10 @ 300 mm | 1260 | | 3000 mm | 1260 | | 300 mm | 1260 | | 150 mm | 1260 | | 300 mm | 1260 |

 Fig. 5.7 Nonlayered elasto-plastic clamped beam.