FEM, CFD, CAE, and all that.

## Tag: mathematica

### Harmonic vibration of a damped system

Unlike in previous post that phase change could only be $0$ or $\pi$, in the damped system, as the ratio between dynamics response and static response is expressed in the same way as the undamped case
$\frac{u_d}{u_{s0}}=R_d\sin(\omega t-\phi),$
the phase lag $\phi$ varies in $[0,\pi]$. The following Mathematica outputs the displacement response factor, phase lag and normalized time history as the function of damping ratio $\zeta$ and frequency ratio $w=\omega/\omega_0$.

Rd[w_, \[Zeta]_] := 1/Sqrt[(1 - w^2)^2 + (2 \[Zeta] w)^2]
\[Phi][w_, \[Zeta]_] := -ArcTan[1 - w^2, -2 \[Zeta] w]
ud[w_, \[Zeta]_, t_] :=
Rd[w, \[Zeta]] Sin[2 Pi (t - \[Phi][w, \[Zeta]]/(2 Pi))]
us[t_] := Sin[2 Pi t]
Manipulate[
GraphicsRow[{GraphicsColumn[{Plot[Rd[w, \[Zeta]], {w, 0, 3},
AxesOrigin -> {1, 1}, PlotRange -> {{0, 3}, {0, 5}},
AxesLabel -> {"\!$$\*FractionBox[\"\[Omega]\", SubscriptBox[\"\ \[Omega]\", \"0\"]]$$", "\!$$\*SubscriptBox[\"R\", \"d\"]$$"},
Epilog -> {PointSize[Large], Red, Point[{w, Rd[w, \[Zeta]]}]},
Ticks -> {{0, 1, 2, 3}, {0, 1, 2, 3, 4, 5}}],
Plot[\[Phi][w, \[Zeta]], {w, 0, 3}, AxesOrigin -> {1, Pi/2},
PlotRange -> {{0, 3}, {0, Pi}},
AxesLabel -> {"\!$$\*FractionBox[\"\[Omega]\", SubscriptBox[\"\ \[Omega]\", \"0\"]]$$", "\[Phi]"},
Epilog -> {PointSize[Large], Red,
Point[{w, \[Phi][w, \[Zeta]]}]},
Ticks -> {{0, 1, 2, 3}, {0, Pi/2, Pi}}]}, Frame -> All],
Plot[{ud[w, \[Zeta], t], us[t]}, {t, 0, 2},
PlotStyle -> {Thick, Dashed},
AxesLabel -> {"t/T",
"\!$$\*SubscriptBox[\"u\", \"d\"]$$/\!$$\*SubscriptBox[\"u\", \ \"st\"]$$"}]}, Frame -> All], {{w, 0.5, "Frequency ratio"}, 0,
3}, {{\[Zeta], 0.2, "Damping ratio"}, 0, 1}]

Special cases include

1. Long period excitation, i.e., $\omega\ll\omega_0$, gives a pseudo-static response. In this case, the system “waits” until it “feels” the excitation completely. $R_d$ is greater but very close to $1$, and the displacement is essentially in phase with excitation force, in other words, dynamic effect is near to none.

2. Short period excitation gives very small $R_d$, though the phase lag is $\pi$. Here the system barely reacts when the load is reversed, thus leads to small displacement.

3. Resonant period, i.e., $\omega\approx\omega_0$, leads to $\phi\approx\pi/2$. Now $R_d$ is very sensitive to damping change, namely, the response is controlled by the damping: a small change of damping ratio $\zeta$ leads to great reaction of the structure. When $\omega=\omega_0$, we have $R_d=1/(2 \zeta)$.

Resonant frequency: small damping

Resonant frequency: medium damping

Resonant frequency: large damping

The last case is the what’s essentially behind the viscous damping devices applied to buildings like this.

### Harmonic vibration of an undamped system

