Beam equations 2023

xyz-system

(* limit α -> Infinity gives Bernoulli solution*)

(*equilibrium equations*)
equilibrium[xyz]:={
Nx'[x]+bx,
Qy'[x]+by,
Qz'[x]+bz,
Tx'[x]+cx,
My'[x]-Qz[x]+cy,
Mz'[x]+Qy[x]+cz};

(*constitutive equations*)
constitutive[xyz]:={
Nx[x]-Ε (A u'[x]-Sz ψ'[x]+Sy θ'[x]),
Qy[x]-G (A α (v'[x]-ψ[x])-Sy ϕ'[x]),
Qz[x]-G (A α (w'[x]+θ[x])+Sz ϕ'[x]),
Tx[x]-G (-Sy α (v'[x]-ψ[x])+Sz α (w'[x]+θ[x])+Irr ϕ'[x]),
My[x]-Ε (Sy u'[x]-Iyz ψ'[x]+Iyy θ'[x]),
Mz[x]-Ε (-Sz u'[x]+Izz ψ'[x]-Izy θ'[x])};

(*quantities*)
displacements[xyz]:={u[x],v[x],w[x],ϕ[x],θ[x],ψ[x]};
forces[xyz]:={Nx[x],Qy[x],Qz[x],Tx[x],My[x],Mz[x]};
functions[xyz]:=Join[displacements[xyz],forces[xyz]];

(*some common boundary conditions*)
clamped[xyz][x_]:={u[x],v[x],w[x],ϕ[x],θ[x],ψ[x]};
free[xyz][x_]:={Nx[x],Qy[x],Qz[x],Tx[x],My[x],Mz[x]};
simply[xyz][x_]:={u[x],v[x],w[x],Tx[x],My[x],Mz[x]};
slide[xyz][x_]:={Nx[x],v[x],w[x],ϕ[x],My[x],Mz[x]};

snb-system

(* limit α -> Infinity gives Bernoulli solution*)

equilibrium[snb]:={
Ns'[s]-Qn[s] κ+bs,
Qn'[s]+Ns[s] κ-Qb[s] τ+bn,
Qb'[s]+Qn[s] τ+bb,
Ts'[s]-Mn[s] κ+cs,
Mn'[s]+Ts[s] κ-Mb[s]*τ-Qb[s]+cn,
Mb'[s]+Mn[s] τ+Qn[s]+cb};

constitutive[snb]:={
Ns[s]-Ε A (u'[s]-v[s] κ)+Ε Sb (ψ'[s]+θ[s] τ)-Ε Sn (θ'[s]+ϕ[s] κ-ψ[s] τ),
Qn[s]-G A α (v'[s]+u[s] κ-w[s] τ-ψ[s])+G Sn (ϕ'[s]-θ[s] κ),
Qb[s]-G A α (w'[s]+v[s] τ+θ[s])-G Sb (ϕ'[s]-θ[s] κ),
Ts[s]-G Irr (ϕ'[s]-θ[s] κ)-G Sb α (w'[s]+v[s] τ+θ[s])+G Sn α (v'[s]+u[s] κ-w[s] τ-ψ[s]),
Mn[s]-Ε Inn (θ'[s]+ϕ[s] κ-ψ[s] τ)-Ε Sn (u'[s]-v[s] κ)+Ε Ibn (ψ'[s]+θ[s] τ),
Mb[s]-Ε Ibb (ψ'[s]+θ[s] τ)+Ε Sb (u'[s]-v[s] κ)+Ε Inb (θ'[s]+ϕ[s] κ-ψ[s] τ)};

clamped[snb][s_]:={u[s],v[s],w[s],ϕ[s],θ[s],ψ[s]};
free[snb][s_]:={Ns[s],Qn[s],Qb[s],Ts[s],Mn[s],Mb[s]};
simply[snb][s_]:={u[s],v[s],w[s],Ts[s],Mn[s],Mb[s]};

displacements[snb]:={u[s],v[s],w[s],ϕ[s],θ[s],ψ[s]};
forces[snb]:={Ns[s],Qn[s],Qb[s],Ts[s],Mn[s],Mb[s]};
functions[snb]:=Join[displacements[snb],forces[snb]];

zxy-system (shear and torsion corrections and symmetry)

(*equilibrium equations*)
equilibrium[zxy]:={
Fz'[z]+bz,
Fx'[z]+bx,
Fy'[z]+by,
Mz'[z]+cz,
Mx'[z]-Fy[z]+cx,
My'[z]+Fx[z]+cy};

(*constitutive equations*)
constitutive[zxy]:={
Fz[z]-Ε A uz'[z],
Fx[z]-κ1 G A (ux'[z]-θy[z])-κ4 G Sqrt[A Irr] θy'[z],
Fy[z]-κ2 G A(uy'[z]+θx[z]),
Mz[z]-κ3 G Irr θz'[z]- κ4 G Sqrt[A Irr]( ux'[z]- θy[z]),
Mx[z]-Ε Iyy θx'[z],
My[z]-Ε Ixx θy'[z]};

(*quantities*)
displacements[zxy]:={uz[z],ux[z],uy[z],θz[z],θx[z],θy[z]};
forces[zxy]:={Fz[z],Fx[z],Fy[z],Mz[z],Mx[z],My[z]};
functions[zxy]:=Join[displacements[zxy],forces[zxy]];

(*some common boundary conditions*)
clamped[zxy][z_]:={uz[z],ux[z],uy[z],θz[z],θx[z],θy[z]};
free[zxy][z_]:={Fz[z],Fx[z],Fy[z],Mz[z],Mx[z],My[z]};

xyz-system (torsion corrections)

(*equilibrium equations Bernoulli model*)
(*area centroid origin, z-down,x-axis, y-to_me*)
(*ey,ez shear centroid position*)

equilibrium[xyz]:={
Fx'[x]+fx,
Fy'[x]+fy,
Fz'[x]+fz,
Mx'[x],
My'[x]-Fz[x],
Mz'[x]+Fy[x]};

(*constitutive equations*)
constitutive[xyz]:={
Fx[x]-EE AA ux'[x],
uy'[x]-θz[x],
uz'[x]+θy[x],
Mx[x]-GG J θx'[x],
My[x]-EE Iyy θy'[x],
Mz[x]-EE Izz θz'[x]};
(*quantities*)
displacements[xyz]:={ux[x],uy[x],uz[x],θx[x],θy[x],θz[x]};
forces[xyz]:={Fx[x],Fy[x],Fz[x],Mx[x],My[x],Mz[x]};
functions[xyz]:=Join[displacements[xyz],forces[xyz]];

(*some common boundary conditions*)
clamped[xyz][x_]:={ux[x],uy[x],uz[x],θx[x],θy[x],θz[x]};
free[xyz][x_]:={Fx[x],Fy[x],Fz[x],Mx[x],My[x],Mz[x]};

sol=DSolve[Thread[Join[equilibrium[xyz],constitutive[xyz],
clamped[xyz][0]-{ux1,uy1,uz1,θx1,θy1,θz1},
clamped[xyz][h0]-{ux2,uy2,uz2,θx2,θy2,θz2}]==0],functions[xyz],x][[1]];

for={Fx[x],Fy[x],Fz[x],Mx[x],My[x],Mz[x]}/.sol;
FF=Simplify[Join[-for/.x->0,for/.x->h0]];
(*MatrixForm[D[FF,{Join[{ux1,uy1,uz1,θx1,θy1,θz1},{ux2,uy2,uz2,θx2,θy2,θz2}]}]]*)

FullSimplify[Simplify[-Join[{δux[1,0],δuy[1,0],δuz[1,0],δθx[1,0],δθy[1,0],δθz[1,0]},{δux[2,0],δuy[2,0],δuz[2,0],δθx[2,0],δθy[2,0],δθz[2,0]}].(0*FF+(FF /.{fx->0,fy->0,fz->0}))]/.
Thread[{ux1,uy1,uz1,θx1,θy1,θz1}->{ux[1,0],uy[1,0],uz[1,0],θx[1,0],θy[1,0],θz[1,0]}]/.
Thread[{ux2,uy2,uz2,θx2,θy2,θz2}->{ux[2,0],uy[2,0],uz[2,0],θx[2,0],θy[2,0],θz[2,0]}]]

Examples

Lecture problem 4.1 (cantilever xyz)

rule={bx->0,by ->b, bz->0,cx->0,cy->0,cz->0,Sy->0,Sz->0,  Iyz->0,Izy->0,Iyy->Ι,Izz->Ι,Irr->2 Ι, α->1};

equs=Join[equilibrium[xyz],constitutive[xyz],clamped[xyz][0],free[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["u[L] = ",Simplify[(u[x] /. sol  /.x->L)]]
Print["v[L] = ",Simplify[(v[x] /. sol /.x->L)]]
Print["w[L] = ",Simplify[(w[x] /. sol /.x->L)]]
Print["ψ[L] = ",Simplify[(ψ[x] /. sol /.x->L)] ]

Lecture problem 4.2 (cantilever xyz)

rule={bx->0,by ->b, bz->0,cx->0,cy->0,cz->0,Sy->0,Sz->0,  Iyz->0,Izy->0,Iyy->Ι,Izz->Ι};

equs=Join[equilibrium[xyz],constitutive[xyz],clamped[xyz][0],clamped[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["u[L] = ",Simplify[Limit[u[x] /. sol , α->1]]]
Print["v[L] = ",Simplify[Limit[v[x] /. sol  ,α->1]]]
Print["ψ[L] = ",Simplify[Limit[ψ[x] /. sol, α->1]] ]

Lecture problem 4.3 (cantilever xyz)

rule={bx->0,by ->0, bz->ρgA,cx->0,cy->0,cz->0};

equs=Join[equilibrium[xyz],free[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],forces[xyz],x][[1]];

Print["Nx(x)= ",FullSimplify[Nx[x] /. sol ]]
Print["Qy(x)= ",FullSimplify[Qy[x] /. sol]]
Print["Qz(x)= ",FullSimplify[Qz[x] /. sol]]
Print["Tx(x)= ",FullSimplify[Tx[x] /. sol]]
Print["My(x)= ",FullSimplify[My[x] /. sol]]
Print["Mz(x)= ",FullSimplify[Mz[x] /. sol]]

Lecture problem 4.4 (cantilever xyz offset)

rule={bx->0,by ->b, bz->0,cx->0,cy->0,cz->0,Sy->0,Sz->-r A,  Iyz->0,Izy->0,Iyy->Ι,Izz->Ι+A r^2};

equs=Join[equilibrium[xyz],constitutive[xyz],clamped[xyz][0],free[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["u[L] = ",Simplify[Limit[u[x] /. sol /.x->L, α->Infinity]]]
Print["v[L] = ",Simplify[Limit[v[x] /. sol /.x->L ,α->Infinity]]]
Print["ψ[L] = ",Simplify[Limit[ψ[x] /. sol/.x->L, α->Infinity]] ]

Lecture problem 4.5 (half-circle snb)

rule={bs->0,bn->0,bb->0,cs->0,cn->0,cb->0,τ->0,κ->1/R};

equs=Join[equilibrium[snb],free[snb][Pi R]-{0,0,P,0,0,0}]  /.rule;

sol=DSolve[Thread[equs==0],forces[snb],s][[1]];

Print["Ns(s)= ",FullSimplify[Ns[s] /. sol ]]
Print["Qn(s)= ",FullSimplify[Qn[s] /. sol]]
Print["Qb(s)= ",FullSimplify[Qb[s] /. sol]]
Print["Ts(s)= ",FullSimplify[Ts[s] /. sol] ]
Print["Mn(s)= ",FullSimplify[Mn[s]  /. sol]]
Print["Mb(s)= ",FullSimplify[Mb[s] /. sol]]

Lecture problem 4.6 (quarter-circle snb)

rule={bs->0,bn->0,bb->0,cs->0,cn->0,cb->0,τ->0,κ->1/R,
Sn->0,Sb->0,Inn->Ι,Ibb->Ι,Ibn->0,Inb->0};

equs=Join[equilibrium[snb],constitutive[snb],clamped[snb][0],free[snb][Pi R/2]-{P,0,0,0,0,0}]  /.rule;

sol=DSolve[Thread[equs==0],functions[snb],s][[1]];

Print["u(s)= ",FullSimplify[Limit[u[s] /. sol /.s->Pi R/2,α->Infinity]]]
Print["v(s)= ",FullSimplify[Limit[v[s] /. sol/.s->Pi R/2,α->Infinity]]]
Print["ψ(s)= ",FullSimplify[Limit[ψ[s] /. sol/.s->Pi R/2,α->Infinity] ]]

Lecture problem 4.7 (half-circle snb)

rule={bs->0,bn->0,bb->ρgA,cs->0,cn->0,cb->0,τ->0,κ->1/R,
Sn->0,Sb->0,Inn->Ι,Ibb->Ι,Ibn->0,Inb->0,Irr->2Ι};

equs=Join[equilibrium[snb],constitutive[snb],clamped[snb][0],clamped[snb][Pi*R]]  /.rule;

sol=DSolve[Thread[equs==0],functions[snb],s][[1]];

Print["ϕ(L/2)= ",FullSimplify[Limit[ϕ[s] /. sol /.s->Pi R/2,α->Infinity]]]

Print["w(L/2)= ",FullSimplify[Limit[w[s] /. sol /.s->Pi R/2,α->Infinity]]]

Example problem 2

rule={bx->0,by ->0, bz->b(1-x/L),cx->0,cy->0,cz->0,Sy->0,Sz->0,  Iyz->0,Izy->0,Iyy->Ι};

equs=Join[equilibrium[xyz],constitutive[xyz],clamped[xyz][0],free[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["Qz[x] = ",Simplify[Limit[(Qz[x] /. sol ),α->Infinity]]]
Print["My[x] = ",Simplify[Limit[(My[x] /. sol ),α->Infinity]]]
Print["w[x] = ",Simplify[Limit[(w[x] /. sol ),α->Infinity]]]

Example problem 3

rule={bx->0,by ->0, bz->0,cx->0,cy->0,cz->0,Sy->0,Sz->0,  Iyz->0,Izy->0,Iyy->Ι};

equs=Join[equilibrium[xyz],constitutive[xyz],clamped[xyz][0],free[xyz][L]-{0,0,P,0,0,0}] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["w[x] = ",Simplify[Limit[(w[x] /. sol ),α->Infinity]]]
Print["θ[x] = ",Simplify[Limit[(θ[x] /. sol ),α->Infinity]]]

Example problem 4

rule={bx->0,by ->0, bz->b,cx->0,cy->0,cz->0,Sy->0,Sz->0,  Iyz->0,Izy->0,Iyy->Ι};

equs=Join[equilibrium[xyz],constitutive[xyz],simply[xyz][0],slide[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["w[x] = ",Simplify[(w[x] /. sol ) /.{ x->L/2,α->1}]]
Print["θ[x] = ",Simplify[(θ[x] /. sol )/. { x->L/2,α->1}] ]

Example problem 5

rule={bx->0,by ->b, bz->b,cx->0,cy->0,cz->0,Sy->0,Sz->0,  Iyz->0,Izy->0,Iyy->Ι,Izz->Ι};

equs=Join[equilibrium[xyz],constitutive[xyz],simply[xyz][0],slide[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["v[x] = ",Simplify[Limit[(v[x] /. sol ), α->Infinity]]]

Example problem 6

rule={bs->0,bn->0,bb->0,cs->0,cn->0,cb->0,τ->4 Pi/L,κ->0};

equs=Join[equilibrium[snb],free[snb][L]-{0,0,P,0,0,0}]/.rule ;

sol=DSolve[Thread[equs==0],forces[snb],s][[1]];

Print["Ns(s)= ",FullSimplify[Ns[s] /. sol ]]
Print["Qn(s)= ",FullSimplify[Qn[s] /. sol]]
Print["Qb(s)= ",FullSimplify[Qb[s] /. sol]]
Print["Ts(s)= ",FullSimplify[Ts[s] /. sol]]
Print["Mn(s)= ",FullSimplify[Mn[s] /. sol]]
Print["Mb(s)= ",FullSimplify[Mb[s] /. sol]]

Example problem 7

rule={bs->0,bn->0,bb->0,cs->0,cn->0,cb->0,τ->0,κ->1/R};

equs=Join[equilibrium[snb],free[snb][Pi R/2]-{0,0,-P,0,0,0}] /.rule;

sol=DSolve[Thread[equs==0] ,forces[snb],s][[1]];

Print["Ns(s)= ",FullSimplify[Ns[s] /. sol ]]
Print["Qn(s)= ",FullSimplify[Qn[s] /. sol]]
Print["Qb(s)= ",FullSimplify[Qb[s] /. sol]]
Print["Ts(s)= ",FullSimplify[Ts[s] /. sol]]
Print["Mn(s)= ",FullSimplify[Mn[s] /. sol]]
Print["Mb(s)= ",FullSimplify[Mb[s] /. sol]]

Example problem 8

rule={bs->0,bn->0,bb->0,cs->0,cn->0,cb->0,τ->0,κ->1/R};

equs=Join[equilibrium[snb],free[snb][0]-{0,0,0,P,0,0}] /.rule;

sol=DSolve[Thread[equs==0] ,forces[snb],s][[1]];

Print["Ns(s)= ",FullSimplify[Ns[s] /. sol ]]
Print["Qn(s)= ",FullSimplify[Qn[s] /. sol]]
Print["Qb(s)= ",FullSimplify[Qb[s] /. sol]]
Print["Ts(s)= ",FullSimplify[Ts[s] /. sol]]
Print["Mn(s)= ",FullSimplify[Mn[s] /. sol]]
Print["Mb(s)= ",FullSimplify[Mb[s] /. sol]]

Example problem 9

rule={bs->0,bn->0,bb->b,cs->0,cn->0,cb->0,τ->0,κ->1/R};

equs=Join[equilibrium[snb],free[snb][Pi R/2]] /.rule;

sol=DSolve[Thread[equs==0] ,forces[snb],s][[1]];

Print["Ns(s)= ",FullSimplify[Ns[s] /. sol ]]
Print["Qn(s)= ",FullSimplify[Qn[s] /. sol]]
Print["Qb(s)= ",FullSimplify[Qb[s] /. sol]]
Print["Ts(s)= ",FullSimplify[Ts[s] /. sol]]
Print["Mn(s)= ",FullSimplify[Mn[s] /. sol]]
Print["Mb(s)= ",FullSimplify[Mb[s] /. sol]]

Example problem 10

rule={bs->0,bn->0,bb->0,cs->0,cn->0,cb->0,τ->0,κ->1/R,
Sn->0,Sb->0,Inn->Ι,Ibb->Ι,Ibn->0,Inb->0};

equs=Join[equilibrium[snb],constitutive[snb],clamped[snb][0],free[snb][Pi R/2]-{0,P,0,0,0,0}]  /.rule;

sol=DSolve[Thread[equs==0],functions[snb],s][[1]];

Print["u(s)= ",FullSimplify[u[s] /. sol /.α->1 ]]
Print["v(s)= ",FullSimplify[v[s] /. sol/.α->1 ]]
Print["ψ(s)= ",FullSimplify[ψ[s] /. sol/.α->1 ]]

Assignment 1

rule={bx->0,by ->b, bz->0,cx->0,cy->0,cz->0};

equs=Join[equilibrium[xyz],free[xyz][L]-{P,0,0,0,0,0}]/.rule;  

sol=DSolve[Thread[equs==0],forces[xyz],x][[1]];

Print["Nx[x] = ",Simplify[(Nx[x] /. sol )]]
Print["Qy[x] = ",Simplify[(Qy[x] /. sol )]]
Print["Mz[x] = ",Simplify[(Mz[x] /. sol )]]

Assignment 2

rule={bs->0,bn->0,bb->0,cs->0,cn->0,cb->0,τ->0,κ->1/R};

equs=Join[equilibrium[snb],free[snb][Pi R]-{-P,0,0,0,0,-H}]/.rule ;

sol=DSolve[Thread[equs==0],forces[snb],s][[1]];

Print["Ns(s)= ",FullSimplify[Ns[s] /. sol ]]
Print["Qn(s)= ",FullSimplify[Qn[s] /. sol]]
Print["Mb(s)= ",FullSimplify[Mb[s] /. sol]]

Assignment 4

rule={bx->0,by ->0, bz->b,cx->0,cy->0,cz->0,Sy->0,Sz->0,  Iyz->0,Izy->0,Iyy->Ι,Izz->Ι,Irr->2Ι};

equs=Join[equilibrium[xyz],constitutive[xyz],clamped[xyz][0],free[xyz][L]] /.rule;  

sol=DSolve[Thread[equs==0],functions[xyz],x][[1]];

Print["N[x] = ",Simplify[(Nx[x] /. sol )]]
Print["Qz[x] = ",Simplify[(Qz[x] /. sol )]]
Print["My[x] = ",Simplify[(My[x] /. sol )]]
Print["u[x] = ",Simplify[(u[x] /. sol )]]
Print["w[x] = ",Simplify[(w[x] /. sol )]]
Print["θ[x] = ",Simplify[(θ[x] /. sol )]]