kill(all);

placeholder for data to be provided based on the FE analysis

xxxFILLMExxx: 123.456$

silence the "rat: replaced ... " warnings

ratprint: false$

profile area

area:xxxFILLMExxx;

## Pure unit z elongation loadcase.

ez=1, kx=0, ky=0.
generic RBE2 master node position.

Reaction force Z (EA axial modulus)

EA:3.50079e+7;

Reaction moments X and Y

uezMx : -712641;

uezMy : -6.67053e+7;

master node shift for alignment with the centre of axial elasticity

tx : -uezMy/EA; 
%,numer;

ty : +uezMx/EA; 
%,numer;

Before proceeding further, the master node has to be shifted
accordingly.

## Pure unit kx and ky curvature loadcases

pure, unit kx curvature (ez=0,kx=1,ky=0).
Centroidal RBE2 master node position.

ukxMx:2.75709e+10;

ukxMy:5.12261e+7;

pure, unit ky curvature (ez=0,kx=0,ky=1).
Centroidal RBE2 master node position.

ukyMx:5.12261e+7;

ukyMy:1.12605e+10;

derived moduli

EJxx:ukxMx$
EJyy:ukyMy$
EJxy:-(ukxMy+ukyMx)/2$

[[vals,mults],vecs]:
eigenvectors(
matrix(
        [ EJxx,-EJxy],
        [-EJxy, EJyy]
    )
)$

eigenvalues are returned in unknown order, let's identify the
higher and the lower.

if vals[1] >= vals[2] then (
i1:1,i2:2    
) else (
i1:2,i2:1            
)$

EJ11: vals[i1],numer;

EJ22: vals[i2],numer;

principal directions

dir1 : vecs[i1][1],numer$ /* eigenvalues are assumed to be distinct*/
dir1 : dir1 / sqrt( dir1[1]^2 + dir1[2]^2);

dir2 : vecs[i2][1],numer$
dir2 : dir2 / sqrt( dir2[1]^2 + dir2[2]^2);

Gz rotation of the PROFILE for having the global X,Y aligned with
with 1,2, respectively

rotGz : -atan2(dir1[2],dir1[1]); /*radians*/

rotGz*180/%pi,numer; /*degrees*/

## Pure 
qz=Ez*(a*y+b*x)
global load

Before proceeding further, the global OXYZ reference system must
be centroidal (but NOT necessarily principal of inertia for the
cross section).
Also, 
- skewsymm. constraints are applied at Z=0
- a free-warping RBE2 link is attached to the profile at Z=l,
whose centroidal master node is simply supported in both the x,y,
in-plane directions. I also locked the z rotation to let the line
of action of the shear forces deviating from x=0,y=0, Z=l.
- a w=0 positioning constraint inhibits the residual tz rigid body
motion.

here, Ez=Eref since the material is homogeneous (not true anymore)

first, "u"nit "a" case, b=0 and, a=1, s.t. a dominant Qy shear
load is presumably produced.

Reaction force Y (defines the Qy shear force), in hindsight, it
equates EJxx by definition.

uaQy : 2.7571e10;

EJxx;

Reaction force X (defines the Qx shear force), in hindsight, it
equates -EJxy by definition.

uaQx : -3.6516;

EJxy;

Reaction moment Z  (defines the Mz transport moment, with respect
to the shear center)

uaMz :  xxxFILLMExxx;

second, "u"nit "b" case, a=0 and b=1, s.t. a dominant Qx shear
load is presumably produced.

Reaction force X, In hindsight, it equates EJyy by definition.

ubQx : xxxFILLMExxx;

EJyy;

Reaction force Y. In hindsight, it equates -EJxy by definition.

ubQy : xxxFILLMExxx;

Reaction moment Z;

ubMz : xxxFILLMExxx;

relative position of the shear centre with respect to the current
(presumably, but not necessarily, centroidal) RBE2 master node
location

/home/ebertocchi/shear_center_coordinates_for_maxima_FE_evaluation.png

with reference to the image above, Mz should equate the transport
moment moving from S to C, thus leading to Mz=Qx*(-dy)+Qy*(+dx)
for any arbitrary Qx,Qy loadings.

linsolve(
    [
        uaMz=uaQx*(-dy)+uaQy*(+dx),
        ubMz=ubQx*(-dy)+ubQy*(+dx)
    ],[dx,dy]
),globalsolve=true;
%,numer;

we now translate accordingly the RBE2 master node, which should
become coincident with the profile shear center.

we also scale "a" and "b" s.t. the Qx, Qy shear values equate the
product of 
* the area (steel only), and of 
* a unit tau shear stress

unit nominal stress Qx (uQx) case

linsolve( 
    [
        ax*uaQx + bx*ubQx = area*1, /* Qx desiderato */
        ax*uaQy + bx*ubQy = 0       /* Qy desiderato */
    ],[ax,bx]
),globalsolve=true;
%,numer;

unit nominal stress Qy (uQy) case

linsolve( 
    [
        ay*uaQx + by*ubQx = 0,
        ay*uaQy + by*ubQy = area*1
    ],[ay,by]
),globalsolve=true;
%,numer;

Shear contribution to the elastic strain energy (lineic) density 

reference shear modulus, for a reasonable normalization

G:210000/2/(1+0.3)$

length of the profile modeled portion

l:40$

general quadratic for m

u(Qx,Qy) := 1/2 * 
matrix( [Qx,Qy] )
.
matrix( 
        [chix, chixy/2],
        [chixy/2, chiy]
)
.
matrix( 
    [Qx],
    [Qy] 
) / (G*area);

expand(u(Qx,Qy)),area=_A,G=_G;

I substituite the FE reported values for the strain energy in the
cases
* unit stress Qx, zero Qy
* unit stress Qy, zero Qx
* unit stress both Qx and Qy (coupled case)

linsolve(
    [
        u(1*area,0     )*l = xxxFILLMExxx,
        u(0     ,1*area)*l = xxxFILLMExxx,
        u(1*area,1*area)*l = xxxFILLMExxx
    ],[chix,chixy,chiy]
),globalsolve=true$
[chix,chixy,chiy],numer;

[[vals,mults],vecs]:
eigenvectors(
    matrix( 
        [chix, chixy/2],
        [chixy/2, chiy]
    )
)$

if vals[1] >= vals[2] then (
i1:1,i2:2    
) else (
i1:2,i2:1            
)$

principal chi factors for shear

chi1sh: vals[i1],numer;

chi2sh: vals[i2],numer;

principal directions for shear

dir1sh : vecs[i1][1],numer;

dir2sh : vecs[i2][1],numer;