From 77f0adf42de492c6d63cab3ac6aa246ab4e8e4d4 Mon Sep 17 00:00:00 2001 From: Colin Caprani <31228546+ccaprani@users.noreply.github.com> Date: Sun, 31 Mar 2024 16:53:52 +1100 Subject: [PATCH] Deflection shape fixes for all element types (#77) --- src/pycba/analysis.py | 2 +- src/pycba/beam.py | 13 ++- src/pycba/inf_lines.py | 78 ++++++++----- src/pycba/load.py | 46 ++++---- src/pycba/results.py | 28 ++++- tests/test_basic.py | 239 +++++++++++++++++++++++++++++++++++++++- tests/test_inf_lines.py | 2 +- 7 files changed, 340 insertions(+), 68 deletions(-) diff --git a/src/pycba/analysis.py b/src/pycba/analysis.py index 718d76c..077545a 100644 --- a/src/pycba/analysis.py +++ b/src/pycba/analysis.py @@ -223,7 +223,7 @@ def _forces(self) -> np.ndarray: for i in range(self._n): dof_i = 2 * i - fmbr = self._beam.get_cnl(i) + fmbr = self._beam.get_ref(i) # Cumulatively apply forces in opposite direction f[dof_i : dof_i + 4] -= fmbr return f diff --git a/src/pycba/beam.py b/src/pycba/beam.py index d9f0670..c4959c5 100644 --- a/src/pycba/beam.py +++ b/src/pycba/beam.py @@ -280,9 +280,10 @@ def get_local_span_coords(self, pos: float) -> (int, float): return ispan, pos_in_span - def get_cnl(self, i_span: int) -> LoadCNL: + def get_ref(self, i_span: int) -> LoadCNL: """ - Returns Consistent Nodal Loads for the member + Returns Released End Forces for the member; that is, the Consistent Nodal Loads + modified for the element type (i.e. releases) Parameters ---------- @@ -291,18 +292,18 @@ def get_cnl(self, i_span: int) -> LoadCNL: Returns ------- - cnl : LoadCNL + ref : LoadCNL The totalled CNL object for the member, considering all loads. """ - cnl = np.zeros(4) + ref = np.zeros(4) L = self.mbr_lengths[i_span] eType = self.mbr_eletype[i_span] for load in self._loads: if load.i_span == i_span: - cnl += load.get_cnl(L, eType) - return cnl + ref += load.get_ref(L, eType) + return ref def get_span_k(self, i_span: int) -> np.ndarray: """ diff --git a/src/pycba/inf_lines.py b/src/pycba/inf_lines.py index 50c4740..fb91738 100644 --- a/src/pycba/inf_lines.py +++ b/src/pycba/inf_lines.py @@ -118,45 +118,67 @@ def get_il(self, poi: float, load_effect: str) -> (np.ndarray, np.ndarray): if not self.vResults: self.create_ils() - x = self.vResults[0].results.x - dx = x[2] - x[1] - idx = np.where(np.abs(x - poi) <= dx * 1e-6)[0][0] - #idx = np.abs(x - poi).argmin() - + x = self.vResults[0].results.x npts = len(self.vResults) eta = np.zeros(npts) - # - # Getting the correct vertical reaction is tricky - # Should be relatively easy to extend to get moment reactions too. + # Preparations for reaction ILs # # Get vector of the node locations node_locations = np.cumsum(np.insert(self.ba.beam.mbr_lengths, 0, 0)) - # The indices of the supported vertical DOFs wrt the node locations vector - vert_sup_dof_idx = np.where(np.array(self.ba._beam.restraints)[::2]==-1)[0] - # The locations then of these vertical supports - vert_sup_locs = node_locations[vert_sup_dof_idx] - # The index of the closest vertical support - closest_vert_sup_idx = np.abs(vert_sup_locs-poi).argmin() - # And its value - closest_vert_sup = vert_sup_locs[closest_vert_sup_idx] - # And now the indixe of this support in the node locations vector - node_idx = np.where(node_locations==closest_vert_sup)[0][0] - # And hence its index in the overall DOFs vector - dof_idx = 2*node_idx - - # Now we must link the supported DOF to the index in the BeamAnalysis reactions vector + # Link the supported DOF to the index in the BeamAnalysis reactions vector idx_mask = np.zeros_like(self.ba._beam.restraints) idx_mask[np.where(np.array(self.ba._beam.restraints)==-1)] = np.arange(self.ba.beam.no_fixed_restraints) - # And finally the index of the vertical support nearest the POI in the reactions vector - vert_sup_idx = idx_mask[dof_idx] - for i, res in enumerate(self.vResults): - if load_effect == "V": + #idx = np.abs(x - poi).argmin() + + if load_effect.upper() == "V": + dx = x[2] - x[1] + idx = np.where(np.abs(x - poi) <= dx * 1e-6)[0][0] + for i, res in enumerate(self.vResults): eta[i] = res.results.V[idx] - elif load_effect == "R": + + elif load_effect.upper() == "R": + # + # Getting the correct reaction is tricky + # + # The indices of the supported DOFs wrt the node locations vector + vert_sups_dof_idx = np.where(np.array(self.ba._beam.restraints)[::2]==-1)[0] + # The locations then of these supports + vert_sups_locs = node_locations[vert_sups_dof_idx] + # The index of the closest support + closest_vert_sup_idx = np.abs(vert_sups_locs-poi).argmin() + # And its value + closest_vert_sup = vert_sups_locs[closest_vert_sup_idx] + # And now the index of this support in the node locations vector + vert_sup_node_idx = np.where(node_locations==closest_vert_sup)[0][0] + # And hence its index in the overall DOFs vector + vert_sup_dof_idx = 2*vert_sup_node_idx + # And finally the index of the support nearest the POI in the reactions vector + vert_sup_idx = idx_mask[vert_sup_dof_idx] + + for i, res in enumerate(self.vResults): eta[i] = res.R[vert_sup_idx] - else: + + elif load_effect.upper() == "MR": + # + # Follows the same logic for the vertical reaction + # + mt_sups_dof_idx = np.where(np.array(self.ba._beam.restraints)[1::2]==-1)[0] + mt_sups_locs = node_locations[mt_sups_dof_idx] + closest_mt_sup_idx = np.abs(mt_sups_locs-poi).argmin() + closest_mt_sup = mt_sups_locs[closest_mt_sup_idx] + mt_sup_node_idx = np.where(node_locations==closest_mt_sup)[0][0] + mt_sup_dof_idx = 2*mt_sup_node_idx+1 + mt_sup_idx = idx_mask[mt_sup_dof_idx] + + for i, res in enumerate(self.vResults): + eta[i] = res.R[mt_sup_idx] + + else: + dx = x[2] - x[1] + idx = np.where(np.abs(x - poi) <= dx * 1e-6)[0][0] + for i, res in enumerate(self.vResults): eta[i] = res.results.M[idx] return (np.array(self.pos), eta) diff --git a/src/pycba/load.py b/src/pycba/load.py index 4450605..bd3894c 100644 --- a/src/pycba/load.py +++ b/src/pycba/load.py @@ -182,7 +182,7 @@ def __init__(self, i_span: int): """ self.i_span = i_span - def get_cnl(self, L, eType): + def get_cnl(self, L): # Enforce virtual base class raise NotImplementedError @@ -215,29 +215,27 @@ def H(self, v: np.ndarray, value: float = 0.0) -> np.ndarray: """ return np.heaviside(v, value) - def released_end_forces(self, cnl: LoadCNL, L: float, eType: int) -> LoadCNL: + + def get_ref(self, L: float, eType: int) -> LoadCNL: """ - The released end forces for each element type: converts the Consistent Nodal - Loads of the applied loading to the correct nodal loading depending on the - element type. + Returns the Released End Forces for a span of length L of element eType: + converts the Consistent Nodal Loads of the applied loading to the correct nodal + loading depending on the element type. Parameters ---------- - cnl : LoadCNL - The nodal loading statically consistent with the externally applied loads. L : float - The length of the member. + The length of the member eType : int - The element type. + The member element type Returns ------- LoadCNL - The nodal loads to be applied in the analysis, consistent with the element - type. - + Released End Forces for this load type: the nodal loads to be applied in + the analysis, consistent with the element type. """ - + cnl = self.get_cnl(L, eType) ref = np.zeros(4) fm = 6 / (4 * L) # flexibility coeff for moment @@ -252,10 +250,10 @@ def released_end_forces(self, cnl: LoadCNL, L: float, eType: int) -> LoadCNL: ref[2] = -fm * cnl.Ma ref[3] = 0.5 * cnl.Ma elif eType == 4: # keep only vertical, remove moments - ref[0] = 0 - ref[1] = 1.0 * cnl.Va - ref[2] = 0 - ref[3] = 1.0 * cnl.Va + ref[0] = -(cnl.Ma+cnl.Mb)/L + ref[1] = 1.0 * cnl.Ma + ref[2] = (cnl.Ma+cnl.Mb)/L + ref[3] = 1.0 * cnl.Mb else: # no nothing if it is FF pass @@ -319,8 +317,7 @@ def get_cnl(self, L: float, eType: int) -> LoadCNL: Ma=w * L**2 / 12.0, Mb=-w * L**2 / 12.0, ) - - return self.released_end_forces(cnl, L, eType) + return cnl def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults: """ @@ -416,8 +413,7 @@ def get_cnl(self, L, eType) -> LoadCNL: Ma=P * a * b**2 / L**2, Mb=-P * a**2 * b / L**2, ) - # implicit conversion to tuple in correct order - return self.released_end_forces(cnl, L, eType) + return cnl def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults: """ @@ -512,7 +508,7 @@ def get_cnl(self, L, eType) -> LoadCNL: Mb=-(w * c / L**2) * (t * s**2 + (t - 2 * s) * c**2 / 12), ) # implicit conversion to tuple in correct order - return self.released_end_forces(cnl, L, eType) + return cnl def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults: """ @@ -608,8 +604,7 @@ def get_cnl(self, L, eType) -> LoadCNL: Ma=Ma, Mb=Mb, ) - # implicit conversion to tuple in correct order - return self.released_end_forces(cnl, L, eType) + return cnl def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults: """ @@ -690,8 +685,7 @@ def get_cnl(self, L, eType) -> LoadCNL: Ma=(m * b / L**2) * (2 * a - b), Mb=(m * a / L**2) * (2 * b - a), ) - # implicit conversion to tuple in correct order - return self.released_end_forces(cnl, L, eType) + return cnl def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults: """ diff --git a/src/pycba/results.py b/src/pycba/results.py index ebeeafe..4738d4a 100644 --- a/src/pycba/results.py +++ b/src/pycba/results.py @@ -8,7 +8,7 @@ import matplotlib.pyplot as plt from scipy import integrate from .beam import Beam -from .load import MemberResults, LoadMaMb +from .load import MemberResults, LoadMaMb, LoadCNL from copy import deepcopy @@ -105,7 +105,7 @@ def _member_analysis(self, beam: Beam, d: np.ndarray) -> List[MemberResults]: fmbr = np.zeros(4) for j in range(4): fmbr[j] = np.sum(kb[j][:] * dmbr[:]) - fmbr += beam.get_cnl(i) + fmbr += beam.get_ref(i) res = self._member_values(beam, i, fmbr, dmbr) # Shift x vals by location of mbr starting point res.x += sumL @@ -138,6 +138,9 @@ def _member_values( """ L = beam.mbr_lengths[i_span] + EI = beam.mbr_EIs[i_span] + etype = beam.mbr_eletype[i_span] + dx = L / self.npts x = np.zeros(self.npts + 3) x[1 : self.npts + 2] = dx * np.arange(0, self.npts + 1) @@ -148,16 +151,31 @@ def _member_values( res = MaMb.get_mbr_results(x, L) # Now get the results for all the applied loads on a simple span + Ma=0 + Mb=0 for load in beam._loads: if load.i_span != i_span: continue res += load.get_mbr_results(x, L) + cnl = load.get_cnl(L, etype) + Ma += cnl.Ma + Mb += cnl.Mb + + # Check element type for any released displacements + R0 = d[1] + theta = 0 + phi_i = 0 + # Account for end release if joint not already rotating + if etype != 1 and abs(R0) < 1e-6 : + theta = (d[2]-d[0])/L + phi_i = (L/(3*EI))*(-(f[1]-0.5*f[3])+(Ma-0.5*Mb)) + + R0 -= phi_i - theta # And superimpose end displacements using Moment-Area - h = L / self.npts - EI = beam.mbr_EIs[i_span] + h = L / self.npts - R = integrate.cumtrapz(res.M[1:-1], dx=h, initial=0) / EI + d[1] + R = integrate.cumtrapz(res.M[1:-1], dx=h, initial=0) / EI + R0 D = integrate.cumtrapz(R, dx=h, initial=0) + d[0] res.R[1:-1] = R diff --git a/tests/test_basic.py b/tests/test_basic.py index 45b4ede..eb187ea 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -7,6 +7,243 @@ import pycba as cba +def test_1span_ee(): + """ + Test fixed-fixed beam with point load in the middle + """ + + P = 10 # kN + L = 10 # m + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, -1, -1, -1] + LM = [[1, 2, P, 0.5*L, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM) + out = beam_analysis.analyze() + assert out == 0 + + Ma = beam_analysis.beam_results.results.M[1] + Mb = beam_analysis.beam_results.results.M[-2] + Mc = beam_analysis.beam_results.results.M[51] + + assert Ma == pytest.approx(-P*L/8) + assert Mb == pytest.approx(-P*L/8) + assert Mc == pytest.approx(P*L/8) + + +def test_1span_ep(): + """ + Test fixed-pinned beam with point load in the middle + """ + + P = 10 # kN + L = 10 # m + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, -1, -1, 0] + LM = [[1, 2, P, 0.5*L, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM) + out = beam_analysis.analyze() + assert out == 0 + + Ma = beam_analysis.beam_results.results.M[1] + Mb = beam_analysis.beam_results.results.M[-2] + Mc = beam_analysis.beam_results.results.M[51] + + assert Ma == pytest.approx(-3*P*L/16) + assert Mb == pytest.approx(0) + assert Mc == pytest.approx(5*P*L/32) + + +def test_1span_ep_eletype2(): + """ + Test fixed-pinned beam with point load in the middle, using eleType 2 + """ + + P = 10 # kN + L = 10 # m + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, -1, -1, -1] # Notice, fixed-fixed supports + LM = [[1, 2, P, 0.5*L, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM, eletype=[2]) + out = beam_analysis.analyze() + assert out == 0 + + Ma = beam_analysis.beam_results.results.M[1] + Mb = beam_analysis.beam_results.results.M[-2] + Mc = beam_analysis.beam_results.results.M[51] + + assert Ma == pytest.approx(-3*P*L/16) + assert Mb == pytest.approx(0) + assert Mc == pytest.approx(5*P*L/32) + + +def test_1span_pe(): + """ + Test pinned-fixed beam with point load in the middle + """ + + P = 10 # kN + L = 10 # m + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, 0, -1, -1] + LM = [[1, 2, P, 0.5*L, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM) + out = beam_analysis.analyze() + assert out == 0 + + Ma = beam_analysis.beam_results.results.M[1] + Mb = beam_analysis.beam_results.results.M[-2] + Mc = beam_analysis.beam_results.results.M[51] + + assert Ma == pytest.approx(0) + assert Mb == pytest.approx(-3*P*L/16) + assert Mc == pytest.approx(5*P*L/32) + + +def test_1span_pe_eletype3(): + """ + Test pinned-fixed beam with point load in the middle, using eleType 3 + """ + + P = 10 # kN + L = 10 # m + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, -1, -1, -1] # Notice, fixed-fixed supports + LM = [[1, 2, P, 0.5*L, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM, eletype=[3]) + out = beam_analysis.analyze() + assert out == 0 + + Ma = beam_analysis.beam_results.results.M[1] + Mb = beam_analysis.beam_results.results.M[-2] + Mc = beam_analysis.beam_results.results.M[51] + + assert Ma == pytest.approx(0) + assert Mb == pytest.approx(-3*P*L/16) + assert Mc == pytest.approx(5*P*L/32) + + +def test_1span_pp(): + """ + Test pinned-pinned beam with point load in the middle + """ + + P = 10 # kN + L = 10 # m + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, 0, -1, 0] + LM = [[1, 2, P, 0.5*L, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM) + out = beam_analysis.analyze() + assert out == 0 + + Ma = beam_analysis.beam_results.results.M[1] + Mb = beam_analysis.beam_results.results.M[-2] + Mc = beam_analysis.beam_results.results.M[51] + + assert Ma == pytest.approx(0) + assert Mb == pytest.approx(0) + assert Mc == pytest.approx(P*L/4) + + +def test_1span_pp_eletype4(): + """ + Test pinned-pinned beam with point load in the middle, using eleType 4 + """ + + P = 10 # kN + L = 10 # m + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, -1, -1, -1] # Notice, fixed-fixed supports + LM = [[1, 2, P, 0.5*L, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM, eletype=[4]) + out = beam_analysis.analyze() + assert out == 0 + + Ma = beam_analysis.beam_results.results.M[1] + Mb = beam_analysis.beam_results.results.M[-2] + Mc = beam_analysis.beam_results.results.M[51] + + assert Ma == pytest.approx(0) + assert Mb == pytest.approx(0) + assert Mc == pytest.approx(P*L/4) + + +def get_1span_beam_def(etype): + P = 10 # kN + L = 10 # m + a = 0.25*L + EI = 30 * 600e7 * 1e-6 # kNm2 + R = [-1, -1, -1, -1] # Notice, fixed-fixed supports + LM = [[1, 2, P, a, 0]] + + beam_analysis = cba.BeamAnalysis([L], EI, R, LM, eletype=[etype]) + beam_analysis.analyze() + + d = beam_analysis.beam_results.results.D + dmax = min(d) + + return P, L, EI, a, dmax + + +def test_1span_def_ff(): + """ + Test fixed-fixed beam deflection for off-centre point load + """ + + P, L, EI, aa, dmax = get_1span_beam_def(etype=1) + + # a>b + b = aa + a = L-b + ymax = -(2*P*a**3*b**2)/(3*(3*a+b)**2*EI) + + assert dmax == pytest.approx(ymax, abs=1e-6) + + +def test_1span_def_fp(): + """ + Test fixed-pinned beam deflection for off-centre point load + """ + P, L, EI, aa, dmax = get_1span_beam_def(etype=2) + + a = aa + b = L-a + ymax = -(P*a**2*b)/(6*EI)*(b/(3*L-a))**0.5 + + assert dmax == pytest.approx(ymax, abs=1e-6) + + +def test_1span_def_pf(): + """ + Test pinned-fixed beam deflection for off-centre point load + """ + P, L, EI, aa, dmax = get_1span_beam_def(etype=3) + + b = aa + ymax = -(P*b)/(3*EI)*(L**2-b**2)**3/(3*L**2-b**2)**2 + + assert dmax == pytest.approx(ymax, abs=1e-6) + + +def test_1span_def_pp(): + """ + Test pinned-pinned beam deflection for off-centre point load + """ + P, L, EI, aa, dmax = get_1span_beam_def(etype=4) + + a = aa + ymax = -3**0.5*P*a*(L**2-a**2)**1.5/(27*EI*L) + + assert dmax == pytest.approx(ymax, abs=1e-6) + + def test_2span_udl(): """ Execute a two-span beam analysis and check the reaction results. @@ -101,7 +338,7 @@ def test_3span_diff_settlement(): dmax = max(beam_analysis.beam_results.results.D) dmin = min(beam_analysis.beam_results.results.D) assert [dmax, dmin] == pytest.approx( - [0.0002778734578837245, -0.004648274177216889], abs=1e-6 + [0.0002778734578837245, -0.004648274177216889], abs=1e-5 ) diff --git a/tests/test_inf_lines.py b/tests/test_inf_lines.py index 64ef774..300a9eb 100644 --- a/tests/test_inf_lines.py +++ b/tests/test_inf_lines.py @@ -177,7 +177,7 @@ def test_parse_beam_notation(): assert eType == [1, 2, 1, 1] -def test_distcretization(): +def test_discretization(): """ Confirm that poi rounding to find the closest idx for ILs works """