File size: 2,117 Bytes
7c0231e
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
"""Grounded linear viscous damper, independent SciPy time-domain reference.

Apache-2.0. Uses the exported solid-component ROM, not a neural network.
"""
import numpy as np
from scipy.linalg import expm
from assembly import solve


def free_decay(chain, rolls=None, supports=None, node=None, axis=2, c=10000.,
               amplitude=.001, mode_count=12, frames=600, duration=None):
    """First-mode release; peak coordinate amplitude in metres, c in N s/m.

    Zero initial velocity, no background damping or imposed forcing. Motion is
    a linear perturbation about equilibrium; thermal prestress is not included.
    Returns exact-in-time samples for the selected truncated modal equations.
    """
    node=len(chain) if node is None else node
    if axis not in (0,1,2) or not 0<=node<=len(chain) or c<0 or not np.isfinite(c):
        raise ValueError('Invalid grounded damper')
    if not 0<amplitude<=.01 or mode_count<1 or frames<1:raise ValueError('Invalid release settings')
    r=solve(chain,rolls=rolls,supports=supports,nmodes=mode_count)
    w=2*np.pi*r['frequency'];n=len(w)
    locations=[(abs(m[j,k,0]),a,j,k) for a,m in enumerate(r['mode_displacements']) for j in range(len(m)) for k in range(3)]
    _,a,j,k=max(locations);probe=r['mode_displacements'][a][j,k]
    z=np.zeros(n);z[0]=amplitude/probe[0]
    b=r['modes'][node*6+axis];C=c*np.outer(b,b)
    A=np.block([[np.zeros((n,n)),np.eye(n)],[-np.diag(w*w),-C]])
    T=6/r['frequency'][0] if duration is None else duration
    if T<=0:raise ValueError('Invalid duration')
    time=np.linspace(0,T,frames+1);P=expm(A*T/frames);state=np.r_[z,np.zeros(n)];states=[]
    for _ in time:states.append(state.copy());state=P@state
    states=np.array(states);q,v=states[:,:n],states[:,n:]
    energy=.5*np.sum(v*v+(q*w)**2,axis=1)
    return dict(time=time,signal=q@probe,reference=amplitude*np.cos(w[0]*time),
                force=-c*(v@b),energy=energy,coordinates=q,velocities=v,
                probe=dict(component=a,sample=j,axis=k),frequency=r['frequency'],
                method='scipy matrix exponential of the modal state-space system')