kill(all); /* dichiaro le quote dimensionali essere ristrette in segno, in modo da avere |x| = x Dimensions are declared positive, in order to simplify abs(x) with x */ assume(b>0,t>0,h>0,c>0); /* considero per semplicità uno spessore costante t For the sake of simplicity, wall thickness is assumed constant. */ /* area */ area : (b+b+h)*t; /* two coordinate systems are employed: OXY is centered at the web midpoint Gxy is centered at the (initially unknown) section center of gravity */ /* coordinate dei punti nodali del profilo profile nodal point coordinates according to the OXY reference system */ XA : b \$ YA: h/2\$ XB : 0 \$ YB: h/2\$ XC : 0 \$ YC:-h/2\$ XD : b \$ YD:-h/2\$ /* derivo lunghezze dei tratti profile segment lengths */ l_AB : sqrt((XA-XB)^2 + (YA-YB)^2)\$ l_BC : sqrt((XB-XC)^2 + (YB-YC)^2)\$ l_CD : sqrt((XC-XD)^2 + (YC-YD)^2)\$ /* coordinate, in funzione di un'ascissa curvilinea adimensionale xi in [0,1] (x,y) coordinates for points along the profile segments are parametrically defined with respect to a normalized arlength coordinate xi that spans from 0 to 1 while moving along the segment from the first to the second endpoint. The Gxy centroidal coordinate system is here employed, with unknown XG,YG. */ x_AB : (XB-XG)*xi + (XA-XG)*(1-xi)\$ y_AB : (YB-YG)*xi + (YA-YG)*(1-xi)\$ x_BC : (XC-XG)*xi + (XB-XG)*(1-xi)\$ y_BC : (YC-YG)*xi + (YB-YG)*(1-xi)\$ x_CD : (XD-XG)*xi + (XC-XG)*(1-xi)\$ y_CD : (YD-YG)*xi + (YC-YG)*(1-xi)\$ /* momento statico del primo ordine First order (or "statical") moments of area, with respect to the Center of Gravity */ Mx : integrate(y_AB * t * l_AB ,xi,0,1) + integrate(y_BC * t * l_BC ,xi,0,1) + integrate(y_CD * t * l_CD ,xi,0,1); My : integrate(x_AB * t * l_AB ,xi,0,1) + integrate(x_BC * t * l_BC ,xi,0,1) + integrate(x_CD * t * l_CD ,xi,0,1); /* derivo baricentro imponendo momento statico nullo se calcolato su sistema di riferimento baricentrico the first order moment of area is known to be null if the reference system is centroidal; such condition is employed to evaluate the center of gravity coordinates with respect to OXY */ linsolve([Mx=0,My=0],[XG,YG]), globalsolve=true; /* ridefinisco coordinate noto il baricentro Multiple assignments that redefine the (x,y) coordinates along the segments according to the just evaluated center of gravity position. */ [x_AB , y_AB ,x_BC ,y_BC ,x_CD ,y_CD ] : ev([x_AB , y_AB ,x_BC ,y_BC ,x_CD ,y_CD ],infeval); /* momenti del secondo ordine second order moments of area, with respect to Gxy */ Jxx : fullratsimp( integrate(y_AB^2 * t * l_AB ,xi,0,1) + integrate(y_BC^2 * t * l_BC ,xi,0,1) + integrate(y_CD^2 * t * l_CD ,xi,0,1) ); Jyy : fullratsimp( integrate(x_AB^2 * t * l_AB ,xi,0,1) + integrate(x_BC^2 * t * l_BC ,xi,0,1) + integrate(x_CD^2 * t * l_CD ,xi,0,1) ); Jxy : integrate(x_AB*y_AB * t * l_AB ,xi,0,1) + integrate(x_BC*y_BC * t * l_BC ,xi,0,1) + integrate(x_CD*y_CD * t * l_CD ,xi,0,1); /* due to the symmetric nature of the cross section, the product (not "mixed") moment of area is zero. */ /* coefficienti alpha, beta we now define the alpha,beta coefficients of the equation d sigma_z --------- = alpha(x,y,...) Sy - beta(x,y,...) *Sx dz Material is assumed homogeneous and isotropic along the section. */ define( alpha(x,y) , fullratsimp ( (-Jxy * x + Jyy*y)/(Jxx*Jyy-Jxy^2) )); define( beta(x,y) , fullratsimp ( (-Jxx * x + Jxy*y)/(Jxx*Jyy-Jxy^2) )); /* inequilibrio dsigmaz_dz for each point along the section, the axial stress component variation with increasing z is defined. Such variation causes a disequilibrium that induces the compensating shear actions. */ define( dsigmaz_dz(x,y) , fullratsimp( alpha(x,y) * Sy - beta(x,y) * Sx ) ); /* tensione tagliante all'estremo A shear stress at the A end, which is a free surface in the case of an open section. In the case of a closed thin wall section, such value is generally nonzero, and it is to be treated as an unknown parameter. */ tau_A : 0; /* tratto AB : tensione tagliante in funzione dell'ascissa curvilinea adimensionalizzata. AB segment : The tau_sz shear stress component is evaluated as a function of the normalized arc length coordinate xi, based on the axial equilibrium of the portion of section spanning from A (s=0) to the current point s=xi*l_AB along the segment */ define( tau_AB(xi), fullratsimp( ( tau_A *t + integrate( dsigmaz_dz(x_AB,y_AB) * t * l_AB, xi ) )/t ) ); /* tratto BC : tensione tagliante in funzione dell'ascissa curvilinea adimensionalizzata BC segment: the shear stress component is evaluated as for segment AB; at the B segment end, the shear stress transmitted from the previous segment is applied. The same normalized arc length coordinate is considered, whereas the absolute arclength coordinate (if needed) may be defined as s = l_AB + xi*l_BC, ds = l_BC * dxi */ define( tau_BC(xi), fullratsimp( ( tau_AB(1) *t + integrate( dsigmaz_dz(x_BC,y_BC) * t * l_BC, xi ) )/t ) ); /* tratto CD : tensione tagliante in funzione dell'ascissa curvilinea adimensionalizzata BC segment: the shear stress component is evaluated as for segments AB and CD. here, s = l_AB + l_BC + xi*l_CD, ds = l_CD * dxi */ define( tau_CD(xi), fullratsimp( ( tau_BC(1) *t + integrate( dsigmaz_dz(x_CD,y_CD) * t * l_CD, xi ) )/t ) ); /* controllo: è tau_D nulla? check : is tau_sz equal to zero at the D end?? */ tau_D : fullratsimp(tau_CD(1)); /* wow. ok. */ /* altro controllo: plot dei grafici per uno specifico dimensionamento further check : let's plot the tau_zs variation along the cross section segments. sample dimensions are required for extracting numerical results to be plotted. */ dim : [h=80,b=40,t=6]; /* applico tagli in direzione x,y t.c. la tensione nominale (taglio/area) sia unitaria */ wxplot2d( ev([tau_AB(xi),tau_BC(xi),tau_CD(xi)],dim,Sy=area*1,Sx=0,infeval), [xi,0,1], [legend,"AB","BC","CD"] ); wxplot2d( ev([tau_AB(xi),tau_BC(xi),tau_CD(xi)],dim,Sy=0,Sx=area*1,infeval), [xi,0,1], [legend,"AB","BC","CD"] ); /* ok, pare credibile. ok, it seems reasonable. */ /* versore associato alla tau sul tratto in order to evaluate resultants, a unit vector is to be defined for each segment that orients the tau_sz actions along the plane. vx = - diff ( x, s) = -diff (x, xi) / l; */ vx_AB : - diff(x_AB,xi)/l_AB; vy_AB : - diff(y_AB,xi)/l_AB; vx_BC : - diff(x_BC,xi)/l_BC; vy_BC : - diff(y_BC,xi)/l_BC; vx_CD : - diff(x_CD,xi)/l_CD; vy_CD : - diff(y_CD,xi)/l_CD; /* risultante delle azioni in direzione x x component of the resultant force of the tau_sx stress distribution. */ Fx : integrate( vx_AB * tau_AB(xi) * t * l_AB + vx_BC * tau_BC(xi) * t * l_BC + vx_CD * tau_CD(xi) * t * l_CD, xi,0,1 ); /* risultante delle azioni in direzione y y component of the resultant force of the tau_sx stress distribution. */ Fy : integrate( vy_AB * tau_AB(xi) * t * l_AB + vy_BC * tau_BC(xi) * t * l_BC + vy_CD * tau_CD(xi) * t * l_CD, xi,0,1 )\$ Fy : fullratsimp(Fy); /* componente z del momento risultante rispetto al baricentro axial z component of the resultant moment of the tau_sx stress distribution. */ MGz : integrate( (vx_AB * (-y_AB) + vy_AB * (+x_AB)) * tau_AB(xi) * t * l_AB + (vx_BC * (-y_BC) + vy_BC * (+x_BC)) * tau_BC(xi) * t * l_BC + (vx_CD * (-y_CD) + vy_CD * (+x_CD)) * tau_CD(xi) * t * l_CD, xi,0,1 )\$ MGz : fullratsimp(MGz); /* posizione del centro di taglio in coordinate (e,f) rispetto al baricentro The evaluation follows of the shear center coordinates (e,f), based on an equal resultant moment condition with respect to the CoG. */ eq : MGz = Sx * (-f) + Sy * (+e); /* vale per Sx=0 e Sy=0, e deve mantenere tale valore al variare arbitrario di Sx e Sy such equation holds for null shear components, we now impose that such null residual condition is also constant for arbitrary varying Sx, Sy */ linsolve( [diff(eq,Sx), diff(eq,Sy)],[e,f] ), globalsolve=true; /* confronto con risultati terzi, tipo https://theconstructor.org/structural-engg/analysis/shear-centre-with-examples/3677/ controllo la differenza check on independent results from literature */ fullratsimp( ( -(XG + b^2*h^2*t/4/Jxx) ) - ( e ) ); /* ok */ /* valuto ora l'energia potenziale elastica associata alla distribuzione di taglio, per unità di lunghezza di trave the shear contribution to the elastic strain energy per unit beam length is evaluated in the following */ dU_dz : integrate( tau_AB(xi)^2/2/G * t* l_AB + tau_BC(xi)^2/2/G * t* l_BC + tau_CD(xi)^2/2/G * t* l_CD, xi,0,1 ); /* confronto con la sua controparte nominale corretta dai fattori chi a nominal counterpart is defined, that contains the three corrective chi factors */ eq : dU_dz = (chi_x * Sx^2 + chi_y * Sy^2 + chi_xy *Sx*Sy)/2/G/area; /* non posso derivare in Sx^2,Sy^2,Sx*Sy, aggiro il problema con un cambio di variabile since it is not allowed to differentiate any expression with respect to non-atomic terms as Sx^2,Sy^2,Sx*Sy, a workaround is applied in the following consisting in a change of variable. */ eq : ratsubst(Sxq , Sx^2 , eq)\$ eq : ratsubst(Syq , Sy^2 , eq)\$ eq : ratsubst(SxSy , Sx*Sy, eq); /* posso ora derivare I'm now allowed to differentiate with respect to such just defined symbols */ linsolve( [diff(eq,Sxq),diff(eq,Syq),diff(eq,SxSy)], [chi_x,chi_y,chi_xy] ), globalsolve=true; /* per farmi un'idea, valuto nel dimensionamento specifico ... */ ev( [chi_x,chi_y,chi_xy],dim,numer); /* noto che non sono funzione dello spessore ... */ ev( [chi_x,chi_y,chi_xy],b=1,h=l); wxplot2d(%,[l,0.1,10],[logx],[legend,"chi_x","chi_y","chi_{xy}"],[y,0-1,5]);