Elastic response to mesh modification

[AdrienWehrle]

Hi everyone,

I have been working on an application based on MASTODON to study the elastic response of a domain on an inclined surface, after the removal of some support at its front.

What I did first was to apply a Pressure boundary condition, with the frontal_pressure function depending on time and taking two different values before and after a given time t (simulating the loss of support in front of the domain).

[BCs]
  [Pressure]
    [front_pressure]  
    boundary = front
    function = frontal_pressure
    displacements = 'disp_x disp_y disp_z'
    []
  []

This worked great and gave interesting results. But this only represents the result of a material removal, not the actual removal itself. Next step was therefore to actively remove some material at the front of the domain and get the pressure release as a result, instead of a preset. Just as it actually happens in a real case.

Based on #20601, I looked into two potential solutions:

1. https://mooseframework.inl.gov/source/userobject/CoupledVarThresholdElementSubdomainModifier.html
2. https://mooseframework.inl.gov/source/meshgenerators/BlockDeletionGenerator.html

My progress (and issues) on those two options:

1. This is on the fly. One input file, and a moving boundary. However, I need to assign the newly deactivated block (considered as removed, chopped off) to an existing subdomain.
   And that is where I got stuck... How to assign it to a subdomain that doesn't have any, literally any, impact on the main domain anymore?
   Here is how far I got into that first direction with my input file: chop_dynamic.i

2. This is the restart option. Two input files, the first one for initialization, and the second one where a new mesh is created with the material chopped off and the solutions remapped based on the first simulation. It also seems like a good idea, because in my case I remove material only once and this is irreversible.
   I think I got pretty far here, but somehow I don't get significant elastic response (only a very slow displacement in x) in my transient simulation although I removed 500m off my domain and initialized with displacements corresponding to the entire loaded domain... So I should have a significant extension-dominated reaction...! Also the solve is very slow, much slower than with a pressure release of similar magnitude, somehow suggesting something is wrong...
   Here is how far I got into that second direction with my input files: chop_init.i (initialization) and chop_restart.i (transient after removal)

Any help/thoughts on one of these options or both would be much appreciated! Thank you a lot in advance!

[GiudGiud]

Hello

For option 1, if you move elements to a domain that is outside the block restrictions of the kernels/variables, then it is effectively removed from the simulation.
There might also be material property values you could set in that domain that negate all form of displacement (like incredibly stiff). This might not be a good idea as I think you may pay a steep price in convergence properties.

For your inputs in 1 and 2, I would set displacement = 'disp_x disp_y ...' in GlobalParams. it s too easy to have forgotten it somewhere.

Guillaume

[GiudGiud]

oh did you block restrict the variables outside of the elements with no physics? This is actually important, otherwise the matrix has a zero for the qps there.

The block restriction of the variables is usually inherited by the kernels btw

[AdrienWehrle]

oh did you block restrict the variables outside of the elements with no physics?

I am not sure I completely understand. Variables are not block-restricted at the moment, but they probably should? Should I use inactive on block 255 (inactive element) for all variables?

[GiudGiud]

They should be. Pass the block = in the variables blocks (each of them) not the kernels.

inactive is not for this. inactive is to turn off an object in the input file, like entirely not just block restrict

[AdrienWehrle]

Thank you a lot, that indeed solved the zero qps issue!

Only remaining problem is I restrict the physics to block 1 and 2 (main domain and chopped area, respectively), then change the block ID of the chopped area from 2 to 255 (to deactivate it). Therefore block 2 doesn't exist anymore, but is passed in block = in variable blocks, resulting in ExodusII_IO_Helper: block id 2 not found in block_ids. at chop time (t=1).
I think one workaround is to have not only the one isolated element (block_id=255), but also a second one with block_id=2 (useless but containing the physics) so that block id 2 doesn't disappear. Do you think that is a good idea?

Thank you a lot for your help!

[GiudGiud]

Yeah I think you ll have to do that.
Or you start with just all elements in block 1, then move some of the chopped elements in block 1 to block 255?

[AdrienWehrle]

Or you start with just all elements in block 1, then move some of the chopped elements in block 1 to block 255?

Indeed, I managed to trade my initial SubdomainBoundingBoxGenerator + CoupledVarThresholdElementSubdomainModifier with threshold on time, with one single CoupledVarThresholdElementSubdomainModifier on a variable depending on space and time! Thank you so much for your help on that!

[AdrienWehrle]

One last thing that I can't get right for this problem (sorry to disturb you with that!):
I added the following function:

[Functions]
  [chop_criterion]
    type = ParsedFunction
    value = 'if((x>29300.)&(t>1.), 1., 0.)'
  []
[]

that I then apply in my subdomain modifier:

[UserObjects]
  [calving_event]
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = 'chop_boolean'
    block = 1
    criterion_type = ABOVE
    threshold = 0.5
    subdomain_id = 255
    moving_boundary_name = downstream 
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

chop_boolean is the variable set by chop_criterion.

I thought I solved my issue, but I get a weird chop where not all expected elements are moved to the 255 subdomain (see screenshots below) although the condition is a clear cut and not at the border of an element row (x=29000 and x=29500 are the two element rows before and after the condition x=29300)... As you can see on the sfcreenshots, the outer band of elements is entirely removed but not issue appears in the second one....

This is independent of the threshold, I always get such elements in and out no matter what value I take...

I was earlier using the same spatial condition with SubdomainBoundingBoxGenerator and I got a nice and clear cut as expected... Any ideas on that? Thank you a lot

[GiudGiud]

I think the centroid of the element is considered. If you move this x=29300 closer to the next layer (3rd) does it cut properly?

[GiudGiud]

actually no scratch that.
This is how the element value is determined:

  Real avg_val = 0;

  for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
    avg_val += _v[qp] * _JxW[qp] * _coord[qp];
  avg_val /= _current_elem_volume;

  return avg_val;

sp you have to make that average below the threshold value.
Either smooth the transition with a linear function in between the jump, or jump to -1 instead of 0, or raise/lower the threshold appropriately

[AdrienWehrle]

Oh woh, I wasn't looking in the right direction, I also initially thought it was about the spatial threshold.

What I have for the moment is my chop_criterion taking -1000 or 1000 and threshold at 0 so there can't be any overrunning.

Thank you a lot for all your help on that! 👍

[AdrienWehrle]

Next step will be to try and have my chop_criterion to depend on a damage value I compute from stresses. Challenge is my damage is a material property and not an (aux)variable, therefore I can't access it directly with a CoupledVarThresholdElementSubdomainModifier... I tried to find a way to convert a material property to a variable, but that sounds like a bad idea. Would you by any chance have an advice on that before I mark this question as answered? Thank you a lot!

The damage implementation is visible here.

[GiudGiud]

Hello

it s fine to convert material property to auxiliary variable in your case. You are only interested about the value.
it s not fine when you care about transferring derivatives.

To do that use MaterialRealAux or one of the other template instantiations depending on the type of your material property

Guillaume