Page 4 - Demo
P. 4
Planar Truss Example for Comrel Add-on RCP Consult, 2021-2026 Page 4 fc[14 + 2*n] := -p[n]: # negative vertical forces [N]end do:for n from upperbound(cc) by -1 to 1 do fc := Remove(fc, cc[n]):end do:## Global stiffness matrixK := Matrix(ndof, ndof): # stiffness matrixT := Array(1..nfe): # transformation matrixfor n to nfe do # sin and cos of the inclination c := cos(theta[n]): s := sin(theta[n]): # coordinate transformation matrix for the bar n T[n] := Matrix([[c, s, 0, 0], [-s, c, 0, 0], [0, 0, c, s], [0, 0, -s, c]]): # local stiffness matrix for the bar n Kloc := k[n]*Matrix([[1, 0, -1, 0], [0, 0, 0, 0], [-1, 0, 1, 0], [0, 0, 0, 0]]): # global stiffness matrix assembly idx:=[entries(LM[n,1..],'nolist')]: K[idx,idx]:=K[idx,idx]+Array(evalm(transpose(T[n]).Kloc.T[n])):end do:## Solve the system of equationsec := [entries(cc, 'nolist')]:ed := [entries(dd, 'nolist')]:ac := Vector([0, 0, 0]):Kcc := SubMatrix(K, ec, ec):Kcd := SubMatrix(K, ec, ed):Kdc := SubMatrix(K, ed, ec):Kdd := SubMatrix(K, ed, ed):ad := LinearSolve(Kdd, fc - Multiply(Kdc, ac)):qd := Multiply(Kcc, ac) + Multiply(Kcd, ad):## Final displacementsa := Vector(ndof):q := Vector(ndof):a[ec] := ac:q[ec] := qd:a[ed] := ad:N := Vector(nfe):## Axial loads in barsfor n to nfe do N[n] := Multiply(Multiply(Multiply(k[n], Array([-1, 0, 1, 0])), T[n]), a(LM[n])):end do:## Visualisationif StrurelPlot then # begin of plot block if StrurelMode = 0 or StrurelMode = 2 then global StrurelPlotCount: if StrurelMode = 0 then tit:=\ elif StrurelMode = 2 then tit:=\# 2 - final stochastic solution end if: # Get a screen resolution getmetric:=define_external(\ ## Plot Function to be used in Maplet TrussPlot := proc(np) # X and Y coordinates in m XY := Matrix([[0, 0],[4, 0],[8, 0],[12, 0],[16, 0],[20, 0],[24, 0],[22, 2],[18, 2],[14, 2],[10, 2],[6, 2],[2, 2]]): if np = 1 then PlotT:=textplot([12.5,6.,\ BarL := Array(1..nfe): # rods as lines BarR := Array(1..nfe): # rod numbers NodN := Array(1..nnp): # nodal numbers # Plot initial state of truss for i to nfe do x1:=XY[IEN[i,1],1]: y1:=XY[IEN[i,1],2]: x2:=XY[IEN[i,2],1]: y2:=XY[IEN[i,2],2]: BarL[i]:=line([x1,y1],[x2,y2],color=blue): # blue color for initial state BarR[i]:=textplot([(x1+x2)/2,(y1+y2)/2,sprintf(\ end do: for i to nnp do NodN[i]:=textplot([XY[i,1],XY[i,2],sprintf(\# nodal numbers end do: BarD := Array(1..nfe): # deformed rods as lines # Plot deformed state of truss - scale factor 2 for i to nfe do n1:=IEN[i,1]: n2:=IEN[i,2]: x1:=XY[n1,1]+a[n1*2-1]*2: y1:=XY[n1,2]+a[n1*2]*2: x2:=XY[n2,1]+a[n2*2-1]*2: y2:=XY[n2,2]+a[n2*2]*2:
