Modelling assignment 2021
Rigidity of Mavic Cxp 700C Iso 622 32H

Step 1: Experimental rigidity [N,mm]

exp={{0,0},{1.89,0.7*g},{3.49,1.2*g},{5.12,1.5*g},{6.20,2*g},{8.32,2.7*g},{9.74,3.0*g}} /.g->9.81;

F=Normal[Fit[exp,{1, w},w]];

Show[ListPlot[exp,Frame->True],Plot[F,{w,Min[exp[[;;,1]]],Max[exp[[;;,1]]]}],FrameLabel->{"w","F"},PlotLabel->F]

Print["kexp=",D[F,w]]

Step 2: Problem parameters [N,mm]

mea={d->634,b->20.1, t->20.8,Ε->70 10^3,ν->0.32,G->Ε/(2+2 ν)};

Grid[mea/.Rule->List,Alignment->Left,Spacings->{5,1.2},Frame->All]

Step 3: Triangle representation

(*image as the starting point, DO NOT DELETE THIS IMAGE*)
 
image =;

(*boundary representation, scaling, and rotation*)

imsh=ImageMesh[image,Method->"LinearSeparable"];
{pz,py}=ImageDimensions[image];
imsh =Scale[imsh,b/pz /.mea,{0,0}];
imsh =Rotate[imsh,Pi];
Show[imsh,Axes ->True,FrameLabel->{"z[mm]","y[mm]"},Frame->True]

(*triangulation*)

msh=TriangulateMesh[imsh];
Show[msh,Axes ->True,FrameLabel->{"z[mm]","y[mm]"},Frame->True]

Step 4: Moments of area

(*moments of area with respect original origin*)
mao={
Ao->RegionMoment[msh,{0,0}],
Szo->RegionMoment[msh,{0,1}],
Syo->RegionMoment[msh,{1,0}],
Izo->RegionMoment[msh,{0,2}],
Iyo->RegionMoment[msh,{2,0}],
Iyzo->RegionMoment[msh,{1,1}]};

zc=(Syo/Ao /.mao);
yc=(Szo/Ao /.mao);

(*moments of area with respect to area centroid*)
mshc =Translate[msh,-{zc,yc}];
mac={
Ac->RegionMoment[mshc,{0,0}],
Szc->RegionMoment[mshc,{0,1}],
Syc->RegionMoment[mshc,{1,0}],
Izc->RegionMoment[mshc,{0,2}],
Iyc->RegionMoment[mshc,{2,0}],
Iyzc->RegionMoment[mshc,{1,1}]};

sym={A,Sz,Sy,Iz,Iy,Iyz};

Show[mshc,Axes ->True,FrameLabel->{"z[mm]","y[mm]"},Frame->True]

Grid[Transpose[{sym,mao[[;;,1]]//.mao,mac[[;;,1]]//.mac}],Alignment->Left,Spacings->{5,1.2},Frame->All]

Step 5: Equation parameters

geo={A->(Ac /.mac),Ι->(Izc /.mac),R->d/2-yc};
val=Join[mea,geo];

Grid[val /.Rule->List,Alignment->Left,Spacings->{5,1.2},Frame->All]

Step 6: Solution to rigidity

(*left hand sides of differential equations*)

equations:={
NN'[x]-QQ[x]/R==0,
QQ'[x]+NN[x]/R==0,
MM'[x]+QQ[x]==0,
NN[x]==Ε A (u'[x]-v[x]/R),
QQ[x]==G A (v'[x]+u[x]/R-ψ[x]),
MM[x]==Ε Ι ψ'[x]};

(*unknown functions*)

functions:={u[x],v[x],ψ[x],NN[x],QQ[x],MM[x]};

(*boundary conditions*)

clamped[x_]:={u[x]==0,v[x]==0,ψ[x]==0};
loaded[x_]:={NN[x]==0,QQ[x]==-FF,MM[x]==0};

problem:=Join[equations,clamped[0],loaded[3 Pi R/2]];

sol=FullSimplify[DSolve[problem,functions,x][[1]]];

Print["u(L)= ",u[x] /. sol/.x->3 Pi R/2]
Print["v(L)= ",vL=v[x] /. sol/.x->3 Pi R/2]
Print["ψ[L]= ",ψ[x] /. sol/.x->3 Pi R/2]

kmod=Simplify[-FF/vL]

Print["kmod=",kmod//.val]