Coupling models together using the Norton-Thevenin method ========================================================= Using the ``frclim.ntfl`` routine to compute coupled system response -------------------------------------------------------------------- This and other jupyter notebooks are available here: https://github.com/twmacro/pyyeti/tree/master/docs/tutorials. The goal of this tutorial is to show that the Norton-Thevenin (NT) method of coupling a “Source” model and a “Load” model together and computing the system response can be equivalent to coupling the models together the “old fashioned” way and computing the system response. The old fashioned way is what is done in a standard Nastran superelement run; that is, it forms coupled system mass, damping, and stiffness matrices by enforcing boundary compatibility, computes system modes (typically), and solves the equations of motion. A typical Coupled Loads Analysis (CLA) couples models together in this fashion. The NT approach also enforces boundary compatibility and computes coupled system response, but does not form system matrices or compute system modes. Instead, it uses the “apparent masses” of the Source and the Load and the “free acceleration” of the Source to compute the system level responses. The NT coupling equations and their derivation are shown here: `pyyeti.frclim.ntfl <../modules/generated/pyyeti.frclim.ntfl.html#pyyeti.frclim.ntfl>`__ (see equations 4 & 5). Note that the NT method is naturally a frequency domain method and, for simplicity, this tutorial will compute only frequency domain responses. Also, these models have a **statically-indeterminate** interface, showing that the method is general. To demonstrate that the NT approach can equal CLA approach for a frequency domain analysis, we’ll pick two models already available in the pyYeti package. The “inboard” model will be the Source, and the “outboard” model will be the Load. These happen to both be Hurty-Craig-Bampton (HCB) models, but that is not necessary for applying the NT method. *Note:* No matter which format the models are in, be sure to include the complete flexibility relative to the boundary in both the Source and Load models when applying the NT method. For free-free models, including the residual vectors relative to the boundary degrees of freedom (DOF) (via RVDOF in Nastran, for example) is a great way to include this flexibility. HCB models will normally contain this required flexibility automatically in the normal and constraint modes. While that is likey the case with our models, for the demonstration here, we will simply consider them to be “exact” and will make no effort to ensure they accurately represent the original FEM models. Steps: 1. Load models 2. Set up a frequency domain forcing function 3. Couple the models together the old fashioned way 1. Compute system responses the old fashioned way 4. Couple the models together via the NT approach 1. Compute apparent masses and free acceleration 2. Compute system responses via the NT method 5. Compare results For this example, we’ll use a simple space-station model that has been reduced into two separate Hurty-Craig-Bampton models: an “inboard” model and an “outboard” model. These model attach at nodes 3, 11, 19, and 27. Notes: - These model use units of kg, mm, s - They are very light-weight trusses: - Mass of inboard model = 1.755 kg - Mass of outboard model = 1.590 kg The inboard model (the Source): .. image:: inboard.png The outboard model (the Load): .. image:: outboard.png -------------- First, do some imports: .. code:: ipython3 import inspect from pathlib import Path from io import StringIO import numpy as np from scipy import linalg as la import matplotlib.pyplot as plt import pandas as pd from pyyeti import cb, frclim, ode, ytools from pyyeti.nastran import bulk, op4, n2p from pyyeti.pp import PP Some settings specifically for the jupyter notebook: .. code:: ipython3 %matplotlib inline plt.rcParams['figure.figsize'] = [6.4, 4.8] plt.rcParams['figure.dpi'] = 150. Need path to data files: .. code:: ipython3 srcdir = Path(inspect.getfile(frclim)).parent / "tests" / "nas2cam_extseout" Load model data ^^^^^^^^^^^^^^^ Also add 2% critical modal damping for both models: .. code:: ipython3 zeta = 0.02 ids = ("out", "in") uset, cords, b = {}, {}, {} mats = {} for id in ids: uset[id], cords[id], b[id] = bulk.asm2uset(srcdir / f"{id}board.asm") mats[id] = op4.read(srcdir / f"{id}board.op4") # add damping: bxx = 0 * mats[id]["kxx"] q = ~b[id] lam = np.diag(mats[id]["kxx"])[q] damp = 2 * np.sqrt(lam) * zeta bxx[q, q] = damp mats[id]["bxx"] = bxx PP(mats, 2); .. code-block:: none [n=2] 'out': [n=30] 'kxx' : float64 ndarray 2116 elems: (46, 46) 'mxx' : float64 ndarray 2116 elems: (46, 46) 'bxx1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'k4xx1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'px1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'gpxx' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'gdxx' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'rvax' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'va' : float64 ndarray 46 elems: (46, 1) 'mug1' : float64 ndarray 1932 elems: (42, 46) 'mug1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mes1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mes1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mee1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mee1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfm' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfb' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfk' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfo' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mef1' : float64 ndarray 1104 elems: (24, 46) 'mef1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqgm' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqgb' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqgk' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqg1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmgm' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmgb' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmgk' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmg1o': float64 ndarray 1 elems: (1, 1) [[ 0.]] 'bxx' : float64 ndarray 2116 elems: (46, 46) 'in' : [n=30] 'kxx' : float64 ndarray 1024 elems: (32, 32) 'mxx' : float64 ndarray 1024 elems: (32, 32) 'bxx1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'k4xx1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'px' : float64 ndarray 960 elems: (32, 30) 'gpxx' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'gdxx' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'rvax' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'va' : float64 ndarray 32 elems: (32, 1) 'mug1' : float64 ndarray 1152 elems: (36, 32) 'mug1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mes1' : float64 ndarray 864 elems: (27, 32) 'mes1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mee1' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mee1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfm' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfb' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfk' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mgpfo' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mef1' : float64 ndarray 512 elems: (16, 32) 'mef1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqgm' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqgb' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqgk' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqg1o' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmgm' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmgb' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmgk' : float64 ndarray 1 elems: (1, 1) [[ 0.]] 'mqmg1o': float64 ndarray 1 elems: (1, 1) [[ 0.]] 'bxx' : float64 ndarray 1024 elems: (32, 32) For convenience and ease of reading, create some shorter names for the mass, damping, and stiffness matrices: .. code:: ipython3 maa = { "in": mats["in"]["mxx"], "out": mats["out"]["mxx"], } kaa = { "in": mats["in"]["kxx"], "out": mats["out"]["kxx"], } baa = { "in": mats["in"]["bxx"], "out": mats["out"]["bxx"], } Model checks ^^^^^^^^^^^^ We’ll use the `pyyeti.cb.cbcheck <../modules/generated/pyyeti.cb.cbcheck.html#pyyeti.cb.cbcheck>`__ routine to run model checks on both models. (There is also a tutorial on running model checks `here `__). For this notebook, the primary reason for running model checks was simply to get the modal effective mass tables since that data may be helpful in interpreting some of the plots below. However, `pyyeti.cb.cbcheck <../modules/generated/pyyeti.cb.cbcheck.html#pyyeti.cb.cbcheck>`__ also turned out to be helpful in diagnosing an numerical issue with the inboard model; more on that in a moment. First, run the model checks and display the modal effective mass tables: .. code:: ipython3 bset_in = b["in"].nonzero()[0] bset_out = b["out"].nonzero()[0] with StringIO() as f: chk_in = cb.cbcheck( f, # "cbcheck_in.txt", maa["in"], kaa["in"], bset_in, bset_in[:6], uset["in"], uref=[600, 150, 150], # this is the center if the 4 i/f nodes rb_norm=True, ) str_chk_in = f.getvalue() with StringIO() as f: chk_out = cb.cbcheck( f, # "cbcheck_out.txt", maa["out"], kaa["out"], bset_out, bset_out[:6], uset["out"], uref=[600, 150, 150], rb_norm=True, ) str_chk_out = f.getvalue() .. code-block:: none Iteration 1 completed Convergence: 4 of 6, tolerance range after 2 iterations is [1.196996633817e-09, 1.190690665692671e-05] Convergence: 6 of 6, tolerance range after 3 iterations is [3.4365399102737635e-12, 2.7573499099751664e-11] Show the modal effective mass tables: .. code:: ipython3 pd.options.display.float_format = '{:,.2f}'.format chk_in.effmass_percent .. raw:: html
T1 T2 T3 R1 R2 R3
Frq (Hz)
6.13 0.00 31.93 6.06 0.00 15.61 82.15
6.13 0.00 6.06 31.93 0.00 82.14 15.60
23.63 0.00 0.00 0.00 80.42 0.00 0.00
70.48 0.00 31.05 0.56 0.00 0.03 1.67
70.79 0.00 0.56 31.25 0.00 1.68 0.03
104.67 44.16 0.00 0.00 0.00 0.00 0.00
188.04 0.63 0.08 0.00 0.00 0.00 0.01
208.60 0.00 0.38 0.00 0.00 0.00 0.01
.. code:: ipython3 chk_out.effmass_percent .. raw:: html
T1 T2 T3 R1 R2 R3
Frq (Hz)
1.65 0.00 0.00 0.00 0.00 0.00 46.18
1.65 11.86 0.00 0.00 0.00 0.00 0.00
1.67 0.00 0.00 0.00 83.95 0.00 0.00
1.67 0.00 0.00 12.24 0.00 12.34 0.00
7.03 0.00 0.00 0.00 0.00 0.88 0.00
7.03 0.00 0.00 0.00 0.00 0.00 0.00
10.91 0.00 0.00 0.00 0.00 0.00 1.52
10.91 2.47 0.00 0.00 0.00 0.00 0.00
13.99 0.00 0.00 0.00 2.69 0.00 0.00
13.99 0.00 0.00 2.76 0.00 2.78 0.00
25.13 0.00 0.00 0.00 0.00 0.00 0.42
25.14 0.99 0.00 0.00 0.00 0.00 0.00
42.17 0.00 0.00 0.00 0.51 0.00 0.00
42.19 0.00 0.00 1.01 0.00 1.02 0.00
46.72 0.00 5.38 0.00 0.00 0.00 2.75
46.89 0.00 0.00 0.00 0.00 0.00 0.00
69.17 0.00 0.12 0.00 0.00 0.00 0.98
87.50 5.08 0.00 0.00 0.00 0.00 0.00
110.85 0.00 81.36 0.44 0.00 0.32 45.80
139.24 0.00 0.69 70.30 0.00 76.14 0.16
152.23 1.14 0.00 0.00 8.66 0.00 0.00
472.62 0.00 0.01 0.67 0.00 2.50 0.03
A quick look at the first few lines of the output: .. code:: ipython3 print(str_chk_in[:str_chk_in.index("\n\n")]) .. code-block:: none Mass matrix is symmetric. Warning: mass matrix is not positive-definite. However, subspace iteration succeeded, which means the mass must be close to positive-definite. Warning: stiffness matrix is not symmetric: abs-max-err = 5.96e-07 .. code:: ipython3 print(str_chk_out[:str_chk_out.index("\n\n")]) .. code-block:: none Mass matrix is symmetric. Mass matrix is positive-definite. Stiffness matrix is symmetric. The model checks showed that the inboard model mass and stiffness were not perfect: the mass was not positive-definite (but close), and the stiffness was not symmetric (but close). However, since they were close, all the model checks ran fine and show that the model is good. Why did the matrices come out this way? That’s how Nastran wrote them and they are fine as far as Nastran is concerned. These Python routines however look for a higher degree of precision. (It’s possible that converting units from millimeters to meters may avoid this issue by providing better relative scaling but, at the time of this writing, that has not been tried.) Because of the numerical imperfections, we’ll get a warning below from the `pyyeti.frclim.calcAM <../modules/generated/pyyeti.frclim.calcAM.html#pyyeti.frclim.calcAM>`__ routine about having trouble using the `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ solver to compute the apparent mass. After issuing the warning, it will automatically switch to the `pyyeti.ode.FreqDirect <../modules/generated/pyyeti.ode.FreqDirect.html#pyyeti.ode.FreqDirect>`__ solver to successfully compute the apparent mass. The reason that `pyyeti.cb.cbcheck <../modules/generated/pyyeti.cb.cbcheck.html#pyyeti.cb.cbcheck>`__ succeeded in solving the eigenproblem while `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ did not is simply because they use different eigen solvers: `pyyeti.cb.cbcheck <../modules/generated/pyyeti.cb.cbcheck.html#pyyeti.cb.cbcheck>`__ uses `pyyeti.ytools.eig_si <../modules/generated/pyyeti.ytools.eig_si.html#pyyeti.ytools.eig_si>`__, which has some leniency built-in, while `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ uses `scipy.linalg.eigh `__, which is very strict. Since it’s illuminating, we’ll show how close the mass matrix is to positive-definite by making a tiny modification to it, similar to what `pyyeti.ytools.eig_si <../modules/generated/pyyeti.ytools.eig_si.html#pyyeti.ytools.eig_si>`__ does. We’ll increase the diagonal values by a tiny amount: .. code:: ipython3 m = maa['in'].copy() print(f"Before: {ytools.mattype(m, 'posdef') = }") d = np.arange(m.shape[0]) m[d, d] *= 1.000000000000001 print(f"After: {ytools.mattype(m, 'posdef') = }") .. code-block:: none Before: ytools.mattype(m, 'posdef') = False After: ytools.mattype(m, 'posdef') = True Define the forcing function ^^^^^^^^^^^^^^^^^^^^^^^^^^^ We’ll use the first TLOAD vector of “inboard” to apply a force to the coupled system. This load vector applies forces to nodes 13 & 21. For simplicity, we’ll keep the forces constant across the frequency range. The analysis will cover the frequency range 0.01 to 100.0 Hz on a logarithmic scale. .. code:: ipython3 freq = np.geomspace(0.01, 100.0, 4 * 167 + 1) # use first TLOAD vector and expand it: force_in = mats["in"]["px"][:, :1] @ np.ones((1, len(freq))) force_in.shape .. code-block:: none (32, 669) Take a quick look at the inboard applied forces: .. code:: ipython3 np.set_printoptions(precision=3, linewidth=130, suppress=True) force_in[:, :5] .. code-block:: none array([[ 11.615, 11.615, 11.615, 11.615, 11.615], [ 0.158, 0.158, 0.158, 0.158, 0.158], [ -6.209, -6.209, -6.209, -6.209, -6.209], [-334.015, -334.015, -334.015, -334.015, -334.015], [-680.302, -680.302, -680.302, -680.302, -680.302], [-417.33 , -417.33 , -417.33 , -417.33 , -417.33 ], [ -0.204, -0.204, -0.204, -0.204, -0.204], [ 6.209, 6.209, 6.209, 6.209, 6.209], [ -9.713, -9.713, -9.713, -9.713, -9.713], [ 619.975, 619.975, 619.975, 619.975, 619.975], [-291.083, -291.083, -291.083, -291.083, -291.083], [-359.031, -359.031, -359.031, -359.031, -359.031], [ -11.615, -11.615, -11.615, -11.615, -11.615], [ 4.454, 4.454, 4.454, 4.454, 4.454], [ -1.912, -1.912, -1.912, -1.912, -1.912], [ -92.337, -92.337, -92.337, -92.337, -92.337], [-739.591, -739.591, -739.591, -739.591, -739.591], [-358.041, -358.041, -358.041, -358.041, -358.041], [ 9.713, 9.713, 9.713, 9.713, 9.713], [ 4.092, 4.092, 4.092, 4.092, 4.092], [ 1.912, 1.912, 1.912, 1.912, 1.912], [ -67.321, -67.321, -67.321, -67.321, -67.321], [ 679.264, 679.264, 679.264, 679.264, 679.264], [-231.794, -231.794, -231.794, -231.794, -231.794], [ -5.2 , -5.2 , -5.2 , -5.2 , -5.2 ], [ 1.575, 1.575, 1.575, 1.575, 1.575], [ -10.157, -10.157, -10.157, -10.157, -10.157], [ 5.768, 5.768, 5.768, 5.768, 5.768], [ -0.037, -0.037, -0.037, -0.037, -0.037], [ 0.603, 0.603, 0.603, 0.603, 0.603], [ -7.554, -7.554, -7.554, -7.554, -7.554], [ -5.275, -5.275, -5.275, -5.275, -5.275]]) Couple models together the old fashioned way: ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We’ll define a matrix ‘S’ that will enforce boundary compatibility. It defines the relationship from the DOF chosen to be the independent set (:math:`p_{ind}`) to all the DOF (:math:`p_{all}`): .. math:: p_{all} = \mathbf{S} \cdot p_{ind} .. math:: \left \{ \begin{array}{c} in_b \\ in_q \\ out_b \\ out_q \end{array} \right \} = \mathbf{S} \left \{ \begin{array}{c} in_b \\ in_q \\ out_q \end{array} \right \} Since :math:`out_b = in_b`: .. math:: \left \{ \begin{array}{c} in_b \\ in_q \\ out_b \\ out_q \end{array} \right \} = \left [ \begin{array}{ccc} I & 0 & 0 \\ 0 & I & 0 \\ I & 0 & 0 \\ 0 & 0 & I \\ \end{array} \right ] \left \{ \begin{array}{c} in_b \\ in_q \\ out_q \end{array} \right \} The following code takes advantage of the fact that the b-set is ordered first for both models. .. code:: ipython3 m = maa["in"].shape[0] n = maa["out"].shape[0] nb = np.count_nonzero(b["in"]) S = {} S["in"] = np.block( [ [np.eye(nb), np.zeros((nb, m + n - 2 * nb))], [np.zeros((m - nb, nb)), np.eye(m - nb), np.zeros((m - nb, n - nb))], ] ) S["out"] = np.block( [ [np.eye(nb), np.zeros((nb, m + n - 2 * nb))], [np.zeros((n - nb, m)), np.eye(n - nb)], ] ) S["tot"] = np.vstack((S["in"], S["out"])) Apply the S matrix to couple the models together: .. code:: ipython3 mc = S["tot"].T @ la.block_diag(maa["in"], maa["out"]) @ S["tot"] kc = S["tot"].T @ la.block_diag(kaa["in"], kaa["out"]) @ S["tot"] bc = S["tot"].T @ la.block_diag(baa["in"], baa["out"]) @ S["tot"] As a check, compare coupled model frequencies to Nastran: .. code:: ipython3 lam, phi = la.eigh(kc, mc) freqsys = np.sqrt(abs(lam)) / 2 / np.pi eigen = bulk.rdeigen(srcdir / "assemble.out") np.allclose(freqsys[6:], eigen[0]["cycles"][6:]) .. code-block:: none True The system frequencies up to 110 Hz: .. code:: ipython3 pv = eigen[0]["cycles"] <= 110.0 eigen[0][["cycles"]][pv] .. raw:: html
cycles
1 0.00
2 0.00
3 0.00
4 0.00
5 0.00
6 0.00
7 1.70
8 1.77
9 1.86
10 3.42
11 7.02
12 7.03
13 10.72
14 10.98
15 13.87
16 14.39
17 14.65
18 15.19
19 25.20
20 25.31
21 29.13
22 42.44
23 43.09
24 46.89
25 47.81
26 69.45
27 87.46
28 97.91
29 100.85
Solve the equations of motion the old fashioned way ''''''''''''''''''''''''''''''''''''''''''''''''''' We’ll solve the equations of motion using `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ and recover responses for both components. Note that upstream data recovery could be done easily with the “mug1”, “mef1” and “mes1” matrices loaded from the op4 files above. .. code:: ipython3 fs = ode.SolveUnc(mc, bc, kc, pre_eig=True) force = S["in"].T @ force_in sol = fs.fsolve(force, freq, incrb="avd") # keep all rb responses for demo # recover displacements, velocities, and accelerations for both components: d, v, a = {}, {}, {} ifforce = {} for id in ids: d[id] = S[id] @ sol.d v[id] = S[id] @ sol.v a[id] = S[id] @ sol.a ifforce[id] = (maa[id] @ a[id] + baa[id] @ v[id] + kaa[id] @ d[id])[:nb] Do some sanity checks on the results. The boundary accelerations, velocities, and displacements should match, and the boundary forces from one component on the other should be equal and opposite. .. code:: ipython3 print(f'{np.allclose(a["in"][:nb], a["out"][:nb]) = }') print(f'{np.allclose(v["in"][:nb], v["out"][:nb]) = }') print(f'{np.allclose(d["in"][:nb], d["out"][:nb]) = }') print() print(f'{np.allclose(-ifforce["out"], ifforce["in"] - force[:nb]) = }') print(f'{abs(ifforce["out"] + (ifforce["in"] - force[:nb])).max() = }') pv = freq >= 1.0 print() print(f'For freq >= 1.0:\n\t{np.allclose(-ifforce["out"][:, pv], ifforce["in"][:, pv] - force[:nb, pv]) = }') .. code-block:: none np.allclose(a["in"][:nb], a["out"][:nb]) = True np.allclose(v["in"][:nb], v["out"][:nb]) = True np.allclose(d["in"][:nb], d["out"][:nb]) = True np.allclose(-ifforce["out"], ifforce["in"] - force[:nb]) = False abs(ifforce["out"] + (ifforce["in"] - force[:nb])).max() = 0.0029273294523193272 For freq >= 1.0: np.allclose(-ifforce["out"][:, pv], ifforce["in"][:, pv] - force[:nb, pv]) = True *Note:* In the lowest frequencies, we lose some numerical precision on the boundary forces if we set ``incrb = "avd"`` in the call to `fsolve <../modules/generated/pyyeti.ode.SolveUnc.fsolve.html#pyyeti.ode.SolveUnc.fsolve>`__ above (the default). That’s because the rigid-body displacements will grow to infinity as frequency approaches zero. This is shown in the next plot; both linear and log scales are shown for the x-axis scales because it is interesting and informative. We’ll explore a similar phenomenon a little more later in the NT context. For now, we’ll note that the forces are still very close and that above 1.0 Hz, even with ``incrb = "avd"``, the ``np.allclose`` check passes. .. code:: ipython3 fig = plt.figure("RB Modal Displacements", clear=True, layout="constrained", figsize=(10, 3)) axes = fig.subplots(1, 2, sharey=True) axes[0].semilogy(freq, abs(d['out'][:6].T)) axes[1].loglog(freq, abs(d['out'][:6].T)) axes[0].set_xlabel("Frequency (Hz)") axes[1].set_xlabel("Frequency (Hz)") axes[0].set_ylabel("RB Modal Displacements (distance)"); .. image:: temp_ntfl_files/temp_ntfl_41_0.png Repeat the above, but with ``incrb = "av"`` in the call to `fsolve <../modules/generated/pyyeti.ode.SolveUnc.fsolve.html#pyyeti.ode.SolveUnc.fsolve>`__: .. code:: ipython3 fs = ode.SolveUnc(mc, bc, kc, pre_eig=True) force = S["in"].T @ force_in sol = fs.fsolve(force, freq, incrb="av") # zero out the rb displacements # recover displacements, velocities, and accelerations for both components: d, v, a = {}, {}, {} ifforce = {} for id in ids: d[id] = S[id] @ sol.d v[id] = S[id] @ sol.v a[id] = S[id] @ sol.a ifforce[id] = (maa[id] @ a[id] + baa[id] @ v[id] + kaa[id] @ d[id])[:nb] print(f'{np.allclose(-ifforce["out"], ifforce["in"] - force[:nb]) = }') print(f'{abs(ifforce["out"] + (ifforce["in"] - force[:nb])).max() = }') fig = plt.figure("RB Modal Displacements", clear=True, layout="constrained", figsize=(10, 3)) axes = fig.subplots(1, 2, sharey=True) axes[0].semilogy(freq, abs(d['out'][:6].T)) axes[1].loglog(freq, abs(d['out'][:6].T)) axes[0].set_xlabel("Frequency (Hz)") axes[1].set_xlabel("Frequency (Hz)") axes[0].set_ylabel("RB Modal Displacements (distance)"); .. code-block:: none np.allclose(-ifforce["out"], ifforce["in"] - force[:nb]) = True abs(ifforce["out"] + (ifforce["in"] - force[:nb])).max() = 1.7281315377588588e-07 .. image:: temp_ntfl_files/temp_ntfl_43_1.png Use the NT approach to couple models together ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ First, we’ll compute the “free acceleration” :math:`A_s`: .. code:: ipython3 fs_in = ode.FreqDirect(maa["in"], baa["in"], kaa["in"]) As = fs_in.fsolve(force_in, freq).a[:nb] Next, we’ll compute the apparent masses and call `pyyeti.frclim.ntfl <../modules/generated/pyyeti.frclim.ntfl.html#pyyeti.frclim.ntfl>`__ to compute coupled system responses. For illustration and as a check on the routine, we’ll use both methods within `pyyeti.frclim.calcAM <../modules/generated/pyyeti.frclim.calcAM.html#pyyeti.frclim.calcAM>`__ to compute the apparent masses. Those two methods are labeled as ``cbtf`` and ``fsolve`` here; for more information, see `pyyeti.frclim.calcAM <../modules/generated/pyyeti.frclim.calcAM.html#pyyeti.frclim.calcAM>`__. Note that we’ll get a warning from `pyyeti.frclim.calcAM <../modules/generated/pyyeti.frclim.calcAM.html#pyyeti.frclim.calcAM>`__ about switching to the `pyyeti.ode.FreqDirect <../modules/generated/pyyeti.ode.FreqDirect.html#pyyeti.ode.FreqDirect>`__ solver. See above for more discussion on this. .. code:: ipython3 AM = {} NT = {} for method in ("cbtf", "fsolve"): AM[method] = {} for id in ids: if method == "cbtf": drm = np.arange(nb) else: drm = np.zeros((nb, maa[id].shape[0])) drm[:, b[id]] = np.eye(nb) AM[method][id] = frclim.calcAM([maa[id], baa[id], kaa[id], drm], freq) NT[method] = frclim.ntfl(AM[method]["in"], AM[method]["out"], As, freq) .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/pyyeti/envs/latest/lib/python3.12/site-packages/pyyeti/frclim.py:231: RuntimeWarning: Switching from `SolveUnc` to `FreqDirect` because complex eigensolver failed; see messages above. Solution may be slow. warnings.warn( Perform some sanity checks and comparisons against the “CLA” (old fashioned) approach. Again, we see some loss of numerical precision at the lowest frequencies as we noted above for the CLA approach when we used ``incrb="dav"``. We’ll discuss this below. .. code:: ipython3 for method in NT: print() print(method) print(f'\t{np.allclose(ifforce["out"][:, pv], NT[method].F[:, pv]) = }') print(f'\t{np.allclose(a["out"][:nb, pv], NT[method].A[:, pv]) = }') print(f'\t{abs(ifforce["out"] - NT[method].F).max() = }') print(f'\t{abs(a["out"][:nb] - NT[method].A).max() = }') .. code-block:: none cbtf np.allclose(ifforce["out"][:, pv], NT[method].F[:, pv]) = True np.allclose(a["out"][:nb, pv], NT[method].A[:, pv]) = True abs(ifforce["out"] - NT[method].F).max() = 4.8351414481562216e-05 abs(a["out"][:nb] - NT[method].A).max() = 2.4382675478551441e-05 fsolve np.allclose(ifforce["out"][:, pv], NT[method].F[:, pv]) = True np.allclose(a["out"][:nb, pv], NT[method].A[:, pv]) = True abs(ifforce["out"] - NT[method].F).max() = 0.00013002828172092862 abs(a["out"][:nb] - NT[method].A).max() = 2.1016948736907182e-05 As noted above, we lose some numerical precision in the lowest frequencies. In the NT approach, this appears to be due to very large apparent mass values. To illustrate, here are the apparent mass curves for both approaches for the ``x`` DOF of the first boundary node for the Source model. It’s plotted on two different scales because it’s interesting and informative. .. code:: ipython3 fig = plt.figure("AM", clear=True, layout="constrained", figsize=(10, 3)) axes = fig.subplots(1, 2, sharey=True) for method in AM: axes[0].semilogy(freq, abs(AM[method]["in"][0, :, 0]), label=method) axes[1].loglog(freq, abs(AM[method]["in"][0, :, 0]), label=method) axes[0].legend() axes[1].legend() axes[0].set_xlabel("Frequency (Hz)") axes[1].set_xlabel("Frequency (Hz)") axes[0].set_ylabel("X Apparent Mass (mass)"); .. image:: temp_ntfl_files/temp_ntfl_51_0.png Similarly, for the Load model: .. code:: ipython3 fig = plt.figure("AM", clear=True, layout="constrained", figsize=(10, 3)) axes = fig.subplots(1, 2, sharey=True) for method in AM: axes[0].semilogy(freq, abs(AM[method]["out"][0, :, 0]), label=method) axes[1].loglog(freq, abs(AM[method]["out"][0, :, 0]), label=method) axes[0].legend() axes[1].legend() axes[0].set_xlabel("Frequency (Hz)") axes[1].set_xlabel("Frequency (Hz)") axes[0].set_ylabel("X Apparent Mass (mass)"); .. image:: temp_ntfl_files/temp_ntfl_53_0.png One way to think about the apparent mass is that it is the set of forces required to enforce a unit acceleration on one boundary DOF while restraining all other boundary DOF. (That would be one column of the apparent mass matrix at a particular frequency.) In other words, in :math:`f = m a`, if :math:`a = \mathbf{I}`, then the apparent mass :math:`m` is equal to the force :math:`f`. As frequency approaches zero, the structure will respond statically to the slowly oscillating sinusoidal forces and the relation between the displacement and the force is: .. math:: F \approx K x Since :math:`a = -\omega^2 x`: .. math:: F \approx \frac{1}{-\omega^2} K a So, as frequency approaches zero, the apparent mass will increase proportionally to :math:`1/\omega^2`, which explains the straight line slope on the log-log plot (see also the next plot). Note that for statically-determinate interfaces, the apparent mass simply approaches the rigid-body mass as frequency approaches 0.0 Hz since the system will move as a rigid-body. However, even down to 0.01 Hz, the results for the NT approaches are still quite close, and above 1.0 Hz, all the ``np.allclose`` checks return True. For illustration, here is a plot showing the slope of the apparent mass is indeed proportional to :math:`1/\omega^2`: .. code:: ipython3 fig = plt.figure(figsize=(5, 3)) axes = fig.subplots(1, 1) axes.semilogy(freq, abs(AM["cbtf"]["out"][0, :, 0]), label="Apparent Mass") axes.loglog(freq, 1/freq**2, label="1/freq**2") axes.legend() axes.set_xlabel("Frequency (Hz)") axes.set_ylabel("X Apparent Mass (mass)"); .. image:: temp_ntfl_files/temp_ntfl_55_0.png Plot interface forces and accelerations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For illustration, we’ll plot some of the interface forces and boundary accelerations. We can see that NT responses are the same as the CLA responses. First, compare the interface forces for the 4th boundary node (as an example): .. code:: ipython3 fig = plt.figure("ifforce", clear=True, layout="constrained", figsize=(8, 11)) frcaxes = fig.subplots(6, 1, sharex=True) ylabels = ["$F_X$", "$F_Y$", "$F_Z$", "$M_X$", "$M_Y$", "$M_Z$"] node = 4 for j, ax in enumerate(frcaxes): row = j + (node - 1) * 6 which = "Force" if ylabels[j].startswith("$F") else "Moment" for method in NT: ax.semilogy( freq, abs(NT[method].F[row]), label=f"NT[{method!r}] I/F {which}" ) ax.semilogy(freq, abs(ifforce["out"][row]), "--", label=f"CLA I/F {which}") ax.legend() ax.set_ylabel(ylabels[j]) frcaxes[0].set_title(f"Magnitude of I/F Force on 'Outboard', Boundary Node #{node}") frcaxes[-1].set_xlabel("Frequency (Hz)"); .. image:: temp_ntfl_files/temp_ntfl_58_0.png Next, compare the interface accelerations for the 4th boundary node. For reference, the free acceleration for this node is also plotted: .. code:: ipython3 fig = plt.figure("ifacce", clear=True, layout="constrained", figsize=(8, 11)) acceaxes = fig.subplots(6, 1, sharex=True) ylabels = ["$A_{X}$", "$A_{Y}$", "$A_{Z}$", "$A_{RX}$", "$A_{RY}$", "$A_{RZ}$"] for j, ax in enumerate(acceaxes): row = j + (node - 1) * 6 for method in NT: ax.semilogy( freq, abs(NT[method].A[row]), label=f"NT[{method!r}] Coupled Acce" ) ax.semilogy(freq, abs(a["out"][row]), "--", label="CLA Coupled Acce") ax.semilogy(freq, abs(As[row]), label="Free Acce") ax.legend() ax.set_ylabel(ylabels[j]) acceaxes[0].set_title(f"Magnitude of I/F Acceleration, Boundary Node #{node}") acceaxes[-1].set_xlabel("Frequency (Hz)"); .. image:: temp_ntfl_files/temp_ntfl_60_0.png