"Nonlinear analysis of a Von-Mises truss
'(solved by four different numerical methods)
'<h4>Input data</h4>
'Half span -'a = 2m', ''Height -'h = -0.5m
'Cross section - circular pipe with diameter'Φ = 100mm'and thickness't = 4mm
'Material - steel. Modulus of elasticity -'E = 210GPa
'Vertical force -'F_y = 2000kN
#include svg_drawing.cpd
#hide
v = 0mm
L = a/2cm
H = h/cm
H_lin = H
H_nl = H
k = L/10
W = 2*(L + 6*k)
Y = H
F_y,max = 2500kN
v_max = abs(2.26*h)
#def svg$ = '<svg viewbox="'-6*k' 'Y' 'W' 'abs(H)'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro; font-size:7px; width:'420'pt; height:'400*abs(H)/L'pt">
#def thin_style$ = style="stroke:darkcyan; stroke-width:1; fill:none;"
#def dash_style$ = style="stroke:#ddd; stroke-width:1; stroke-dasharray: 7 3; fill:none;"
#def ghost_style$ = style="stroke:#ddd; stroke-width:1; fill:white;"
#def lin_style$ = style="stroke:#ccc; stroke-width:2; fill:none;"
#def pin_style$ = style="stroke:red; stroke-width:1; fill:white;"
#def thick_style$ = style="stroke:darkcyan; stroke-width:3; fill:none;"
#def load_style$ = style="stroke:deepskyblue; stroke-width:1; fill:deepskyblue;"
#show
#def frame$
	line$(0; 0; L; H; dash_style$)
	line$(L; H; 2*L; 0; dash_style$)
	line$(0; 0; L; H_lin; lin_style$)
	line$(L; H_lin; 2*L; 0; lin_style$)
	circle$(L; H_lin; 2.5; ghost_style$)
	line$(0; 0; L; H_nl; main_style$)
	line$(L; H_nl; 2*L; 0; main_style$)
	#hide
	x1 = 0
	y1 = k/10
	#show
	#for i = 1 : 2
		#hide
		δ = 0.8*k
		x2 = x1 - δ
		y2 = y1 - abs(δ)
		x3 = x1 + δ
		y3 = y1 + abs(δ)
		#show
		line$(x2; y3; x3; y3; thick_style$)
		line$(x2; y3; x1; y1; thin_style$)
		line$(x3; y3; x1; y1; thin_style$)
		#hide
		x1 = x1 + 2*L
		#show
	#loop
	'<polygon points="'L','H_nl - 4*δ' 'L','H_nl - 0.32*δ' 'L - 0.32*δ','H_nl - 1.4*δ' 'L + 0.32*δ','H_nl - 1.4*δ' 'L','H_nl - 0.32*δ'" load_style$ />
	texth$(L+20;H_nl-3.2*δ;Fy='F_y'kN)
	dimh$(0; L; (1-2*sign(H_nl))*k; 'a'm)
	dimh$(L; 2*L; (1-2*sign(H_nl))*k; 'a'm)
	#if H_nl < 0
		dimv$(L - k; H_nl; 0; 'h + v'm)
	#else
		dimv$(L - k; 0; H_nl; 'h + v'm)
	#end if
	circle$(0; 0; 2.5; pin_style$)
	circle$(L; H_nl; 2.5; pin_style$)
	circle$(2*L; 0; 2.5; pin_style$)
#end def
#val
svg$
frame$
'</svg>
#equ
'<h4>Linear Solution</h4>
'Strut length -'l = srss(a; h)
'Area -'A = π*(Φ^2 - (Φ - 2*t)^2)/4|cm^2
'Axial stiffness -'EA = E*A
'Axial forces in bars -'N(F) = F*l/(2*h)', 'N(F_y)
'Vertical displacement -'v_l(F) = -sqrt((l + (N(F)*l)/EA)^2 - a^2) - h|mm
v_l(F_y)
'<h4>Nonlinear Solution</h4>
'The vertical displacement is the only unknown - <var>v</var> = ?
'We will use 3-rd order geometric nonlinearity theory for the solution. The equilibrium equations are then derived for the deformed state of the structure, as follows:
'Length and elongation in deformed state
l′(v) = sqrt(a^2 + (h + v)^2)','Δl(v) = l′(v) - l
'Horizontal reaction -'F_x(v) = EA*(Δl(v)/l)*(a/l′(v))
'Vertical reaction -'F_y(v) = EA*(Δl(v)/l)*((h + v)/l′(v))
'Vertical reaction derivative -'F′_yv(v) = EA*(1/l - a^2/l′(v)^3)
'<h4>1. Fixed point iteration method</h4>
'Relative strain -'ε = F_y/(2*EA)
'Relative precision -'δ_max = 10^-5
'Initial value -'v_0 = 1mm
'We express the unknown vertical displacement at the middle joint as a function of the vertical force:
v = sqrt(1/(1/l - ε/(h + v_0))^2 - a^2) - h|mm
#hide
n = 0
#repeat 100
	n = n + 1
	v_0 = v
	v = sqrt(1/(1/l - ε/(h + v_0))^2 - a^2) - h|mm
	δ = abs(v - v_0)
	#if δ < δ_max*abs(v)
		#break
	#end if
