OneEndFixedBC didn't work for applying a external torque about pyelastica HOT 4 CLOSED

Hi @Junang-Wang

It seems the problem not the boundary condition but the selected parameters might be causing a stiff and unstable simulation. Have you consider trying out the example parameters here and changing one by one to your parameters? You can debug this way.

from pyelastica.

The following is my code snippet:

#--------simulator----------
class UnitMagneticRodSimulator(BaseSystemCollection, Constraints, Forcing, Damping, Connections,CallBacks):
pass
Uniform_M_Sim = UnitMagneticRodSimulator()
#--------mmGS unit----------
density = 2.273 #mg/mm^3
base_length = 3 #mm
base_radius = 0.3 #mm
scale_E = 1e-4 #scale dowing Young's modulus
E = 1.4e9 * scale_E #Young's modulus
shear_modulus = E/3 * scale_E # shear modulus

dt = 1e-4 # time step
nu = 5
endtime = 5

#--------rod definition-------
n_elem = 20
start = np.array([0.0,0.0,0.0])
direction = np.array([0.0,1.0,0.0])
normal = np.array([1.0,0.0,0.0])
M_rod = CosseratRod.straight_rod(n_elem, start, direction, normal, base_length, base_radius, density, youngs_modulus=E, shear_modulus=shear_modulus)

Uniform_M_Sim.append(M_rod)
#--------magnetic properties-----
magnetization_density = 1.28e5 * 1e-3 #A/mm
magnetic_field_angle = np.pi/2
magnetic_field = 90e3scale_E # mg/(s^2A)

magnetization_direction = np.ones((n_elem)) * direction.reshape(3, 1)

#------set the constant magnetic field object-----
magnetic_field_amplitude = magnetic_field* np.array([0.0, np.cos(magnetic_field_angle), np.sin(magnetic_field_angle)])
magnetic_field_object = ConstantMagneticField(
magnetic_field_amplitude, ramp_interval = 0.5, start_time = 0.5, end_time = endtime
)

#--------constrain----------
Uniform_M_Sim.constrain(M_rod).using(
OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
)

#--------damping------------
Uniform_M_Sim.dampen(M_rod).using(
AnalyticalLinearDamper, damping_constant = nu, time_step = dt
)
#--------force--------------
Uniform_M_Sim.add_forcing_to(M_rod).using(
MagneticForces,
external_magnetic_field = magnetic_field_object,
magnetization_density = magnetization_density,
magnetization_direction = magnetization_direction,
rod_volume = M_rod.volume,
rod_director_collection = M_rod.director_collection.copy(),
)
#------callback function------
class MagneticRodCallBack(CallBackBaseClass):
def init(self, step_skip:int, callback_params:dict):
super().init()
self.step_skip = step_skip
self.callback_params = callback_params

def make_callback(self, system, time, current_step: int):
    if current_step % self.step_skip == 0:
        self.callback_params["time"].append(time)
        self.callback_params["position"].append(system.position_collection.copy())
        self.callback_params["velocity"].append(system.velocity_collection.copy())
        return

MR_list = defaultdict(list)
Uniform_M_Sim.collect_diagnostics(M_rod).using(
MagneticRodCallBack, step_skip = 100, callback_params = MR_list
)

Uniform_M_Sim.finalize()
#--------time integration-----
timestepper = PositionVerlet()
final_time =endtime
total_steps = int(final_time / dt)
integrate(timestepper, Uniform_M_Sim, final_time, total_steps)

from pyelastica.

Hi @armantekinalp, it seems that I have chosen an unreasonable shear modulus. Thank you for your advise.

from pyelastica.

No problem.

from pyelastica.