Large flunctuations of velocities on the fixed boundary

[lyyc199586]

Hi everyone,

I am conducting large deformation elastodynamics simulations involving a compressive wave propagating from the left end of a bar. At the left boundary, I apply a displacement boundary condition that ramps up to a fixed displacement value, holding it constant afterward.

While the displacement profile appears correct, I am observing significant fluctuations in the velocity at the left boundary. I expected these velocities to remain relatively small due to the fixed displacement condition following the ramp-up phase. Could these fluctuations be attributed to numerical artifacts, or are there other potential causes for this behavior?

Thank you for any insights.

animation of disp_x and vel_x:

disp_vel.mp4

my input file:

rho = 1e-6
hht_alpha = -0.00
beta = '${fparse (1-hht_alpha)^2/4}'
gamma = '${fparse 1/2-hht_alpha}'

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = './mesh/bar.msh'
  []
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  alpha = ${hht_alpha}
  beta = ${beta}
  gamma = ${gamma}
  large_kinematics = true
  use_displaced_mesh = false
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[AuxVariables] 
  [accel_x]
  []
  [vel_x]
  []
  [accel_y]
  []
  [vel_y]
  []
  [accel_z]
  []
  [vel_z]
  []
  [stress_xx]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = RankTwoAux
      rank_two_tensor = cauchy_stress
      index_i = 0
      index_j = 0
      execute_on = 'INITIAL TIMESTEP_END'
    []
  []
[]

[Kernels]
  [solid_x]
    type = TotalLagrangianStressDivergence
    variable = disp_x
    displacements = 'disp_x disp_y disp_z'
    component = 0
  []
  [inertia_x] # M*accel + eta*M*vel
    type = ADInertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    density = density
  []
  [solid_y]
    type = TotalLagrangianStressDivergence
    variable = disp_y
    displacements = 'disp_x disp_y disp_z'
    component = 1
  []
  [inertia_y] # M*accel + eta*M*vel
    type = ADInertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    density = density
  []
  [solid_z]
    type = TotalLagrangianStressDivergence
    variable = disp_z
    displacements = 'disp_x disp_y disp_z'
    component = 2
  []
  [inertia_z] # M*accel + eta*M*vel
    type = ADInertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    density = density
  []
[]

[AuxKernels]
  [accel_x] # Calculates and stores acceleration at the end of time step
    type = NewmarkAccelAux
    variable = accel_x
    displacement = disp_x
    velocity = vel_x
    execute_on = timestep_end
  []
  [vel_x] # Calculates and stores velocity at the end of the time step
    type = NewmarkVelAux
    variable = vel_x
    acceleration = accel_x
    execute_on = timestep_end
  []
  [accel_y] # Calculates and stores acceleration at the end of time step
    type = NewmarkAccelAux
    variable = accel_y
    displacement = disp_y
    velocity = vel_y
    execute_on = timestep_end
  []
  [vel_y] # Calculates and stores velocity at the end of the time step
    type = NewmarkVelAux
    variable = vel_y
    acceleration = accel_y
    execute_on = timestep_end
  []
  [accel_z] # Calculates and stores acceleration at the end of time step
    type = NewmarkAccelAux
    variable = accel_z
    displacement = disp_z
    velocity = vel_z
    execute_on = timestep_end
  []
  [vel_z] # Calculates and stores velocity at the end of the time step
    type = NewmarkVelAux
    variable = vel_z
    acceleration = accel_z
    execute_on = timestep_end
  []
[]

[Functions]
  [load]
    type = ADParsedFunction
    expression = '0.4*(1-exp(-2.3983e5*t))'
  []
[]

[BCs]
  [left_x]
    type = ADFunctionDirichletBC
    boundary = '1'
    variable = disp_x
    function = load
    preset = true
  []
[]

[Materials]
  [bulk]
    type = ADGenericConstantMaterial
    prop_names = 'density'
    prop_values = '${rho}'
  []
  [elasticity]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 30000
    poissons_ratio = 0.2
  []
  [stress]
    type = ComputeLagrangianLinearElasticStress
    # outputs = 'vtk'
  []
  [strain]
    type = ComputeLagrangianStrain 
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  start_time = 0
  end_time = 3e-4
  dt = 1e-6
  dtmin = 5e-9
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8
  l_tol = 1e-6
  # r_tol = 1e-6
  # nl_max_its = 20
  nl_max_its = 500
  l_max_its = 500
  line_search = None
  [TimeIntegrator]
    type = NewmarkBeta
    inactive_tsteps = 1
  []
  []
[]

[maxnezdyur]

This is very similar to this thread: #22574

I believe it would be fixed with ramping up your BC or applying an initial condition on your velocity.

[lyyc199586]

I will try with IC on velocity

[lyyc199586]

humm, the displacement boundary is indeed ramped up from 0 to maximum value using an exponential function, and I tried to set initial velocity to 0, but still see those fluctuations.

[maxnezdyur]

If there is no change to your problem definition, I believe you will have these artifacts.
Look at the velocity of your disp bc, if you take the derivative at 0, the velocity is close to 1. So if you imagine the acceleration needed to going from v= 0 to v =1 for any dt-> 0 you can see that the acceleration needed is unbounded. So a ramped BC for your displacement would look like some sigmoid function.

[lyyc199586]

ok, I get this part. But why the fluctuation of velocity keeps to exist after the displacement is almost fixed?

[maxnezdyur]

It's been a while since I looked at this problem, but when you do the math for a simple problem, I think you will be able to see why. Do the first few time steps manually for a problem with IC = 0 for displacement and velocity. The whole block then has a BC on displacement such that disp = t. When you write out these equations, you should see that how the acceleration and velocity is calculated the velocity will oscillate.