I face a problem where I know the beam's local stiffness matrix and I want to find the global stiffness matrix.
My problem involves the rotation of this matrix and I did not find any solution to this over the internet. Can anyone help me?
Question 1: How to find the expanded rotation matrix $\left[R_{exp}\right]$?
$$
\left[K_{glo}\right]_{12 \times 12} = \left[R_{exp}\right]_{12 \times 12} \cdot \left[K_{loc}\right]_{12 \times 12} \cdot \left[R_{exp}\right]_{12 \times 12}^{T}
$$
Question 2: If $\left[R_{exp}\right]$ has the same properties of the rotation matrix $\left[R\right]_{3 \times 3}$, we would have
$$\left[R_{exp}\right]^{T} = \left[R_{exp}\right]^{-1}$$
But if $\left[R_{exp}\right]$'s inverse it's not the same as it's transpose, which must I use?
$$
\left[K_{glo}\right] = \left[R_{exp}\right] \cdot \left[K_{loc}\right] \cdot \left[R_{exp}\right]^{T}
$$
or
$$
\left[K_{glo}\right] = \left[R_{exp}\right] \cdot \left[K_{loc}\right] \cdot \left[R_{exp}\right]^{-1}
$$
Problem description
Let's say that we have two system of coordinates:
- global: $(\vec{x}, \vec{y}, \ \vec{z})$
- local: $(\vec{n}, \ \vec{r}, \ \vec{t})$
I have a beam with the extremities at $\vec{0} = (0, 0, 0)$ and $\vec{p} = (p_x, p_y, p_z)$.
- Compute $R$
My idea was inspired by the rotations in 3D. My first objective
is to bring $\vec{p}$ to the $x$-axis. That means, get $R$ such that$$ \left(L, \ 0, \ 0\right) = L \cdot \vec{x} = \left[R\right] \cdot \vec{p} = \left[R\right] \cdot L \cdot \vec{n} $$
where $L = \|\vec{p}\| = \sqrt{p_x^2 + p_y^2 + p_z^2}$.
As the general equation for rotation, around the unit vector $\vec{u}$
by an amount of $\theta$ is given by$$ \left[R\right] = \left[I\right] \cdot \cos \theta + \left[u\right]_{\times} \cdot \sin \theta + \vec{u} \otimes \vec{u} \cdot \left(1-\cos \theta\right) $$
To bring $\vec{p}$ into $L \cdot \vec{x}$, we just need to define
$\vec{u}$ the vector perpendicular to $\vec{x}$ and $\vec{p}$$$ \vec{u} = \dfrac{\vec{p} \times \vec{x}}{L \cdot \sin \theta} = \dfrac{\left(0, \ p_z, -p_y\right)}{L \cdot \sin \theta} $$
And the rotation angle $\theta$ is the angle between $\vec{p}$ and
$\vec{x}$$$ \cos \theta = \dfrac{\vec{p} \cdot \vec{x}}{L} = \dfrac{p_{x}}{L} $$
Calculating $\left[R\right]$:
$$ \left[u\right]_{\times} =
\begin{bmatrix} & -u_{z} & u_{y} \\ u_{z} & & -u_{x} \\ -u_{y} & u_{x} & \end{bmatrix}
= \dfrac{1}{L \cdot \sin \theta} \begin{bmatrix} & p_{y} & p_{z} \\
-p_{y} & & \\
-p_{z} & & \end{bmatrix} $$$$ \left[R(\vec{p})\right] = \dfrac{1}{L} \begin{bmatrix} p_{x} &
p_{y} & p_{z} \\
-p_{y} & p_{x} & \\
-p_{z} & & p_{x} \end{bmatrix} + \dfrac{1}{p_y^2 + p_z^2} \left(1-\dfrac{p_{x}}{L}\right) \begin{bmatrix} & & \\ & p_{y}^2 &
-p_y p_z \\ & -p_y p_z & p_{z}^2 \end{bmatrix} $$
As we can see $R(\vec{p})^{T} \ne R(-\vec{p})$
- Local Stiffness Matrix:
The local stiffness matrix in the space $(\vec{n}, \ \vec{r}, \ \vec{t})$ is well-known and it's like
$$ \underbrace{\begin{bmatrix}
\square & & & & & & \square & & & & & \\
& \square & & & & \square & & \square & & & & \square \\
& & \square & & \square & & & & \square & & \square & \\
& & & \square & & & & & & \square & & \\
& & \square & & \square & & & & \square & & \square & \\
& \square & & & & \square & & \square & & & & \square \\
\square & & & & & & \square & & & & & \\
& \square & & & & \square & & \square & & & & \square \\
& & \square & & \square & & & & \square & & \square & \\
& & & \square & & & & & & \square & & \\
& & \square & & \square & & & & \square & & \square & \\
& \square & & & & \square & & \square & & & & \square \\
\end{bmatrix}}_{\left[K_{loc}\right]}
\underbrace{\begin{bmatrix}
u_{n1} \\ u_{r1} \\ u_{t1} \\ \theta_{n1} \\ \theta_{r1} \\ \theta_{t1} \\
u_{n2} \\ u_{r2} \\ u_{t2} \\ \theta_{n2} \\ \theta_{r2} \\ \theta_{t2}
\end{bmatrix}}_{\left[U_{loc}\right]} =
\underbrace{\begin{bmatrix}
F_{n1} \\ F_{r1} \\ F_{t1} \\ M_{n1} \\ M_{r1} \\ M_{t1} \\
F_{n2} \\ F_{r2} \\ F_{t2} \\ M_{n2} \\ M_{r2} \\ M_{t2}
\end{bmatrix}}_{\left[F_{loc}\right]}
$$
$u$ in this case is the displacement, $\theta_{ij}$ is the tangent's direction, $F$ is the force and $M$ the momentum.
- Attempt to compute $\left[K_{glo}\right]$
My attempt to compute $\left[K_{glo}\right]$ is mounting the matrix $\left[R_{exp}\right]$ with the matrix $\left[R\left(\vec{p}\right)\right]$ in the diagonal, as we can transform
$$ \begin{bmatrix}
u_{x1} \\ u_{y1} \\ u_{z1}
\end{bmatrix}
= \left[R\left(\vec{p}\right)\right] \cdot
\begin{bmatrix}
u_{n1} \\ u_{r1} \\ u_{t1}
\end{bmatrix} $$
$$
\begin{bmatrix}
\theta_{x1} \\ \theta_{y1} \\ \theta_{z1}
\end{bmatrix}
= \left[R\left(\vec{p}\right)\right] \cdot
\begin{bmatrix}
\theta_{n1} \\ \theta_{r1} \\ \theta_{t1}
\end{bmatrix} $$
$$ \begin{bmatrix}
u_{x1} \\ u_{y1} \\ u_{z1} \\ \theta_{x1} \\ \theta_{y1} \\ \theta_{z1} \\
u_{x2} \\ u_{y2} \\ u_{z2} \\ \theta_{x2} \\ \theta_{y2} \\ \theta_{z2}
\end{bmatrix} = \left[R_{exp}\right]
\begin{bmatrix}
u_{n1} \\ u_{r1} \\ u_{t1} \\ \theta_{n1} \\ \theta_{r1} \\ \theta_{t1} \\
u_{n2} \\ u_{r2} \\ u_{t2} \\ \theta_{n2} \\ \theta_{r2} \\ \theta_{t2}
\end{bmatrix}$$
$$ \left[R_{exp}\right] =
\begin{bmatrix}
\left[R(\vec{p})\right] & & & \\
& \left[R(\vec{p})\right] & & \\
& & \left[R(\vec{p})\right] & \\
& & & \left[R(\vec{p})\right] \\
\end{bmatrix}$$
Problem 1: When I test it with $\vec{p} = (p_x, \ p_y, \ 0)$ it works fine. $K_{glo}$ is correct. But for any $p_z \ne 0$ the results are wrong.
Problem 2: If $p_z = L \cdot \vec{z}$, the received matrix $\left[K_{glo}\right]$ it's the negative of it should be.
Problem 3: When I rotate $\vec{p}$ to $L \vec{x}$, it's not guaranteed that $\vec{r}$ becomes $\vec{y}$ and $\vec{t}$ becomes $\vec{z}$.
Does anyone know how to solve it?
PS: I'm testing with python with random values. As I don't know a formula to the global matrix, to verify if the global matrix is correct, I solve a system: Given the known values of $U_{glo}$ and $F_{glo}$:
- if $\left[K_{glo}\right] \cdot \left[U_{glo}\right] = \left[F_{glo}\right]$,
- $\left[K_{glo}\right]$ is correct
- else
- $\left[K_{glo}\right]$ is not correct
Best Answer
After a while, I found the solution. I will divide it into two parts: