I'm trying to write a program to calculate the steady state response of a beam structure to harmonic loading. The system is damped not only by 'structural damping' (which is linear) but also by air resistance (which is nonlinear, and calculated using $F_{li}=-\frac{1}{2}\rho A_i C_{Di}u_{li}|u_{li}|$... the minus sign and absolute value resulting in an air resistance force that acts to oppose the motion of the beam). Although the air resistance doesn't vary sinusoidally when the motion is sinusoidal, it isn't far off so I want to represent it as a harmonic force acting to oppose the loading force.
I have the mass and stiffness matrices for the beam ($\boldsymbol{[M]}$ and $\boldsymbol{[K]}$) and so I am able to obtain the eigenvectors and eigenvalues ($\boldsymbol{[V]}$ and $\boldsymbol{[D]}$). The eigenvalues yields a vector containing the natural frequencies of the system (the square root of the diagonal entries of $\boldsymbol{[D]}$, I'll call it $\boldsymbol{\omega}$.
The system is excited by a force vector in the global reference frame $\boldsymbol{F_0}$ at an excitation frequency $\omega_F$.
There is a rotation matrix which rotates the beams from the local to the global reference frame $\boldsymbol{[R]}$.
I have successfully used the procedure to calculate the steady state response without aerodynamic damping:
- Get eigenvalues and eigenvectors $\boldsymbol{[V]}$ and $\boldsymbol{[D]}$
- Get generalised force vector $\boldsymbol{Q}=\boldsymbol{[V]}^T \boldsymbol{F_0}$
- Calculate modal participation factor vector $\boldsymbol{q}$. The individual elements are calculated as follows:
$$q_i=\frac{Q_i}{{\omega_i}^2}\frac{1}{\sqrt{[(1-(\omega_F/\omega_i)^2)^2+(2\zeta_i\omega_F/\omega_i)^2]}}$$
We can but the elements of $\boldsymbol{q}$ without the scaling given by $\boldsymbol{Q}$ on the diagonal of a matrix $\boldsymbol{[A]}$, which I'll use later.
Calculate response from $\boldsymbol{x}=\boldsymbol{[V]q}$
The problem is that in the aerodynamic forcing equation above, $u_i$ (the nodal velocity in the local reference frame) is dependent on $\boldsymbol{x}$, which is calculated as $\boldsymbol{u_l}=\omega\boldsymbol{[R]^Tx}$ so the equation cannot be solved directly for $\boldsymbol{x}$.
I've tried solving it by using Newton-Raphson iterations but I've not had much luck... for some reason it works when the beams are lying in the original reference plane but as soon as $\boldsymbol{[R]}$ is anything other than identity it stops working. Does anyone have any pointers on how how to solve this:
$$\boldsymbol{x}=\boldsymbol{[V][A][V]^TF_0+[V][A][V]^TF_d}$$
where
$$\boldsymbol{F_d}=\boldsymbol{[R]F_l}$$
The components of $\boldsymbol{F_l}$ are calculated from $F_{li}=-\frac{1}{2}\rho A_i C_{Di}u_{li}|u_{li}|$ in which $\boldsymbol{u_l}=\omega\boldsymbol{[R]^Tx}$.
Any help will be very much appreciated... although even if I don't get any I've definitely improved my Tex skills!