Page 3 - Demo
P. 3
Planar Truss Examplefor ComrelAdd-onRCPConsult, 2021-2025Page 3%% Code based on Max Ehre work https://github.com/maxehre/polynomial_surrogates%% Used with his kind permission as version from 2017 that adapted to Strurel.functionuout=PlanarTruss(XP)%{============================================================================Planar Truss Example for Matlab============================================================================Inputs:X = [A1, A2, E1, E2, P1, P2, P3, P4, P5, P6, u_y]E1,E2: Youngs modulus of horizontal and inclined bars resp. -log-normally distributedA1,A2: cross section of horizontal and inclined bars resp. lognormally distributedP1-P6: negative vertical forces attacking in nodes 8-13. -Gumbel distributedu_y: maximal allowed vertical deflection in node 4. constantOutput:uout: vertical truss deflection at bottom center (N04)============================================================================Truss description taken from Lee, S.H., Kwak, B.M. (2006).Response surface augmented moment method for efficient reliability analysis.Structural Safety 28(3), 261 -272.https://www.sciencedirect.com/science/article/abs/pii/S0167473005000421Based on Diego Andr%u00e9s Alvarez Mar%u00edn work https://github.com/diegoandresalvarez/elementosfinitos/tree/master/codigo/repaso_matricial/cercha_2dN: NodesR: Rods============================================================================N13__R4___N12__R8___N11__R12__N10__R16__N09__R20__N08/\\/\\/\\/\\/\\/\\/ \\/ \\/ \\/ \\/ \\/ \\R1 R3 R5 R7 R9 R11 R13 R15 R17 R19 R21 R23/ \\/ \\/ \\/ \\/ \\/ \\/___R2___\\/___R6___\\/__R10___\\/__R14___\\/__R18___\\/__R22___\\N01 N02 N03 N04 N05 N06 N07||\\/uout => u_y @ N04============================================================================%}%%---------------------------------------------------------------------------%% Planar truss model%%---------------------------------------------------------------------------%% Global flag to manage plot facilities of engineStrurelPlot=true;%% Vector inputs order from Stochastic ModelA1=XP(1);% cross section of horizontal barsA2=XP(2);% cross section of inclined barsE1=XP(3);% Youngs modulus of horizontal barsE2=XP(4);% Youngs modulus of inclined barsP=XP(510);% negative vertical forces attacking in nodes 813. (1x6 vector)%% Initial input variablesnfe =23;% number of elements (bars)ned =2;% number of dof per nodennp =13;% number of nodal pointsndof =ned*nnp;% number of degrees of freedom (dof)%% Element, nodes and dofs association%% IEN: connectivity matrix, nfe x nodes bar 1 has nodes 1 and 13, bar 2 has nodes 1 and 2 ...IEN =[113;12;132;1312;212;23;123;1211;311;34;114;1110;410;45;105;109;59;56;96;98;68;67;87];%% LM: localization matrix, nfe x dof LM =cell(nfe,1);fore =1nfeLM{e}=[IEN(e,1)*ned-1,IEN(e,1)*ned,IEN(e,2)*ned-1,IEN(e,2)*ned];% element 1 has dof 1-2 & 25-26 ...end%% Deterministic rod propertiesGrid =4;% size of grid side in [m]ng =12;% count of grid cellsang =atan2(Grid,Grid)*180/pi;% inclination angle of the truss [deg]theta =repmat([ang 0-ang 0],1,ng/2);theta(end)=[];% inclination angle [deg]leng =repmat(4*[1/sqrt(2)1],1,ng);leng(end)=[];% bar length [m]%% Stochastic rod propertiesArea =repmat([A2A1],1,ng);Area(end)=[];% bar cross sectional area [m2]E =repmat([E2E1],1,ng);E(end)=[];% Young's modulus [N/m^2]%% Material propertiesk =E.*Area./leng;% stiffness of each bar%% Stiffness matrix assemblyK =zeros(ndof);% stiffness matrixT =cell(nfe,1);% transformation matrixfore =1nfe % sin and cos of the inclinationc =cosd(theta(e));s =sind(theta(e));% coordinate transformation matrix for the bar eT{e}=[c s 00;s c 00;00c s;00s c ];% local stiffness matrix for the bar e Kloc =k(e)*[1010;0000;1010;0000];% global stiffness matrix assemblyK(LM{e},LM{e})=K(LM{e},LM{e})+T{e}'*Kloc*T{e};end