This can be found in virtually every dynamics textbook. The undamped system follows the frequency of the excitation, while the displacement response factor $R_d$.  Response is out of phase when excitation frequency is greater than the natural frequency:  $\omega>\omega_0$. If there is only one thing to keep in mind of harmonic SDF vibration,  I guess it’s the plot for the response against the frequency ratio. The dynamic object in Mathematica can be built by:

Manipulate[GraphicsRow[
{GraphicsColumn[{Plot[Rd[r], {r, 0, 3}, AxesOrigin -> {1, 0},
Frame -> True, PlotRange -> {{0, 3}, {-5, 5}},
AxesLabel -> {"\!$$\*FractionBox[\"\[Omega]\", SubscriptBox[\"\ \[Omega]\", \"0\"]]$$", "\!$$\*SubscriptBox[\"R\", \"d\"]$$"},
Epilog -> {Red, PointSize[Large], Point[{w, Rd[w]}]}],
Plot[Abs[Rd[r]], {r, 0, 3}, AxesOrigin -> {1, 1}, Frame -> True,
PlotRange -> {{0, 3}, {0, 5}},
AxesLabel -> {"\!$$\*FractionBox[\"\[Omega]\", SubscriptBox[\"\ \[Omega]\", \"0\"]]$$", "|\!$$\*SubscriptBox[\"R\", \"d\"]$$|"},
Epilog -> {Red, PointSize[Large], Point[{w, Abs[Rd[w]]}]}]},
Alignment -> Right, Frame -> All],
Plot[{1/(1 - w^2) Sin[w t], Sin[w t]}, {t, 0, 4 Pi}, Frame -> True,
PlotStyle -> {Thick, Dashed},
PlotRange -> {{0, 4 Pi}, {-10, 10}}, AxesLabel -> {"t", "u"}]
}, Frame -> All]
, {{w, 0, "Frequency ratio"}, 0, 3}]

In the right, the dashed curve is the excitation history, while the thick one is for response, which would be out of phase to the excitation when frequency ratio is greater than 1.0 ($\omega/\omega_0>1.0$).

### Behind Watson and Alpha

A diagram from Wolfram blog explains what’s behind Watson and Alpha

### Cantor set by Mathematica

Essentially one can not really show the complete Cantor set, here by its definition is what’s left after first n remove of “one-third” open sets. The Mathematica script is then:

