kill(all);

1 = x
2 = y

relative representation of b and h

h is kept dimensional
b is defined through the upsilon = b/h ratio
t is defined through the epsilon = t/h ratio

b : h * upsilon;

t : h * epsilon;

Dimensions are declared positive, 
in order to simplify abs(x) with x

assume(
 upsilon>0,
 epsilon>0, 
 h>0
);

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

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$

(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)$

profile segment length functions

define(
    dl_dxi_AB(xi),
    sqrt( diff(x_AB(xi),xi)^2 + diff(y_AB(xi),xi)^2   )
);
define(
    dl_dxi_BC(xi),
    sqrt( diff(x_BC(xi),xi)^2 + diff(y_BC(xi),xi)^2   )
);
define(
    dl_dxi_CD(xi),
    sqrt( diff(x_CD(xi),xi)^2 + diff(y_CD(xi),xi)^2   )
);

Area calculation by integration along the section
(not robust for curved segments - but why?)

area :  integrate(t * dl_dxi_AB(xi) ,xi,0,1)
      + integrate(t * dl_dxi_BC(xi) ,xi,0,1)
      + integrate(t * dl_dxi_CD(xi) ,xi,0,1);

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

Mx  :   integrate(y_AB(xi) * t * dl_dxi_AB(xi) ,xi,0,1)
      + integrate(y_BC(xi) * t * dl_dxi_BC(xi) ,xi,0,1)
      + integrate(y_CD(xi) * t * dl_dxi_CD(xi) ,xi,0,1);

My  :   integrate(x_AB(xi) * t * dl_dxi_AB(xi) ,xi,0,1)
      + integrate(x_BC(xi) * t * dl_dxi_BC(xi) ,xi,0,1)
      + integrate(x_CD(xi) * t * dl_dxi_CD(xi) ,xi,0,1);

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;

second order moments of area, with respect to Gxy

Jxx : fullratsimp(
        integrate(y_AB(xi)^2 * t * dl_dxi_AB(xi) ,xi,0,1)
      + integrate(y_BC(xi)^2 * t * dl_dxi_BC(xi) ,xi,0,1)
      + integrate(y_CD(xi)^2 * t * dl_dxi_CD(xi) ,xi,0,1)
    );

Jyy : fullratsimp(
        integrate(x_AB(xi)^2 * t * dl_dxi_AB(xi) ,xi,0,1)
      + integrate(x_BC(xi)^2 * t * dl_dxi_BC(xi) ,xi,0,1)
      + integrate(x_CD(xi)^2 * t * dl_dxi_CD(xi) ,xi,0,1)
    );

Jxy :   integrate(x_AB(xi)*y_AB(xi) * t * dl_dxi_AB(xi) ,xi,0,1)
      + integrate(x_BC(xi)*y_BC(xi) * t * dl_dxi_BC(xi) ,xi,0,1)
      + integrate(x_CD(xi)*y_CD(xi) * t * dl_dxi_CD(xi) ,xi,0,1);

due to the symmetric nature of the cross section, the product (not "mixed") moment
of area is zero.

we now define the alpha,beta coefficients of the equation

d sigma_z
---------  = alpha(x,y,...) Qy - beta(x,y,...) *Qx
    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) ));

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) * Qy - beta(x,y) * Qx ) 
);

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;

AB segment : The tau_iz 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 (xi=0) to the current point along the segment

define(
    tau_AB(xi),
    fullratsimp(
         ( 
              tau_A *t 
            + integrate(
                dsigmaz_dz(x_AB(dummy),y_AB(dummy)) * t * dl_dxi_AB(dummy),
                dummy,0,xi   
              )   
     )/t
    )
);

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. 

define(
    tau_BC(xi),
    fullratsimp(
         ( 
              tau_AB(1) *t 
            + integrate(
                dsigmaz_dz(x_BC(dummy),y_BC(dummy)) * t * dl_dxi_BC(dummy),
                dummy,0,xi   
              )   
  )/t
    )
);

BC segment:  the shear stress component is evaluated as for segments AB and CD.

define(
    tau_CD(xi),
    fullratsimp(
         ( 
              tau_BC(1) *t 
            + integrate(
             dsigmaz_dz(x_CD(dummy),y_CD(dummy)) * t * dl_dxi_CD(dummy),
                dummy,0,xi   
              )   
          )/t
    )
);

check : is tau_iz equal to zero at the D end??

tau_D : fullratsimp(tau_CD(1));

wow. ok.

further check : let's plot the tau_iz variation along the cross 
section segments.
sample dimensions are required for extracting numerical 
results to be plotted.

dim : [h=80,upsilon=40/80,epsilon=6/80];

Shear resultants Qx,Qy 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, Qy=area*1, Qx=0, infeval
    ),
    [xi,0,1],
    [legend,"AB","BC","CD"]
    );

I evaluate the peak value for the actual shear stress vs. 
nominal stress ratio along the BC segment, and along the section.

here, tmp is the evaluation of the stress  

tmp : ev(
            tau_BC(xi), Qy=area*1,  Qx=0 ,infeval
        )$
linsolve(
    diff(
        tmp,
        xi
    )=0,
    xi
)$
fullratsimp(ev(tmp , %));
%,dim,numer;

