Page 9 - Demo

Planar Truss Examplefor ComrelAdd-onRCPConsult, 2023-2025Page 9!! 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)doi=1,nbdkeT(i)=0.d0doj=1,nbdkeT(i)=keT(i)+TT(i,j)*ken(j)enddoenddodoi=1,nbdde(i)=ac(LM(i,n))enddoNN(n)=0.d0doi=1,nbdNN(n)=NN(n)+keT(i)*de(i)enddoenddo!!* Visualizationif(StrurelPlot)thenif(StrurelMode.eq.0.or.StrurelMode.eq.2)then! Local variables: width, heigth, ouput window number, etc.wo=-2.d0;we=26.d0;ho=-2.d0;he=10.d0