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.
*/

/*
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(xi) := (XB-XG)*xi + (XA-XG)*(1-xi)$
y_AB(xi) := (YB-YG)*xi + (YA-YG)*(1-xi)$
x_BC(xi) := (XC-XG)*xi + (XB-XG)*(1-xi)$
y_BC(xi) := (YC-YG)*xi + (YB-YG)*(1-xi)$
x_CD(xi) := (XD-XG)*xi + (XC-XG)*(1-xi)$
y_CD(xi) := (YD-YG)*xi + (YC-YG)*(1-xi)$

/*
calcolo aree

Area calculation by integration along the section
*/

area :  integrate(t * l_AB ,xi,0,1)
      + integrate(t * l_BC ,xi,0,1)
      + integrate(t * l_CD ,xi,0,1);

/*
momento statico del primo ordine

First order (or "statical") moments of area, with respect to the Center of Gravity
*/

Mx  :   integrate(y_AB(xi) * t * l_AB ,xi,0,1)
      + integrate(y_BC(xi) * t * l_BC ,xi,0,1)
      + integrate(y_CD(xi) * t * l_CD ,xi,0,1);

My  :   integrate(x_AB(xi) * t * l_AB ,xi,0,1)
      + integrate(x_BC(xi) * t * l_BC ,xi,0,1)
      + integrate(x_CD(xi) * 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(xi)^2 * t * l_AB ,xi,0,1)
      + integrate(y_BC(xi)^2 * t * l_BC ,xi,0,1)
      + integrate(y_CD(xi)^2 * t * l_CD ,xi,0,1)
    );

Jyy : fullratsimp(
        integrate(x_AB(xi)^2 * t * l_AB ,xi,0,1)
      + integrate(x_BC(xi)^2 * t * l_BC ,xi,0,1)
      + integrate(x_CD(xi)^2 * t * l_CD ,xi,0,1)
    );

Jxy :   integrate(x_AB(xi)*y_AB(xi) * t * l_AB ,xi,0,1)
      + integrate(x_BC(xi)*y_BC(xi) * t * l_BC ,xi,0,1)
      + integrate(x_CD(xi)*y_CD(xi) * 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(dummy),y_AB(dummy)) * t * l_AB,
                dummy,0,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(dummy),y_BC(dummy)) * t * l_BC,
                dummy,0,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(dummy),y_CD(dummy)) * t * l_CD,
                dummy,0,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

Shear resultants Sx,Sy are applied for plotting purposes s.t. the nominal shear stress, 
avaluated as the load divided by the area, is unitary.
*/

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;

those vectors are a function of xi (and should be defined accordingly) in the
case of curved segments.
*/

vx_AB : - diff(x_AB(xi),xi)/l_AB;
vy_AB : - diff(y_AB(xi),xi)/l_AB;

vx_BC : - diff(x_BC(xi),xi)/l_BC;
vy_BC : - diff(y_BC(xi),xi)/l_BC;

vx_CD : - diff(x_CD(xi),xi)/l_CD;
vy_CD : - diff(y_CD(xi),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(xi)) + vy_AB * (+x_AB(xi))) * tau_AB(xi) * t * l_AB
        + (vx_BC * (-y_BC(xi)) + vy_BC * (+x_BC(xi))) * tau_BC(xi) * t * l_BC
        + (vx_CD * (-y_CD(xi)) + vy_CD * (+x_CD(xi))) * 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 based on comparison with known 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 the 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

a numerical value may then be derived for the specific dimensions.
*/

ev( [chi_x,chi_y,chi_xy],dim,numer);

/*
noto che non sono funzione dello spessore

it may be observed that those coefficients do not depend on the thickness;
the only remaining geometrical ratio they may dependo on is hence h/b.
*/

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]);