C_R

3537 Reputation

21 Badges

6 years, 60 days

MaplePrimes Activity


These are Posts that have been published by C_R

This post discusses a solution for modeling a traveling load on Maplesim's Flexible Beam component and provides an example of a bouncing load.

The idea for the above example came from an attempt to reproduce a model of a mass sliding on a beam from MapleSim's model gallery. However, reproducing it using contact components in combination with the Flexible Beam component turned out to be not straightforward, and this will be discussed in the following.

To simulate a traveling load on the Flexible Beam component, one could apply forces at discrete locations for a certain duration. However, the fidelity of this approach is limited by the number of discrete locations, which must be defined using the Flexible Beam Frame component, as well as the way in which the forces are activated.

One potential solution to address the issue of temporal activation of forces is to attach contact elements (such as Rectangle components) at distinct locations along the beam, which are defined by Flexible Beam Frame components, and make contact using a spherical or toroidal contact element. However, this approach also introduces two new problems:

  • An additional bending moment is generated when the load is not applied at the center of the contact element's attachment point, the Flexible Beam Frame component. Depending on the length of the contact element, deformations caused by this moment can be greater than the deformation caused by the force itself when the force is applied at the ends of the contact element. Overall, this unwanted moment makes the simulation unrealistic and must be avoided.
  • When the beam bends, a gap (see below) or an overlap is created between adjacent Rectangle components. If there is a gap, the object exerting a force on the beam can fall through it. Overlaps can create differences in dynamic behavior when the radius of curvature of the beam is on the opposite side of the point of contact.

To avoid these problems, the solution presented here uses an intermediate kinematic chain (encircled in yellow below) that redistributes the contact force on the Rectangle component on two support points (ports to attach Flexible Beam Frame components) in a linear fashion.

 

 

To address gaps, the contact element (Rectangle) attached to the kinematic chain has the same width as the chain and connects to the adjacent contact elements via multibody frames. In the image below, 10 contact elements are laid on top of a single Flexible Beam component, like a belt made out of tiles. The belt has to be pinned to the flexible beam at one location (highlighted in yellow). The location of this fixed point determines how the flexible beam is loaded by tangential contact forces (friction forces) and should be selected carefully.

 

 

Some observations on the attached model:

  • Low damping and high initial potential energy of the ball can result in a failed simulation (due to constraint projection failure). Increasing the number of elastic coordinates has a similar effect. Constraint projection can be turned off in the simulation settings to continue simulation.
  • The bouncing ball excites several eigenmodes at once, causing the beam to wiggle chaotically in combination with the varying bouncing frequency of the ball. A similar looking effect can also be achieved with special initial conditions, as demonstrated with Maple in this excellent post on Euler beams and partial differential equations.
  • Repeated simulations with low damping lead to different results (an indication of chaotic behavior; see three successive simulations below (gold) compared to a saved solution(red)). The moment in the animation when the ball travels backward represents a metastable equilibrium point of the simulation. This makes predictions beyond this point difficult, as the behavior of the system is highly dependent on the model parameters. Whether the reversal is a simulation artifact or can happen in reality remains to be seen. Overall, this example could evolve into a nice experimental fun project for students.

  • Setting the gravitational constant for Mars, everything is different. I could not reproduce the fun factor on Earth. A reason more to stay ,-)

Ball_bouncing_on_a_flexible_beam.msim

 

 

Transfer functions are normally not used with units. Involving units when deriving transfer functions can help identify unit inconsistencies and reduce the likelihood of unit conversion errors.

Maple is already a great help in not having to do this manually. However, the final step of simplification still requires manual intervention, as shown in this example.

Given transfer function

H(s) = 60.*Unit('m'*'kg'/('s'^2*'A'))/(.70805*s^2*Unit('kg'^2*'m'^2/('s'^3*'A'^2))+144.*s*Unit('kg'^2*'m'^2/('s'^4*'A'^2))+0.3675e-4*s^3*Unit('kg'^2*'m'^2/('s'^2*'A'^2)))

