From 54ac71a40414bf4d7cbad99e8fb6a1bfe3fdb962 Mon Sep 17 00:00:00 2001 From: Petr Krysl Date: Sat, 23 Dec 2023 11:22:33 -0800 Subject: [PATCH] update twisted_beam_tut --- tutorials/index.md | 12 +++++----- tutorials/twisted_beam_tut.jl | 41 ++++++++++++++++++++--------------- tutorials/twisted_beam_tut.md | 7 +++--- 3 files changed, 33 insertions(+), 27 deletions(-) diff --git a/tutorials/index.md b/tutorials/index.md index d202378..d369033 100644 --- a/tutorials/index.md +++ b/tutorials/index.md @@ -6,18 +6,18 @@ - [Beam vibration](beam_modal_tut.md) A good place to start, it goes into the basics. - [Aircraft frame test rig: Geometry](garteur_geometry_tut.md) The next four tutorials develop a dynamic investigation of a test rig. -- [Aircraft frame test rig: Visualization](garteur_geometry_vis_tut.md) +- [Aircraft frame test rig: Visualization](garteur_geometry_vis_tut.md) - [Aircraft frame test rig: Modal analysis](garteur_modal_tut.md) - [Aircraft frame test rig: Harmonic Vibration Analysis](garteur_hva_tut.md) -- [Prestressed column: Fundamental natural frequency](prestressed_column_modal_tut.md) -- [L-shaped aluminum frame: Vibration under applied loading](argyris_frame_modal_tut.md) -- [Free vibration of steel ring: Introduction](circle_modal_tut.md) +- [Prestressed column: Fundamental natural frequency](prestressed_column_modal_tut.md) Illustration of the effect of pre-stress on the natural frequencies. +- [L-shaped aluminum frame: Vibration under applied loading](argyris_frame_modal_tut.md) Interaction of two modes of buckling. +- [Free vibration of steel ring: Introduction](circle_modal_tut.md) The next three tutorials discuss a well known benchmark. - [Free vibration of steel ring: Beam model convergence](circle_modal_conv_tut.md) - [Free vibration of steel ring: Solid model convergence](circle_3d_modal_conv_tut.md) -- [Fast Lagrange top](fast_top_tut.md) +- [Fast Lagrange top](fast_top_tut.md) Illustration of nonlinear rotation dynamics. - Beam with time-dependent axial force (Mathieu problem) ## Shell Statics -- [Twisted beam deflection](twisted_beam_tut.md) +- [Twisted beam deflection](twisted_beam_tut.md) A well-known benchmark problem solved with triangular shell elements. diff --git a/tutorials/twisted_beam_tut.jl b/tutorials/twisted_beam_tut.jl index d90f9a9..7d85945 100644 --- a/tutorials/twisted_beam_tut.jl +++ b/tutorials/twisted_beam_tut.jl @@ -2,6 +2,8 @@ # Source code: [`twisted_beam_tut.jl`](twisted_beam_tut.jl) +# Last updated: 12/23/23 + # ## Description # The initially twisted cantilever beam is one of the standard test @@ -25,7 +27,10 @@ # Test Finite Element Accuracy,” Finite Elements in Analysis Design, vol. 11, pp. # 3–20, 1985. -# [2] Simo, J. C., D. D. Fox, and M. S. Rifai, “On a Stress Resultant Geometrically Exact Shell Model. Part II: The Linear Theory; Computational Aspects,” Computational Methods in Applied Mechanical Engineering, vol. 73, pp. 53–92, 1989. +# [2] Simo, J. C., D. D. Fox, and M. S. Rifai, “On a Stress Resultant +# Geometrically Exact Shell Model. Part II: The Linear Theory; Computational +# Aspects,” Computational Methods in Applied Mechanical Engineering, vol. 73, +# pp. 53–92, 1989. # [3] Zupan D, Saje M (2004) On "A proposed standard set of problems to test # finite element accuracy": the twisted beam. Finite Elements in Analysis @@ -87,20 +92,24 @@ params = params_thicker_dir_3 # The mesh is initially generated for a rectangular 2d domain. The element size # can be controlled with these two variables: -nL = 48; -nW = 8; -nL, nW = 8, 4 +nL, nW = 8, 2 +# This represents a fairly coarse mesh. We may try something more refined: +# nL, nW = 4*8, 4*2 fens, fes = T3block(L, W, nL, nW, :a); # The 2d mesh is expanded into a three dimensional domain, and the -# locations of the nodes are tweaked to produce the pre twisted shape. +# locations of the nodes are modified to produce the pre-youtwisted shape. fens.xyz = xyz3(fens) -for i in 1:count(fens) - a = fens.xyz[i, 1] / L * (pi / 2) - y = fens.xyz[i, 2] - (W / 2) - z = fens.xyz[i, 3] - fens.xyz[i, :] = [fens.xyz[i, 1], y * cos(a) - z * sin(a), y * sin(a) + z * cos(a)] +fens = let + for i in 1:count(fens) + x = fens.xyz[i, 1] + a = x / L * (pi / 2) + y = fens.xyz[i, 2] - (W / 2) + z = fens.xyz[i, 3] + fens.xyz[i, :] = [x, y * cos(a) - z * sin(a), y * sin(a) + z * cos(a)] + end + fens end @@ -112,8 +121,7 @@ t3ffm = FEMMShellT3FFModule # The shell elements have a type `FESetShellT3`. The type of the mesh is the # plain isoparametric triangle, which is instructed to delegate to the shell # element. -sfes = FESetShellT3() -accepttodelegate(fes, sfes) +accepttodelegate(fes, FESetShellT3()) # Use the convenience function to make an instance of the finite element model # machine for the T3FF shell. @@ -137,7 +145,6 @@ l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance) for i in 1:6 setebc!(dchi, l1, true, i) end -applyebc!(dchi) numberdofs!(dchi); # Associate the finite element model machine with geometry. The shell @@ -149,7 +156,7 @@ t3ffm.associategeometry!(femm, geom0) # Assemble the system stiffness matrix. K = t3ffm.stiffness(femm, geom0, u0, Rfield0, dchi); -# Though load is a concentrated force applied at the center of the beam. First +# The load is a concentrated force applied at the center of the beam. First # we select the node. nl = selectnode(fens; box = Float64[L L 0 0 0 0], tolerance = tolerance) # Next we create a mesh of the loaded boundary (i.e. the selected node), so that @@ -165,8 +172,8 @@ fi = ForceIntensity(v); # Finally we computed the load vector corresponding to the force intensity. F = distribloads(lfemm, geom0, dchi, fi, 3); -# Extract the free-free block of the matrix. And the free block for the right -# and side vector. +# Extract the free-free block of the matrix, and the free block for the +# right-hand side vector. K_ff = matrix_blocked(K, nfreedofs(dchi))[:ff] F_f = vector_blocked(F, nfreedofs(dchi))[:f] @@ -177,7 +184,7 @@ scattersysvec!(dchi, U[:]) # The deflection in the correct direction at the loaded node is now extracted. tipdefl = dchi.values[nl, params.dir][1] -@show tipdefl / params.uex * 100 +@info "Normalized deflection: $(tipdefl / params.uex * 100)%" # Generate a graphical display of resultants. The resultants only make sense in # a coordinate system aligned with the shell surface. Here we construct such a diff --git a/tutorials/twisted_beam_tut.md b/tutorials/twisted_beam_tut.md index df92356..a86c406 100644 --- a/tutorials/twisted_beam_tut.md +++ b/tutorials/twisted_beam_tut.md @@ -176,7 +176,6 @@ l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance) for i in 1:6 setebc!(dchi, l1, true, i) end -applyebc!(dchi) numberdofs!(dchi); ```` @@ -194,7 +193,7 @@ Assemble the system stiffness matrix. K = t3ffm.stiffness(femm, geom0, u0, Rfield0, dchi); ```` -Though load is a concentrated force applied at the center of the beam. First +The load is a concentrated force applied at the center of the beam. First we select the node. ````julia @@ -229,8 +228,8 @@ Finally we computed the load vector corresponding to the force intensity. F = distribloads(lfemm, geom0, dchi, fi, 3); ```` -Extract the free-free block of the matrix. And the free block for the right -and side vector. +Extract the free-free block of the matrix, and the free block for the +right-hand side vector. ````julia K_ff = matrix_blocked(K, nfreedofs(dchi))[:ff]