[Study] Shell-to-solid coupling

[WeiyiVirtonomy]

This PR aims to couple a shell body and a solid body (tie-contact). The coupling is achieved by the consistent mapping of velocity and conservative mapping of force on an overlapping layer.
For details, please check the README file: tests/3d_examples/test_3d_solid_to_shell_coupling/doc/README.md
An example test_3d_solid_to_shell_coupling is used to prove the feasibility of the method. However, the validation is not finished yet due to high computational cost, and will be left for future work.

[WeiyiVirtonomy]

I'll compare it with the fully volumetric solid results.

Example: a solid plate is reinforced by two shell plates on top and bottom surface. The stiffness ratio is 0.1.
The 1st load case is extension, while the 2nd load case is shear.
The height of the solid plate is 20mm, while the thickness of shell is selected as 4mm for now.
For the fully volumetric solid test, a composite material is used.

With a resolution of thickness/4.0:

• Extension case
  Left: solid, right: solid-shell coupling

Grey: solid-shell coupling

The position discrepancy is larger near the coupled surfaces:

• Shear case

• Stability
  I started from dp = thickness, but the coupled case soon failed. The stability of the coupling algorithm and the resolution needs further investigation.

• Discussion
  For the next step, I think there are several directions:

1. Replace the lattice particle generator with the relaxed particles and identify the surface layers as shell
2. Test on different stiffness ratios
3. Test on different shell thickness / length ratios
4. Use v_mid+omegaxr_rel instead of the mid-surface velocity as the constraint velocity
5. Discuss how to deal with the normal direction at the corners

I hope we can discuss the priority of these problems tomorrow.

[WeiyiVirtonomy]

I tried to compute the coupling force with ghost shell particles, but it doesn't work. I'm not sure if I understood your idea correctly, so I listed the procedure here:

1. Solid and shell are modelled as bodies without overlapping

2. For solid, the shell particles and ghost particles are considered a part of the solid body, which means that:
   2.1 The linear correction matrix includes the contribution of shell & dummy particles
   2.2 The corrected PK1 stress of shell & dummy particles needs to be included in 1st integration
   2.3 The velocity of shell & dummy particles needs to be included in 2nd integration

3. For shell, the coupling force is computed as the contact force acting on solid particles

The contribution of dummy shell particles is added to the contact neighbors using the fluid-shell and shell-fluid contact relations, which implies that the properties (e.g., vel and stress) of ghost particles are the same as the mid-surface particles.

Based on this assumption, 2.3 is easy to implement by:

Matd deformation_gradient_change_rate = Matd::Zero();
        for (size_t k = 0; k < contact_configuration_.size(); ++k)
        {
            const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
            const auto *Vol_k = contact_Vol_[k];
            const auto *vel_k = contact_vel_[k];
            for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
            {
                size_t index_j = contact_neighborhood.j_[n];

                SPH::Vecd gradW_ij = contact_neighborhood.dW_ij_[n] * Vol_k[index_j] * contact_neighborhood.e_ij_[n];
                deformation_gradient_change_rate -= (vel_[index_i] - vel_k[index_j]) * gradW_ij.transpose();
            }
        }
        deformation_gradient_change_rate *= B_[index_i];

The velocity of the contact shell particles are just the real velocity of the shell, which is equivalent to applying a velocity boundary condition to the solid.

The problem is how to calculate the corrected PK1 stress of the shell particles. The method I'm using now assumes that the deformation tensor of a shell particle j used in the 1st integration of solid is the same as the real deformation tensor of the shell. The linear correction matrix is also the one of the single shell. The material is one of shell. The strain and total stress corrections are neglected. Hence, the PK1 stress of shell i is computed by:

Matd F_i = F_[index_i] + dF_dt_[index_i] * dt * 0.5;
Matd stress_PK1_B_i = elastic_solid_.StressPK1(F_i, index_i) * B_[index_i].transpose();
// transform the PK1 stress in initial local coordinate to the global coordinate
stress_PK1_B_[index_i] = transformation_matrix0_[index_i].transpose() * stress_PK1_B_i * transformation_matrix0_[index_i];

I didn't consider the numerical damping when computing the contribution of shells:

Vecd force = Vecd::Zero();
        for (size_t k = 0; k < contact_configuration_.size(); ++k)
        {
            const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
            const Matd *stress_PK1_B_k = contact_stress_PK1_B_[k];
            const auto *Vol_k = contact_Vol_[k];
            for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
            {
                size_t index_j = contact_neighborhood.j_[n];
                Vecd e_ij = contact_neighborhood.e_ij_[n];
                force += mass_[index_i] * inv_rho0_ * contact_neighborhood.dW_ij_[n] * Vol_k[index_j] *
                         (stress_PK1_B_[index_i] + stress_PK1_B_k[index_j]) * e_ij;
            }
        }

The external forces on the shell are computed in the same way:

force += mass_k[index_j] * inv_rho0_k * contact_neighborhood.dW_ij_[n] * Vol_[index_i] *
                         (stress_PK1_B_[index_i] + stress_PK1_B_k[index_j]) * e_ij;

However, the solid exploded after several steps:

