Equações de Lagrange

T = (1/2)m (((θ'[t])^2)(d^2)+(Ω^2) (d^2)((Sin[θ[t]])^2)); 
V = (1/2)k (d^2)(((1-Cos[θ[t]]))^2)+ m g d (1-Cos[θ[t]]);
EL = D[D[T,θ'[t]],t]-D[T-V,θ[t]]-(-b θ'[t])//TrigReduce

EM = Solve[EL==0, θ''[t]]//First//Collect[#, {Sin[θ[t]], Sin[2θ[t]]}, FullSimplify]&

TimeScalesRule = {g->((ωp)^2) d, k-> ((ωn)^2) m};

EMS = EM/.TimeScalesRule //FullSimplify

Posição de equilíbrio dinâmico do sistema

Solve[-((ωn)^2)-((ωp)^2)+((Ω^2)+((ωn)^2)) Cos[θ[t]]==0, θ[t]]

Ruleθ = {θ[t]->ArcCos[((((ωn)^2)+((ωp)^2))/((Ω^2)+((ωn)^2)))]+ϵ δ[t]}

RuleSin ={Sin[θ[t]]->((Series[Sin[θ[t]]/.Ruleθ, {ϵ,0,1}]//Normal)/.{ϵ->1})}

RuleCos ={Cos[θ[t]]->((Series[Cos[θ[t]]/.Ruleθ, {ϵ,0,1}]//Normal)/.{ϵ->1})}

(EMS)/.Join[RuleSin,RuleCos]/.{θ''[t]-> δ''[t], θ'[t]-> δ'[t]}//FullSimplify//Collect[#,δ[t],FullSimplify]&

Simulação numérica

EME =Join[EM/.{Rule->Equal}, {θ[0]==θ0, θ'[0]==ω0}]

PAR = {m->1.0, g->9.81, k->100000, Ω->1000*2*π/60, d-> 0.10, θ0->π/6, ω0-> 0, b->0.001};

SOL = NDSolve[EME/.PAR, θ, {t,0,100}]

Plot[θ[t]/.SOL, {t,0,100}, PlotRange->All]