Page 4 - Demo

Planar Truss Examplefor ComrelAdd-onRCPConsult, 2023-2025Page 4E2=v_XP_(4)! Youngs modulus of inclined bars!! LM: localization matrix, dof x nfe don=1,nfeLM(1,n)=IEN(1,n)*ned1;LM(2,n)=IEN(1,n)*ned LM(3,n)=IEN(2,n)*ned-1;LM(4,n)=IEN(2,n)*nedenddo!! Deterministic rod properties! inclination angle of the truss [deg]ang=atan2(200.d0,200.d0)*180.d0/PIdon=1,6m=(n-1)*4! inclination angles [deg]theta(m+1)=ang;theta(m+2)=0.d0;theta(m+3)=-ang;theta(m+4)=0.d0enddol1=4.d0/dsqrt(2.d0);l2=4.d0don=1,12m=(n-1)*2leng(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]enddo!! Material propertiesdon=1,nfeke(n)=E(n)*Area(n)/leng(n)enddo!! Global stiffness matrixcallMATNUL(K,ndof,ndof)! cleandon=1,nferad=theta(n)*PI/180.d0! sin and cos of the inclinationc=dcos(rad)s=dsin(rad)! coordinate transformation matrix for the bar ncallMATNUL(T,nbd,nbd)T(1,1)=c;T(1,2)=sT(2,1)=s;T(2,2)=cT(3,3)=c;T(3,4)=sT(4,3)=s;T(4,4)=c! local stiffness matrix Kloc for the bar n callMATNUL(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 TTcallMATRAN(T,TT,nbd,nbd)! get product T'*Kloc*T in KGcallMATMLT(TT,KL,KT,nbd,nbd,nbd)callMATMLT(KT,T,KG,nbd,nbd,nbd)! collect coefficients in a global stiffness matrix Kdoi=1,nbddoj=1,nbdK(LM(i,n),LM(j,n))=K(LM(i,n),LM(j,n))+KG(i,j)enddoenddoenddo!! Boundary conditions! Applied loadscallMATNUL(fc,ndof,1)doi=8,nnpfc(i*ned)=-P(i-7)! negative vertical forces [N]enddo! Supports and free nodes! Remove supports from force vectordoi=1,ncfc(bc(i))=0.d0enddo! Remove supports from K matrix all coefficents on row and col to 0. and diagonal to 1.doi=1,ncdoj=1,ndofK(bc(i),j)=0.d0! clean a row membersenddodoj=1,ndofK(j,bc(i))=0.d0! clean a col membersenddoK(bc(i),bc(i))=1.d0! diagonal memberenddo!! Solve the system of equations K*ac=fccallLINSOL(K,ac,fc,ndof)!! Axial loads in barsdon=1,nfeken(1)=ke(n);ken(2)=0.d0;ken(3)=ke(n);ken(4)=0.d0rad=theta(n)*PI/180.d0! sin and cos of the inclinationc=dcos(rad)s=dsin(rad)! coordinate transformation matrix for the bar ncallMATNUL(T,nbd,nbd)T(1,1)=c;T(1,2)=sT(2,1)=s;T(2,2)=cT(3,3)=c;T(3,4)=sT(4,3)=s;T(4,4)=c! transpose T in TTcallMATRAN(T,TT,nbd,nbd)! get stiffness of the bardoi=1,nbdkeT(i)=0.d0doj=1,nbdkeT(i)=keT(i)+TT(i,j)*ken(j)enddoenddo! get deformations in bar nodesdoi=1,nbdde(i)=ac(LM(i,n))enddo! compute a normal force in the barNN(n)=0.d0