Cantor[x__, n_] :=
Nest[Join[Partition[#[[All, 1]], 1],
Partition[#[[All, 1]] + (#[[All, 2]] - #[[All, 1]])/3, 1],
2] \[Union]
Join[Partition[#[[All, 1]] + 2 (#[[All, 2]] - #[[All, 1]])/3, 1],
Partition[#[[All, 2]], 1], 2] &, x, n]

Now

Cantor[{{0, 1}}, 2]

gives

{{0, 1/9}, {2/9, 1/3}, {2/3, 7/9}, {8/9, 1}}

where the first argument is the interval to be sliced, and the second argument gives the level (iteration times) of generation. The output is the intervals left after the remove action. Try

Cantor[{{0, 1}}, 3]

we have

{{0, 1/27}, {2/27, 1/9}, {2/9, 7/27}, {8/27, 1/3}, {2/3, 19/27}, {20/
27, 7/9}, {8/9, 25/27}, {26/27, 1}}

and

Cantor[{{0, 1}}, 4]

{{0, 1/81}, {2/81, 1/27}, {2/27, 7/81}, {8/81, 1/9}, {2/9, 19/
81}, {20/81, 7/27}, {8/27, 25/81}, {26/81, 1/3}, {2/3, 55/81}, {56/
81, 19/27}, {20/27, 61/81}, {62/81, 7/9}, {8/9, 73/81}, {74/81, 25/
27}, {26/27, 79/81}, {80/81, 1}}

and

Cantor[{{0, 1}}, 5]

{{0, 1/243}, {2/243, 1/81}, {2/81, 7/243}, {8/243, 1/27}, {2/27, 19/
243}, {20/243, 7/81}, {8/81, 25/243}, {26/243, 1/9}, {2/9, 55/
243}, {56/243, 19/81}, {20/81, 61/243}, {62/243, 7/27}, {8/27, 73/
243}, {74/243, 25/81}, {26/81, 79/243}, {80/243, 1/3}, {2/3, 163/
243}, {164/243, 55/81}, {56/81, 169/243}, {170/243, 19/27}, {20/27,
181/243}, {182/243, 61/81}, {62/81, 187/243}, {188/243, 7/9}, {8/9,
217/243}, {218/243, 73/81}, {74/81, 223/243}, {224/243, 25/27}, {26/
27, 235/243}, {236/243, 79/81}, {80/81, 241/243}, {242/243, 1}}

Plotting could be done using

PlotCantor[{x1_, x2_}] :=
Plot[x = 1, {x, x1, x2}, Filling -> Bottom,
PlotRange -> {{0, 1}, {0, 1}}, Axes -> {True, False}]
Show[PlotCantor /@ Cantor[{{0, 1}}, 2]]
Show[PlotCantor /@ Cantor[{{0, 1}}, 4]]

### Linear shape function, mass matrix, and element volume

In unsteady CFD problem, under certain time integration scheme, L2 inner product of shape functions gives what so called mass matrix, i.e. the matrix with its entries in the form $M_{ij}=\int_{\Omega}{N_i N_j}\,d{\Omega}$, where $N_i$ is shape function. When shape function is linear, there is a better way to demonstrate $M_{ij}$ is proportional to local element volume (area) than trivially plugging in and calculating or using Gaussian integral. Assume in 2D triangular element ABC whose area is $A^{ABC}$,  and $N_i(x,y), i=1,2,3$ stands for the shape function whose unit value is reached at A, B, C correspondingly.

If we interpret local mass entry $M_{13}^{ABC}=(N_1, N_3)$ over element ABC shown above as the mass of a piece of wood splinter in the shape of tetrahedron A1-A-B-C, where $N_1(x,y)$ is the thickness function on the element, while $N_3(x,y)$ is the density function, then it’s not difficult to see that $M_{13}=(N_1,N_3)=(N_2,N_3)=M_{23}$, because they are both the mass of tetrahedron with the same volume and density distribution from bottom (A1-A-B or B1-A-B) to top (C). The next step is to show that the mass of tetrahedron C1-A-B-C ($M_{33}$) is the half of the total mass of the prism A1-B1-C1-A-B-C with the same density distribution. The argument is based on linearity, which allows us to assign a single parameter for above integrals. Set this single coordinate $s$ as the line parameter from C to line A-B, i.e. $s=0$ on AB, $s=1$ at C. Then, for mass of prism and tetrahedron, we both have

mass $=\int{L(s) H(s) D(s)}\,ds$

where $L\times H$ gives the area for cross section parallel to A1-A-B-B1 at $s$ . $H$ is the thickness of ABC, while $D$ is the density distribution. We have

$L\propto (1-s),\quad D\propto s$,
and $H=1$ for prism, $H\propto s$ for tetrahedron. Then we have

$M_p\propto \int{(1-s) s}\,ds$

$M_{33}\propto \int{(1-s) s^2}\,ds$

which proves

$2 M_{33}=M_p=\frac{1}{3} A_{ABC}$

Since $M_{33}+M_{13}+M_{23}$ is nothing but the sum of the mass of three tetrahedrons, which is the mass of the prism. This sum is then $\frac{1}{3}A_{ABC}$. Using $M_{13}=M_{23}$, finally we have

$M_{11}=M_{22}=M_{33}=\frac{1}{6}A_{ABC}$

$M_{12}=M_{13}=M_{23}=\frac{1}{12}A_{ABC}$

Note: for 3D tetrahedron element, linear shape function gives similar relation, in this case, above ratio over the volume of a tetrahedron would be 1/10 and 1/20:

$M_{11}=M_{22}=M_{33}=M_{44}=\frac{1}{10}V_{ABCD}$

$M_{12}=M_{13}=M_{14}=M_{23}=M_{24}=M_{34}=\frac{1}{20}V_{ABCD}$