wxplot2d(
    ev(
        [tau_AB(xi),tau_BC(xi),tau_CD(xi)],
        dim,  Qy=0,  Qx=area*1 ,infeval
    ),
    [xi,0,1],
    [legend,"AB","BC","CD"]
    );

I evaluate the peak value for the actual shear stress vs. 
nominal stress
ratio along the AB segment, and along the section. 

tmp : ev(
            tau_AB(xi), Qx=area*1,  Qy=0 ,infeval
        )$
linsolve(
    diff(
        tmp,
        xi
    )=0,
    xi
)$
fullratsimp(ev(tmp , %));
%,dim,numer;

in order to evaluate resultants, a unit vector is to be defined 
for each
segment that orients the tau_iz actions along the plane.

vx = - diff ( x, s) = -diff (x, xi) / diff (l, xi);
vy = - diff ( y, s) = -diff (y, xi) / diff (l, xi);

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

the negative sign is due to the opposite "inward" and "growing s" directions; 
tau_zi are oriented inward, whereas the

[ diff ( x, s) , diff ( y, s) ]

vector is oriented towards growing s

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

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

vx_CD : - diff(x_CD(xi),xi)/dl_dxi_CD(xi);
vy_CD : - diff(y_CD(xi),xi)/dl_dxi_CD(xi);

x component of the resultant force of the tau_iz stress distribution. 

Fx : integrate(
          vx_AB * tau_AB(xi) * t * dl_dxi_AB(xi)
        + vx_BC * tau_BC(xi) * t * dl_dxi_BC(xi)
        + vx_CD * tau_CD(xi) * t * dl_dxi_CD(xi),
        xi,0,1
    );

y component of the resultant force of the tau_iz stress distribution. 

Fy : integrate(
          vy_AB * tau_AB(xi) * t * dl_dxi_AB(xi)
        + vy_BC * tau_BC(xi) * t * dl_dxi_BC(xi)
        + vy_CD * tau_CD(xi) * t * dl_dxi_CD(xi),
        xi,0,1
    )$
Fy : fullratsimp(Fy);

axial z component of the resultant moment of the tau_iz stress 
distribution.

MGz : integrate(
          (vx_AB * (-y_AB(xi)) + vy_AB * (+x_AB(xi))) * tau_AB(xi) * t * dl_dxi_AB(xi)
        + (vx_BC * (-y_BC(xi)) + vy_BC * (+x_BC(xi))) * tau_BC(xi) * t * dl_dxi_BC(xi)
        + (vx_CD * (-y_CD(xi)) + vy_CD * (+x_CD(xi))) * tau_CD(xi) * t * dl_dxi_CD(xi),
        xi,0,1
    )$
MGz : fullratsimp(MGz);

The evaluation follows of the shear center coordinates (e,f), based on 
an equal resultant moment condition with respect to the CoG.

eq : MGz = Qx * (-f) + Qy * (+e);

eq_res : rhs(eq) - lhs(eq);

such equation holds for null shear components, 

we now impose that such null residual
condition is also constant for arbitrary varying Qx, Qy

ev( eq_res, Qx=0, Qy=0);

linsolve( 
    [
        diff(eq_res,Qx)=0, 
        diff(eq_res,Qy)=0
]
,[e,f] ), globalsolve=true;

check based on comparison with known results from literature, e.g.
https://theconstructor.org/structural-engg/analysis/shear-centre-with-examples/3677/

fullratsimp(
      ( -(XG + b^2*h^2*t/4/Jxx) ) 
    - ( e ) 
);

ok

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* dl_dxi_AB(xi)   
      + tau_BC(xi)^2/2/G * t* dl_dxi_BC(xi)   
      + tau_CD(xi)^2/2/G * t* dl_dxi_CD(xi),
      xi,0,1  
    );

a nominal counterpart is defined, 
that contains the three corrective chi factors

dU_dz_nom : (chi_x * Qx^2 + chi_y * Qy^2 + chi_xy *Qx*Qy)/2/G/area;

eq : dU_dz = dU_dz_nom;

eq_res : lhs(eq) - rhs(eq);

such residual - which is null for null Qx,Qy - has to remain invariant 
with varying Qx, Qy.
However, it is simpler to set it as an invariant in Qx^2, Qy^2, Qx*Qy

i.e.

diff(eq_res, Qx^2 )=0
diff(eq_res, Qy^2 )=0
diff(eq_res, Qx*Qy)=0

since it is not allowed to differentiate any expression with respect to 
non-atomic terms as Qx^2,Qy^2,Qx*Qy, a workaround is applied in the following
consisting in a change of variable.

eq_res : ratsubst( Qxq  , Qx^2  , eq_res  )$
eq_res : ratsubst( Qyq  , Qy^2  , eq_res  )$
eq_res : ratsubst( QxQy , Qx*Qy , eq_res  )$

eq_res;

check: no residual Qx, Qy terms within eq_res

diff( eq_res , Qx);
diff( eq_res , Qy);

ok

final value for the three corrective shear strain energy coefficients: their value
has to be an invariant with respect to the value of the shear components.

linsolve(
    [
        diff(eq_res, Qxq) =0,
        diff(eq_res, Qyq) =0,
        diff(eq_res, QxQy)=0
    ],
    [ chi_x, chi_y, chi_xy]
), globalsolve=true;

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

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