分数质量-弹簧-阻尼系统研究（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客  

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

近年来，国内外专家开始将环形弹簧引入到结构减震领域，对其开展了初步的研究。Filiatrault 等设计了一种含有常规环形弹簧的减震装置，并基于振动台试验考察了其减震性能。Khoo 等研制了一种钢框架基于环形弹簧的钢框架结构耗能节点。Fang 等提出了一种含有超弹性环形弹簧的钢框架结构自复位耗能型梁柱节点。Issa 等研发了一种活塞型环形弹簧自复位耗能支撑。以上研究表明，环形弹簧是一种兼具耗能和自复位特性的减震部件，可用于新型自复位控制系统的研制。

尽管国内外学者对基于环形弹簧的建筑结构减震控制开展了初步的研究，但是，环形弹簧减震器的研究基础仍然十分薄弱，相关的理论和试验研究偏少。本文研究分数质量-弹簧-阻尼器系统（FMSDS）数学模型的数值解。

📚2 运行结果

 

 

 

部分代码：

function[t,x,v]= FMSDS(m,b,k,alpha,beta,tau,dtau,x0,v0,f)
%
% DESCRIPTION:
%
%  Function returns numerical solution of the 
%  Fraction Mass-Spring-Damper System (FMSDS)
%  using Grunwald-Letnikov derivative definition.
%
%                    |
%  /| Damper      --------
%  /|   --       |        | 
%  /|--- |-------| Mass   |  Force
%  /|   -- b     |      m |<--------> f
%  /| Spring     |        |
%  /|-/\/\/\/\---|        |
%  /|       k     --------  
%  /|                | Displacement,  x 
%  /|<------------------------------>
%  /|                | Velocity, v 
%                    0
%
%  Fractional differential equation:
%
%      d^alpha             d^beta
%   m ---------- x(t) + b -------- x(t) + k x(t) = f     for 0 <= t <= tau 
%      d t^alpha           d^beta                            alpha > beta
%
%  Initial conditions:
%                       x(0) = x0                          
%               x'(0) = v(0) = v0                          
%
%  Grunwald-Letnikov definition:
%
%                    Nf                                
%                   sum b_j.y(t-j.dtau)               
%  d^alpha          j=0                
%  --------- y(t) = -------------------
%  d t^alpha        dtau^alpha         
%   
%  Difference equation :
%
%     p                   p
%    sum bc_j.x(p-j)     sum cc_j.x(p-j)  
%    j=0                 j=0    
%  m --------------- + b --------------- + k.x(p) = f
%    dtau^alpha          dtau^beta
%
%  or
%                                  p   
%   m.dtau^(-alpha) ( bc_0.x(p) + sum bc_j.x(p-j) ) +
%                                 j=1  
%
%                                  p   
% + b.dtau^(-beta)  ( cc_0.x(p) + sum cc_j.x(p-j) ) + k.x(p) = f
%                                 j=1  
%
%  or
%                           p                              p
%  x(p)=(f-m.dtau^(-alpha) sum bc_j.x(p-j)-b.dtau^(-beta) sum cc_j.x(p-j))/
%                          j=1                            j=1
%      
%       (m.dtau^(-alpha).bc_0 + b.dtau^(-beta).cc_0 + k)
%      
% INPUTS:
%  m       - mass                                [ kg     ]
%  b       - damper coefficient                  [ kg/s   ]
%  k       - stiffness (or spring constant)      [ kg/s^2 ]
%  alpha   - fractional-order derivative         [ -      ]
%  beta    - fractional-order derivative         [ -      ]
%  tau     - simulation time                     [ s      ]
%  dtau    - time step                           [ s      ]
%  x0      - initial displacement                [ m      ]
%  v0      - initial velocity                    [ m/s    ]
%  f       - force                               [ N      ]
%
% OUTPUTS:
%  t       - vector of time                      [ s      ]
%  x       - vector of distance                  [ m      ]
%  v       - vector of velocity                  [ m/s    ]
%
% AUXILIARY VARIABLES:
%  n       - length of time vector
%  p, i, j - index
%  bc, cc  - binomial coefficients
%  s1, s2  - sum

%
  n = tau/dtau+1;

  t(1) = 0.0;
  x(1) = x0; 
  v(1) = v0; 

  t(2) = dtau;
  x(2) = x(1) + dtau*v(1); 
  v(2) = ( x(2) - x(1) )/dtau; 

  p = 1;
  for i = 3:n
    p = p + 1;  
    t(i) = (i-1)*dtau;  
    
    bc = zeros(1,p+1);
    cc = zeros(1,p+1);
    bc(1) = 1;
    cc(1) = 1;
    for j = 2:p+1
        bc(j) = (1 - (1+alpha)/(j-1))*bc(j-1);
        cc(j) = (1 - (1+beta)/(j-1))*cc(j-1);
    end
       
    s1 = 0;
    s2 = 0;
    for j = 2:p+1
        s1 = s1 + bc(j)*x(p+2-j);
        s2 = s2 + cc(j)*x(p+2-j);
    end
        
    x(i) = ( f - m*dtau^(-alpha)*s1    - b*dtau^(-beta)*s2       )/...
           (     m*dtau^(-alpha)*bc(1) + b*dtau^(-beta)*cc(1) + k); 
    v(i) = ( x(i)  -  x(i-1) )/dtau; 
  end
end

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1]庄鹏,徐蒙,王尉.环形弹簧阻尼器的力学性能研究[J].建筑技术开发,2020,47(10):101-103.

[2]徐蒙,庄鹏,韩淼.环形弹簧阻尼器的滞回性能试验研究[J].动力学与控制学报,2020,18(05):79-85.

🌈4 Matlab代码实现