Page 4 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2023-2025 Page 4 LM(3,n)=IEN(2,n)*ned-1; LM(4,n)=IEN(2,n)*ned end do!! Deterministic rod properties! inclination angle of the truss [deg] ang=atan2(200.d0,200.d0)*180.d0/PI do n=1,6 m=(n-1)*4! inclination angles [deg] theta(m+1)=ang; theta(m+2)=0.d0; theta(m+3)=-ang; theta(m+4)=0.d0 end do l1=4.d0/dsqrt(2.d0); l2=4.d0 do n=1,12 m=(n-1)*2 leng(m+1)=l1; leng(m+2)=l2 ! bar length [m] Area(m+1)=A2; Area(m+2)=A1 ! bar cross sectional area [m2] E(m+1)=E2; E(m+2)=E1 ! Young's modulus [N/m^2] end do!! Material properties do n=1,nfe ke(n)=E(n)*Area(n)/leng(n) end do!! Global stiffness matrix call MATNUL(K,ndof,ndof) ! clean do n=1,nfe rad=theta(n)*PI/180.d0 ! sin and cos of the inclination c=dcos(rad) s=dsin(rad)! coordinate transformation matrix for the bar n call MATNUL(T,nbd,nbd) T(1,1)= c; T(1,2)=s T(2,1)=-s; T(2,2)=c T(3,3)= c; T(3,4)=s T(4,3)=-s; T(4,4)=c! local stiffness matrix Kloc for the bar n call MATNUL(KL,nbd,nbd) KL(1,1)= ke(n); KL(1,3)=-ke(n) KL(3,1)=-ke(n); KL(3,3)= ke(n)! transpose T in TT call MATRAN(T,TT,nbd,nbd)! get product T'*Kloc*T in KG call MATMLT(TT,KL,KT,nbd,nbd,nbd) call MATMLT(KT, T,KG,nbd,nbd,nbd)! collect coefficients in a global stiffness matrix K do i=1,nbd do j=1,nbd K(LM(i,n),LM(j,n))=K(LM(i,n),LM(j,n))+KG(i,j) end do end do end do!! Boundary conditions! Applied loads call MATNUL(fc,ndof,1) do i=8,nnp fc(i*ned)=-P(i-7) ! negative vertical forces [N] end do! Supports and free nodes! Remove supports from force vector do i=1,nc fc(bc(i))=0.d0 end do! Remove supports from K matrix - all coefficents on row and col to 0. and diagonal to 1. do i=1,nc do j=1,ndof K(bc(i),j)=0.d0 ! clean a row members end do do j=1,ndof K(j,bc(i))=0.d0 ! clean a col members end do K(bc(i),bc(i))=1.d0 ! diagonal member end do!! Solve the system of equations K*ac=fc call LINSOL(K,ac,fc,ndof)!! Axial loads in bars do n=1,nfe ken(1)=-ke(n); ken(2)=0.d0; ken(3)=ke(n); ken(4)=0.d0 rad=theta(n)*PI/180.d0 ! sin and cos of the inclination c=dcos(rad) s=dsin(rad)! coordinate transformation matrix for the bar n call MATNUL(T,nbd,nbd) T(1,1)= c; T(1,2)=s T(2,1)=-s; T(2,2)=c T(3,3)= c; T(3,4)=s T(4,3)=-s; T(4,4)=c! transpose T in TT call MATRAN(T,TT,nbd,nbd)! get stiffness of the bar do i=1,nbd keT(i)=0.d0 do j=1,nbd keT(i)=keT(i)+TT(i,j)*ken(j) end do end do! get deformations in bar nodes do i=1,nbd de(i)=ac(LM(i,n)) end do! compute a normal force in the bar NN(n)=0.d0 do i=1,nbd NN(n)=NN(n)+keT(i)*de(i) end do end do