[WeiyiVirtonomy]

I tried to compute the coupling force with ghost shell particles, but it doesn't work. I'm not sure if I understood your idea correctly, so I listed the procedure here:

1. Solid and shell are modelled as bodies without overlapping

2. For solid, the shell particles and ghost particles are considered a part of the solid body, which means that: 2.1 The linear correction matrix includes the contribution of shell & dummy particles 2.2 The corrected PK1 stress of shell & dummy particles needs to be included in 1st integration 2.3 The velocity of shell & dummy particles needs to be included in 2nd integration 3. For shell, the coupling force is computed as the contact force acting on solid particles

The contribution of dummy shell particles is added to the contact neighbors using the fluid-shell and shell-fluid contact relations, which implies that the properties (e.g., vel and stress) of ghost particles are the same as the mid-surface particles.

Based on this assumption, 2.3 is easy to implement by:

Matd deformation_gradient_change_rate = Matd::Zero();
        for (size_t k = 0; k < contact_configuration_.size(); ++k)
        {
            const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
            const auto *Vol_k = contact_Vol_[k];
            const auto *vel_k = contact_vel_[k];
            for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
            {
                size_t index_j = contact_neighborhood.j_[n];

                SPH::Vecd gradW_ij = contact_neighborhood.dW_ij_[n] * Vol_k[index_j] * contact_neighborhood.e_ij_[n];
                deformation_gradient_change_rate -= (vel_[index_i] - vel_k[index_j]) * gradW_ij.transpose();
            }
        }
        deformation_gradient_change_rate *= B_[index_i];

The velocity of the contact shell particles are just the real velocity of the shell, which is equivalent to applying a velocity boundary condition to the solid.

The problem is how to calculate the corrected PK1 stress of the shell particles. The method I'm using now assumes that the deformation tensor of a shell particle j used in the 1st integration of solid is the same as the real deformation tensor of the shell. The linear correction matrix is also the one of the single shell. The material is one of shell. The strain and total stress corrections are neglected. Hence, the PK1 stress of shell i is computed by:

Matd F_i = F_[index_i] + dF_dt_[index_i] * dt * 0.5;
Matd stress_PK1_B_i = elastic_solid_.StressPK1(F_i, index_i) * B_[index_i].transpose();
// transform the PK1 stress in initial local coordinate to the global coordinate
stress_PK1_B_[index_i] = transformation_matrix0_[index_i].transpose() * stress_PK1_B_i * transformation_matrix0_[index_i];

I didn't consider the numerical damping when computing the contribution of shells:

Vecd force = Vecd::Zero();
        for (size_t k = 0; k < contact_configuration_.size(); ++k)
        {
            const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
            const Matd *stress_PK1_B_k = contact_stress_PK1_B_[k];
            const auto *Vol_k = contact_Vol_[k];
            for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
            {
                size_t index_j = contact_neighborhood.j_[n];
                Vecd e_ij = contact_neighborhood.e_ij_[n];
                force += mass_[index_i] * inv_rho0_ * contact_neighborhood.dW_ij_[n] * Vol_k[index_j] *
                         (stress_PK1_B_[index_i] + stress_PK1_B_k[index_j]) * e_ij;
            }
        }

The external forces on the shell are computed in the same way:

force += mass_k[index_j] * inv_rho0_k * contact_neighborhood.dW_ij_[n] * Vol_[index_i] *
                         (stress_PK1_B_[index_i] + stress_PK1_B_k[index_j]) * e_ij;

However, the solid exploded after several steps:

Sorry, I should use the material of the solid instead of the shell.

It can run to the end, but the solid particles are clustered with larger shell thickness

[WeiyiVirtonomy]

Using the internal force of coupled solid particles as the external force of shell looks good so far:

[WeiyiVirtonomy]

Shall I move the coupling classes to the source file of extra_source_and_tests?

[Xiangyu-Hu]

Shall I move the coupling classes to the source file of extra_source_and_tests?

I think that you can put it into the main library.

[Xiangyu-Hu]

@WeiyiVirtonomy is this pull request is ready for review?

[Xiangyu-Hu]

@WeiyiVirtonomy is this ready for a real pull request?

[WeiyiVirtonomy]

@WeiyiVirtonomy is this ready for a real pull request?

I'm still waiting for the solid-discretized results for comparison. I'll mark it as ready for review after adding the validation results to README file

[Xiangyu-Hu]

@WeiyiVirtonomy could it possible using a distributor between soft and stiff material so that it works for arbitrary shell solid pair?

[WeiyiVirtonomy]

@WeiyiVirtonomy could it possible using a distributor between soft and stiff material so that it works for arbitrary shell solid pair?

Could you explain a bit more what the distributor between different material means here?
Now, the algorithm can deal with solid and shell with different materials (soft shell-stiff solid or soft solid-stiff shell), although the limitation is not tested yet

[WeiyiVirtonomy]

@WeiyiVirtonomy is this ready for a real pull request?

The validation case with solid discretization still hasn't finished. I don't want to hold this PR further, so I have added comments in the README and code and leave the validation as future work.
If this makes sense to you, then this PR is ready for review