Hello,
I’m trying to implement a 2D lumped beam model with modelingtoolkit. Using the small angle approximation, I managed to have something working with 5 elements. But increasing the number of elements to 6 make ODEProblem last forever …any clues ?
Plus when I try to use exact version with cos and sine (you can try by uncommenting the other rigidlink definition) I received the following error msg:
ERROR: MethodError: no method matching var_derivative!(::TearingState{NonlinearSystem}, ::Int64)
Any ideas on the two topics greatly appreciated !
using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit: t_nounits as t, D_nounits as D
import DifferentialEquations
using Plots
#-----------------------------------custom 2D mech library----------------------------------
@connector Frame begin
x(t), [description = "x position of frame"]
y(t), [description = "y position of frame"]
phi(t), [description = "Rotation angle of frame"]
tau(t), [connect=Flow, description = "Cut torque in frame"]
fx(t), [connect=Flow, description = "Cut force x in frame"]
fy(t), [connect=Flow, description = "Cut force y in frame"]
end
@mtkmodel Fixed begin
@components begin
frame = Frame()
end
@parameters begin
x0 = 0.0, [description = "Fixed x position"]
y0 = 0.0, [description = "Fixed y position"]
phi0 = 0.0, [description = "Fixed y position"]
end
@equations begin
frame.x ~ x0
frame.y ~ y0
frame.phi ~ phi0
end
end
@mtkmodel PartialCompliantRotational begin
@components begin
frame_a = Frame()
frame_b = Frame()
end
@variables begin
phi_rel(t), [description = "Relative rotation angle between frames", guess = 0.0]
w_rel(t), [description = "Relative angular velocity between flanges", guess = 0.0]
tau(t), [description = "Torque between frames", guess = 0.0]
end
@equations begin
phi_rel ~ frame_b.phi - frame_a.phi
D(phi_rel) ~ w_rel
frame_b.tau ~ tau
frame_a.tau ~ -tau
frame_a.x ~ frame_b.x
frame_a.y ~ frame_b.y
frame_a.fx ~ -frame_b.fx
frame_a.fy ~ -frame_b.fy
end
end
@mtkmodel SpringDamperRotational begin
@extend phi_rel, w_rel, tau = partial_comp = PartialCompliantRotational()
@variables begin
tau_c(t), [description = "Spring torque"]
tau_d(t), [description = "Damper torque"]
end
@parameters begin
d, [description = "Damping constant"]
c, [description = "Spring constant"]
phi_rel0 = 0.0, [description = "Unstretched spring angle"]
end
@equations begin
tau_c ~ c * (phi_rel - phi_rel0)
tau_d ~ d * w_rel
tau ~ tau_c + tau_d
end
end
@mtkmodel MassInertia begin
@components begin
frame = Frame()
end
@parameters begin
m = 1.0, [description = "Mass kg"]
i = 10.0, [description = "Inertia kg⋅m²"]
end
@variables begin
x(t), [description = "Absolute x position", guess = 0.0]
vx(t), [description = "Absolute x velocity", guess = 0.0]
ax(t), [description = "Absolute x acceleration", guess = 0.0]
y(t), [description = "Absolute y position", guess = 0.0]
vy(t), [description = "Absolute y velocity", guess = 0.0]
ay(t), [description = "Absolute y acceleration", guess = 0.0]
phi(t), [description = "Absolute angular position", guess = 0.0]
w(t), [description = "Absolute angular velocity", guess = 0.0]
a(t), [description = "Absolute angular acceleration", guess = 0.0]
end
@equations begin
frame.x ~ x
D(x) ~ vx
D(vx) ~ ax
frame.y ~ y
D(y) ~ vy
D(vy) ~ ay
frame.phi ~ phi
D(phi) ~ w
D(w) ~ a
m*ax ~ frame.fx
m*ay ~ frame.fy
i*a ~ frame.tau
end
end
# @mtkmodel RigidLink begin
# @parameters begin
# l = l
# phi_0 = 0, [description = "link angle"]
# end
# @components begin
# frame_a = Frame()
# frame_b = Frame()
# end
# @equations begin
# # geometric
# frame_b.x ~ frame_a.x + l*cos(phi_0+frame_a.phi)
# frame_b.y ~ frame_a.y + l*sin(phi_0+frame_a.phi)
# frame_b.phi ~ frame_a.phi+phi_0
# # torque
# frame_a.tau ~ -(-frame_b.fx*l*sin(phi_0+frame_a.phi) +
# frame_b.fy*l*cos(phi_0+frame_a.phi)) - frame_b.tau
# # forces
# frame_a.fx ~ -frame_b.fx
# frame_a.fy ~ -frame_b.fy
# end
# end
@mtkmodel RigidLink begin
@parameters begin
l = l
end
@components begin
frame_a = Frame()
frame_b = Frame()
end
@equations begin
# geometric
frame_b.x ~ frame_a.x + l*(1)
frame_b.y ~ frame_a.y + l*(frame_a.phi)
frame_b.phi ~ frame_a.phi
# torque
frame_a.tau ~ -(-frame_b.fx*l*(frame_a.phi) +
frame_b.fy*l*(1)) - frame_b.tau
# forces
frame_a.fx ~ -frame_b.fx
frame_a.fy ~ -frame_b.fy
end
end
@mtkmodel VariableForce begin
@components begin
fx = RealInput()
fy = RealInput()
frame = Frame()
end
@equations begin
frame.fx ~ -fx.u
frame.fy ~ -fy.u
frame.tau ~ 0.0
end
end
@mtkmodel BeamElement begin
@parameters begin
L = 50e-3, [description = "length of element"] # m
B = 0.8e-3, [description = "width of beam"] # m
H = 20.9e-3, [description = "height of beam"] # m
E = 200.0e9 , [description = "Modulus"]
ρ = 8000.0, [description = "Density"]
d = 0.01, [description = "damping"]
end
@components begin
spring = SpringDamperRotational(name=:spring,phi_rel0=0.0,c=E*(H*B^3/12)/L,d=d)
beam_1 = RigidLink(;name=:beam_1,l=L/2)
beam_2 = RigidLink(;name=:beam_2,l=L/2)
mass = MassInertia(;name=:mass,m=L*B*H*ρ,i=(L*B*H*ρ)*(L^2+B^2)/12)
end
@equations begin
connect(spring.frame_b,beam_1.frame_a)
connect(beam_1.frame_b,beam_2.frame_a,mass.frame) #mass.frame
end
end
#-------------------------------System definition------------------------------------------#
N=5
@mtkmodel system begin
@components begin
fixed = Fixed(;name=:fixed)
elements = [BeamElement(;name=Symbol(string("el_",i))) for i in 1:N]
force = VariableForce(;name=:force)
step = Step(;name=:step, height=1, offset=0, start_time=0.1, duration=Inf, smooth=true)
end
@equations begin
# for i=1:1 # use below syntax
# connect(elements[i].beam_2.frame_b,elements[i+1].spring.frame_a)
# end
[connect(elements[i].beam_2.frame_b,elements[i+1].spring.frame_a) for i in 1:(N-1)]
connect(fixed.frame,elements[1].spring.frame_a)
connect(elements[N].beam_2.frame_b,force.frame)
connect(step.output,force.fy)
force.fx.u ~ 0.0
end
end
@mtkbuild sys = system()
@time prob = ODEProblem(sys,[[eval(Meta.parse("sys.el_$i.spring.phi_rel => 0.0")) for i in 1:N]
[eval(Meta.parse("sys.el_$i.spring.w_rel => 0.0")) for i in 1:N];]
, (0, 1.0)) #[eval(Meta.parse("sys.el_$i.spring.w_rel => 0.0")) for i in 1:3]]
sol = solve(prob)
plot(sol,idxs=[eval(Meta.parse("sys.el_$i.beam_2.frame_b.y")) for i in 1:N])