"Elastic Analysis of Deep Beam
'<hr />
#rad
'<h3>Input data</h3> 
'<table><tr><td>
'Beam length -'l = ? {4}m
'Height -'h = ? {2}m
'Load -'q = ? {100}*(kN/m)
'Length of load -'b = ? {0.8}m
'Distance to load start -'a = ? {1.6}m
'Length of support -'d = ? {0.4}m
'Number of iterations -'N = ? {21}
'</td><td><img src="https://calcpad.eu/media/mechanics/elastic/deep-beam-fea.png" style="height:160pt; margin-left:50pt;"></td></tr></table>
#hide
c = a + b
q(x) = if((x < a) + (x > c); 0kN/m; q)
α(n) = n*π/l
αh(n) = n*π*h/l
q_n(n) = 2*q/(α(n)*l)*(cos(α(n)*a) - cos(α(n)*c))
r = q*b/(2*d)
e = l - d
r(x) = if((x < d) + (x > e); r; 0kN/m)
r_n(n) = 2*r/(α(n)*l)*(1 - cos(α(n)*l) - (cos(α(n)*d) - cos(α(n)*e)))
q_series(x) = $Sum{q_n(n)*sin(α(n)*x) @ n = 1 : N}
r_series(x) = $Sum{r_n(n)*sin(α(n)*x) @ n = 1 : N}
PlotWidth = 400
PlotHeight = 150
#post
'Original and Fourier functions for the load
$Plot{q(x) & q_series(x) @ x = 0m : l}
'Original and Fourier functions for the reaction in supports
$Plot{r(x) & r_series(x) @ x = 0m : l}
#hide
sh(n) = sinh(αh(n))
cth(n) = coth(αh(n))
β(n) = αh(n)/sh(n)
Δ(n) = α(n)*(1 - β(n))*(1 + β(n))
A(n) = r_n(n)/α(n)^2
B(n) = q_n(n)*(1/sh(n) + β(n)*cth(n))/(α(n)*Δ(n)) - r_n(n)*(cth(n) + β(n)/sh(n))/(α(n)*Δ(n))
C(n) = -B(n)*α(n)
D(n) = (q_n(n)*β(n) - r_n(n))/Δ(n)
ch(n; x) = cosh(α(n)*x)
th(n; x) = tanh(α(n)*x)
Y(n; y) = ch(n; y)*(A(n) + B(n)*th(n; y) + C(n)*y + D(n)*y*th(n; y))
Y_1(n; y) = ch(n; y)*(α(n)*(A(n)*th(n; y) + B(n)) + C(n)*(1 + α(n)*y*th(n; y)) + D(n)*(th(n; y) + α(n)*y))
Y_2(n; y) = ch(n; y)*α(n)*(A(n)*α(n) + B(n)*α(n)*th(n; y) + C(n)*(2*th(n; y) + α(n)*y) + D(n)*(2 + α(n)*y*th(n; y)))
Y(N; 0m) - r_n(N)/α(N)^2
Y_1(N; 0m)
Y_1(N; h)
Y(N; h) - q_n(N)/α(N)^2
PlotStep = 5
PlotHeight = PlotWidth*h/l
#post
'<h3>Calculation of stresses</h3> 
σ_x(x; y) = $Sum{Y_2(n; y)*sin(α(n)*x) @ n = 1 : N}
'<!--'k = round(abs(σ_x(l/2; h))/20)*50'-->
$Map{σ_x(x; y) @ x = 0m : l & y = 0m : h}
'Plot for σ<sub>x</sub> at <var>x</var> = <var>l</var>/2
$Plot{σ_x(l/2; y)|y & -k|0m & k|0m @ y = 0m : h}
'Bottom value -'σ_x(l/2; 0m)
'Top value -'σ_x(l/2; h)
σ_y(x; y) = -$Sum{α(n)^2*Y(n; y)*sin(α(n)*x) @ n = 1 : N}
$Map{σ_y(x; y) @ x = 0m : l & y = 0m : h}
'Top value -'σ_y(l/2; h)
τ_xy(x; y) = -$Sum{α(n)*Y_1(n; y)*cos(α(n)*x) @ n = 1 : N}
$Map{τ_xy(x; y) @ x = 0m : l & y = 0m : h}
'Value at 3/4 of span -'τ_xy(3*l/4; h/2)
σ_m(x; y) = 0.5*(σ_x(x; y) + σ_y(x; y))
Δσ(x; y) = 0.5*sqr((σ_x(x; y) - σ_y(x; y))^2 + 4*τ_xy(x; y)^2)
σ_1(x; y) = σ_m(x; y) + Δσ(x; y)
$Map{σ_1(x; y) @ x = 0m : l & y = 0m : h}
σ_2(x; y) = σ_m(x; y) - Δσ(x; y)
$Map{σ_2(x; y) @ x = 0m : l & y = 0m : h}