#loop
#show
'After calculating the above expression iteratively'n'times, we get:
v
'Relative error -'δ = abs(v - v_0)/abs(v)
'<h4>2. Newton-Raphson′s method</h4>
'Initial value -'v_0 = 0mm
'We repeatedly calculate the following expression:
v = v_0 - (2*F_y(v_0) - F_y)/F′_yv(v_0)
#hide
n = 0
#repeat 100
	n = n + 1
	v_0 = v
	v = v_0 - (2*F_y(v_0) - F_y)/(2*F′_yv(v_0))
	δ = abs(v - v_0)
	#if δ < δ_max*abs(v)
		#break
	#end if
#loop
#show
'After'n'iterations we get:'v
'Relative error -'δ = abs(v - v_0)/abs(v)
'<h4>3. Secant method</h4>
'Slope reduction factor -'α = 1
'Initial value -'v_0 = 0mm
'Force value -'F_y0 = 2*F_y(v_0)
'We calculate the first approximation using Newton-Raphson′s method
v_1 = v_0 - α*((F_y0 - F_y)/(2*F′_yv(v_0)))
'Force value -'F_y1 = 2*F_y(v_1)
'The next approximation is evaluated by the formula:
v_2 = v_1 - α*(F_y1 - F_y)*((v_1 - v_0)/(F_y1 - F_y0))
'We continue the calculations iteratively until we reach convergence.
#hide
n = 0
#repeat 100
	n = n + 1
	v_0 = v_1
	v_1 = v_2
	F_y0 = 2*F_y(v_0)
	F_y1 = 2*F_y(v_1)
	v_2 = v_1 - α*(F_y1 - F_y)*((v_1 - v_0)/(F_y1 - F_y0))
	δ = abs(v_2 - v_1)
	#if δ < δ_max*abs(v_2)
		#break
	#end if
#loop
#show
'After'n'iterations we get:'v_2
'Relative error -'δ = abs(v_2 - v_1)/abs(v_2)
'<h4>4. Solution with Calcpad <small>(modified Anderson-Bjork′s method)</small></h4>
v = $Root{2*F_y(v) = F_y @ v = 0mm : 200m}
'System behavior graph (force-displacement)
'<!--'PlotWidth = 400','PlotHeight = 180','ReturnAngleUnits = 1'-->
#def plot$
	#hide
	v_2 = v
	F′ = 2*F′_yv(v)
	dv = min(2000kN/F′; v_max/6)
	v_1(v) = min(max(0mm; v_2 + (1 - 2*v/v_max)*dv); v_max)
	F_1(v) = F_y + F′*(v_1(v) - v_2)
	#show
	$Plot{F_y & 2*F_y(v) & v_1(v)|F_1(v) & v_l(F_y)|F_y & v_2|F_y @ v = 0mm : v_max}
#end def
plot$
'<h4>Results</h4>
'Axial forces in bars -'N = Δl(v)/l*EA
'Rotation angle -'α = atan2(l; v)
'Reactions in supports
'Horizontal -'R_x = F_x(v)'='N*cos(α)
'Vertical -'R_y = F_y(v)'='N*sin(α)
#hide
L = a/2cm
k = L/10
δ = 0.7*k
#show
#val
#def results$
	#hide
	N = Δl(v)/l*EA
	α = atan2(l; v)
	R_x = F_x(v)'='N*cos(α)
	R_y = F_y(v)'='N*sin(α)
	H_lin = (h + v_l(F_y))/cm
	H_nl = (h + v)/cm
	Y = H + 2*k
	#show
	svg$
	'<polygon points="'0','0' '0','-4*δ' '-0.32*δ','-2.8*δ' '0.32*δ','-2.8*δ' '0','-4*δ'" load_style$/>
	'<polygon points="'2*L','0' '2*L','-4*δ' '2*L - 0.32*δ','-2.8*δ' '2*L + 0.32*δ','-2.8*δ' '2*L','-4*δ'" load_style$/>
	'<polygon points="'-0.32*δ','0' '-4*δ','0' '-2.8*δ','-0.32*δ' '-2.8*δ','0.32*δ' '-4*δ','0'" load_style$/>
	'<polygon points="'2*L + 0.32*δ','0' '2*L + 4*δ','0' '2*L + 2.8*δ','-0.32*δ' '2*L + 2.8*δ','0.32*δ' '2*L + 4*δ','0'" load_style$/>
	texth$(-4*δ;-δ;Rx='R_x'kN)
	texth$(2*L+4*δ;-δ;Rx='R_x'kN)
	texth$(3.2*δ;-3.2*δ;Ry='R_y'kN)
	texth$(2*L-3.2*δ;-3.2*δ;Ry='R_y'kN)
	texth$(L/2;0;N='N'kN)
	texth$(3*L/2;0;N='N'kN)
	frame$
	'</svg>
#end def
results$
#hide
n = 51
#show
'Animation&emsp;<button id="start">Start</button>
'<style>#start{height:1.8em; padding-top:0;} .fr{display:none;}</style>
#for j = 1 : n
	#hide
	F_y = F_y,max*(j - 1)/(n - 1)
	v = $Root{2*F_y(v) = F_y @ v = 0mm : 200m}
	#show
	'<div class="fr" id="fr-'j'">
	plot$
	results$
	'</div>
#loop
'<script>$("#fr-1").show();$("#start").click(function(){$(this).hide();var i='n';var fr=$("#fr-"+i);setInterval(function(){fr.hide();if(++i>'n')i=1;fr=$("#fr-"+i);fr.show();}, 100);});</script>
#equ