H(s) = 60.*Units:-Unit(m*kg/(s^2*A))/(.70805*s^2*Units:-Unit(kg^2*m^2/(s^3*A^2))+144.*s*Units:-Unit(kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units:-Unit(kg^2*m^2/(s^2*A^2)))

(1)

Desired output (derived by hand) where the transfer function is separated in a dimensionless expression and a gain that can be attributed to units with a physical meaning in the context of an application (here displacement per voltage).

H(s) = 60.*Unit('m'/'V')/(.70805*s^2*Unit('s'^2)+144.*s*Unit('s')+0.3675e-4*s^3*Unit('s'^3))

H(s) = 60.*Units:-Unit(m/V)/(.70805*s^2*Units:-Unit(s^2)+144.*s*Units:-Unit(s)+0.3675e-4*s^3*Units:-Unit(s^3))

(2)

is(simplify((H(s) = 60.*Units[Unit](m*kg/(s^2*A))/(.70805*s^2*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*s*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units[Unit](kg^2*m^2/(s^2*A^2))))-(H(s) = 60.*Units[Unit](m/V)/(.70805*s^2*Units[Unit](s^2)+144.*s*Units[Unit](s)+0.3675e-4*s^3*Units[Unit](s^3)))))

true

(3)

Units to factor out in the denominator are Unit('kg'^2*'m'^2/('s'^5*'A'^2)). Quick check:

Unit('m'*'kg'/('s'^2*'A'))/Unit('kg'^2*'m'^2/('s'^5*'A'^2)) = Unit('m'/'V')

Units:-Unit(m*kg/(s^2*A))/Units:-Unit(kg^2*m^2/(s^5*A^2)) = Units:-Unit(m/V)

(4)

simplify(Units[Unit](m*kg/(s^2*A))/Units[Unit](kg^2*m^2/(s^5*A^2)) = Units[Unit](m/V))

Units:-Unit(s^3*A/(m*kg)) = Units:-Unit(s^3*A/(m*kg))

(5)

"Simplification" attempts with the denominator

denom(rhs(H(s) = 60.*Units[Unit](m*kg/(s^2*A))/(.70805*s^2*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*s*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units[Unit](kg^2*m^2/(s^2*A^2)))))

s*(.70805*s*Units:-Unit(kg^2*m^2/(s^3*A^2))+144.*Units:-Unit(kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units:-Unit(kg^2*m^2/(s^2*A^2)))

(6)

collect(s*(.70805*s*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units[Unit](kg^2*m^2/(s^2*A^2))), Unit('kg'^2*'m'^2/('s'^5*'A'^2)))

s*(.70805*s*Units:-Unit(kg^2*m^2/(s^3*A^2))+144.*Units:-Unit(kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units:-Unit(kg^2*m^2/(s^2*A^2)))

(7)

is not effective because all units are wrapped in Unit commands. Example:

Unit('kg'^2*'m'^2/('s'^2*'A'^2))

Units:-Unit(kg^2*m^2/(s^2*A^2))

(8)

Expand does not expand the argument of Unit commands.

expand(Units[Unit](kg^2*m^2/(s^2*A^2))); lprint(%)

Units:-Unit(kg^2*m^2/(s^2*A^2))

 

Units:-Unit(kg^2*m^2/s^2/A^2)

 

NULL

C1: Expanding Unit command

An expand facility could be a solution that expands a Unit command with combined units to a product of separate Unit commands.

When all units are expanded in a separate Unit command, collect or factor can be used to collect units:

.70805*s*Unit('kg')^2*Unit('m')^2/(Unit('A')^2*Unit('s')^3)+144.*Unit('kg')^2*Unit('m')^2/(Unit('A')^2*Unit('s')^4)+0.3675e-4*s^2*Unit('kg')^2*Unit('m')^2/(Unit('A')^2*Unit('s')^2)

.70805*s*Units:-Unit(kg)^2*Units:-Unit(m)^2/(Units:-Unit(A)^2*Units:-Unit(s)^3)+144.*Units:-Unit(kg)^2*Units:-Unit(m)^2/(Units:-Unit(A)^2*Units:-Unit(s)^4)+0.3675e-4*s^2*Units:-Unit(kg)^2*Units:-Unit(m)^2/(Units:-Unit(A)^2*Units:-Unit(s)^2)

(9)

collect(.70805*s*Units[Unit](kg)^2*Units[Unit](m)^2/(Units[Unit](A)^2*Units[Unit](s)^3)+144.*Units[Unit](kg)^2*Units[Unit](m)^2/(Units[Unit](A)^2*Units[Unit](s)^4)+0.3675e-4*s^2*Units[Unit](kg)^2*Units[Unit](m)^2/(Units[Unit](A)^2*Units[Unit](s)^2), [Unit('A'), Unit('kg'), Unit('m'), Unit('s')])

(.70805*s/Units:-Unit(s)^3+144./Units:-Unit(s)^4+0.3675e-4*s^2/Units:-Unit(s)^2)*Units:-Unit(m)^2*Units:-Unit(kg)^2/Units:-Unit(A)^2

(10)

factor(.70805*s*Units[Unit](kg)^2*Units[Unit](m)^2/(Units[Unit](A)^2*Units[Unit](s)^3)+144.*Units[Unit](kg)^2*Units[Unit](m)^2/(Units[Unit](A)^2*Units[Unit](s)^4)+0.3675e-4*s^2*Units[Unit](kg)^2*Units[Unit](m)^2/(Units[Unit](A)^2*Units[Unit](s)^2))

0.3675e-4*Units:-Unit(kg)^2*Units:-Unit(m)^2*(19266.66666*s*Units:-Unit(s)+3918367.346+.9999999999*s^2*Units:-Unit(s)^2)/(Units:-Unit(A)^2*Units:-Unit(s)^4)

(11)

C2: Using the Natural Units Environment

In this environment, no Unit commands are required and the collection of units should work with Maple commands.
However, for the expressions discussed here, this would lead to a naming conflict with the complex variable s of the transfer function and the unit symbol s for seconds.

NULL

C3: A type declaration or unit assumptions on names

A type declaration as an option of commands like in

Units[TestDimensions](s*(.70805*s*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^2*Units[Unit](kg^2*m^2/(s^2*A^2))), {s::(Unit(1/s))})

true

(12)

could help Maple in simplification tasks (in its general meaning of making expressions shorter or smaller).
Alternatively, assumptions could provide information of which "unit type" a name is

`assuming`([simplify(H(s) = 60.*Units[Unit](m*kg/(s^2*A))/(.70805*s^2*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*s*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units[Unit](kg^2*m^2/(s^2*A^2))))], [s::(Unit(1/s))]); `assuming`([combine(H(s) = 60.*Units[Unit](m*kg/(s^2*A))/(.70805*s^2*Units[Unit](kg^2*m^2/(s^3*A^2))+144.*s*Units[Unit](kg^2*m^2/(s^4*A^2))+0.3675e-4*s^3*Units[Unit](kg^2*m^2/(s^2*A^2))), 'units')], [s::(Unit(1/s))])

Error, (in assuming) when calling 'property/ConvertProperty'. Received: 'Units:-Unit(1/s) is an invalid property'

 

On various occasions (beyond transfer functions) I have looked for such a functionality.

 

C4: DynamicSystems Package with units

C4.1: The complex variable s could be attributed to the unit 1/s (i.e. Hertz) either by default or as an option. This could enable using units within the dynamic system package which is not possible in Maple 2022. An example what the package provides currently can be found here: help(applications, amplifiergain)
The phase plot shows that the package is already implicitly assuming that the unit of s is Hertz. A logical extension would be to have magnitude plots with units (e.g. m/V, as in this example).

 

C4.2: A dedicated "gain" command that takes units into account and that could potentially simplify the transfer function to an expression like (2) in SI units. In such a way the transfer function is separated into a dimensionless (but frequency depended) term and a gain term with units.
This would make the transfer of transfer functions to MapleSim easy and avoid unit conversion errors.

 

Download Collecting_and_expanding_units.mw

A failing slinky is another intriguing physics phenome that can be easily reproduced with MapleSim.

The bottom of a vertically suspended slinky does not move when the top is released until the slinky is fully collapsed.

 

 

To model this realistically in MapleSim, it is necessary to

  • Establish a stretched equilibrium state at the start of the fall
  • Avoid penetration of windings when windings collapse (i.e. get into contact)

The equilibrium state is achieved with the snapshot option. Penetration is avoided with the Elasto Gap component. Details can be found in the attached model.

A good overview of “Slinky research” is given here. The paper provides a continuous description of the collapse process (using an inhomogenous wave equation combined with contact modeling!!!) and introduces a finite time for the collapse of all windings. Results for a slinky are presented that collapses after 0.27s. The attached model has sufficient fidelity to collapse at the same time.

Real Slinkies also feature a torsional wave that precedes the compression wave and disturbs an ideal collapse. This can be seen on slowmo footage and advanced computer models. With a torsion spring constant at hand (are there formulas for coil springs?), it could also be modeled with MapleSim.

Falling_slinky.msim

Maple allows to extract, manipulate, and optimize equations from a MapleSim model. Code can be generated from the equations in various programming languages. To verify the code, C code can be imported back into the original MapleSim model and compared to the model.

This verification step is not an everyday task, but it is advisable before the code is used elsewhere (e.g., in a controller). This post summarizes helpfull links and provides an additional example with equations that are too large to be efficiently verified by code review.

Comparison to a physical model is demonstrated here on an older version of MapleSim (~2015). In newer versions the import has changed (basics are described in Tutorial 6.6: Using the External C Code/DLL Custom Component App). An external C compiler must be set-up to make the import work.

The attached MapleSim model verifies against an optimized custom component. Instead of manually entering and modifying the code as described in the Tutorial 6.6, the model uses a Maple worksheet that programmatically generates C code from Maple equations and modifies the C code (sets C definitions and parameters) to be usable for MapleSim’s External C/Library Block App.

The Maple worksheet to generate and modify C code has been improved in many details with support from MaplePrime users for which I would like to express my thanks.

C_code_generation_of_optimised_code_for_MapleSim.mw

C_code_generation_of_optimised_code.msim

 

 

I have polished up findings with custom components to share it here:

Optimized code generated with Maple’s codegen package cannot be used in the same way as it was possible with older versions of MapleSim’s Custom Component Template.

Intermediate variables `tx` (where x is an integer) of the optimized code are interpreted as physical parameters in the current template version and not as variables. This makes sense and is more consistent with MapleSim’s definition of variables and parameters, but leads to errors in MapleSim.

The attached model shows how optimized code can be generated for the current template and compares an older, still working (!) template with the new one.

The attached worksheet contains commands to programmatically generate optimized code for the current Custom Component Template.

CustomComponentTemplates_comparision.msim

Optimized_code_for_custom_component_template.mw

1 2 3 4 5 6 Page 4 of 6