Channel flow with clamped elastic walls

[rubensamarojr]

Hello there,

Is possible to simulate a channel flow with clamped elastic walls using the current master version ? Similar to that shown in the video below:

https://www.youtube.com/watch?v=Bfhe6MGZ3eI

Also, is it possible to attach two elastic solids or model an elastic solid with two young's modulus ? I'd like to run a case similar to the draft below:

Until now, I've simulated only elastic solids attached to one fixed wall. Maybe if I change a specific class related to elastic solids in SPHinXsys I can do that.

Any help is very welcome.

[rubensamarojr]

Hi there,

After a more careful reading of the examples, I verified that the definition of clamped ends are quite simple. All we need is to define the shape of the constrain regions as MultiPolygonShape class.

The question I still have is: Is it possible to attach two elastic solids with different Young's modulus ?

Best

[FabienPean-Virtonomy]

Hi @rubensamarojr, AFAIK attachment/welding boundary condition is currently not implemented.

[Xiangyu-Hu]

@rubensamarojr @FabienPean-Virtonomy Hi, one of my PhD students, Yaru, is working on this. She is designing a special relationship, we called dual to handle two or more bodies with different materials welded together. However, we need wait some time for it ready.

[rubensamarojr]

Hello all,

Thank you @Xiangyu-Hu and @FabienPean-Virtonomy .

Glad to hear that this feature is under development.

For my specific case, the Young's modulus for the left region is E_left ~10^7, while for the right region is E_right ~10^4. I mean, the Young's modulus depends only on the spatial region. Do you think that if I add a hard-code condition (if pos < XX) in the class that calculated the constitutive equation, then I get this, at least thinking about the numerical code ? Maybe I can get some weird physical results regardless if the code is ok.

See the code below that I added to the elastic_solid.cpp source file. I'm using a linear elastic solid and here the limit position is x=1.0m. I'm not sure if it is enough :confused:.

OBS. Since sound_speed will be calculated based on the higher Young Modulus, E_left ~10^7, I assume the time step condition ensures elastic solid stability. Also, the horizontal displacement is not large, so that the spatial IF condition is OK for this case.

Matd LinearElasticSolid::ConstitutiveRelation(Matd &F, size_t particle_index_i)
    {
        Matd strain = 0.5 * (~F * F - Matd(1.0));
        Matd sigmaPK2 = lambda0_ * strain.trace() * Matd(1.0) + 2.0 * G0_ * strain;
        // New sigmaPK2 based on the position X
        Real posX = 1.0;
        if(pos_n_[particle_index_i][0] > posX) {
            // New Poisson ratio
            Real poisson_ratio_right = 0.3;
            // New Young's modulus
            Real youngs_modulus_right = 10.0e+4.0;
            // New lambda
            Real lambda0_right = nu_ * youngs_modulus_right / (1.0 + poisson_ratio_right) / (1.0 - 2.0 * poisson_ratio_right);
            // New bulk modulus
            Real G0_right = 0.5 * youngs_modulus_right / (1.0 + poisson_ratio_right);
            sigmaPK2 = lambda0_right * strain.trace() * Matd(1.0) + 2.0 * G0_right * strain;
        }
        return sigmaPK2;
    }

[FabienPean-Virtonomy]

Hey there, I think that should be a sufficient workaround for your scenario at the moment

[rubensamarojr]

Hi there,

Thanks for the reply.

Below, you can see the files where I made minor changes. The code is running well after the changes.

SPHINXsys/src/shared/materials/elastic_solid.h

class ElasticSolid : public Solid
...
virtual Matd ConstitutiveRelationPosX(Matd &deformation, size_t particle_index_i) = 0;

...

class LinearElasticSolid : public ElasticSolid
...
virtual Matd ConstitutiveRelationPosX(Matd &deformation, size_t particle_index_i) override;

SPHINXsys/src/shared/materials/elastic_solid.cpp

Matd LinearElasticSolid::ConstitutiveRelationPosX(Matd &F, size_t particle_index_i)
{
    Matd strain = 0.5 * (~F * F - Matd(1.0));
    // New Poisson
    Real poisson_ratio_posx = 0.3;
    // New Young's modulus
    Real youngs_modulus_posx = 1.3e7;
    // New lambda
    Real lambda0_posx = nu_ * youngs_modulus_posx / (1.0 + poisson_ratio_posx) / (1.0 - 2.0 * poisson_ratio_posx);
    // New bulk modulus
    Real G0_posx = 0.5 * youngs_modulus_posx/ (1.0 + poisson_ratio_posx);
    Matd sigmaPK2 = lambda0_posx * strain.trace() * Matd(1.0) + 2.0 * G0_posx * strain;

    return sigmaPK2;
}

SPHINXsys/src/shared/particle_dynamics/solid_dynamics/solid_dynamics.h

class BaseElasticRelaxation : public ParticleDynamics1Level, public ElasticSolidDataInner
...
StdLargeVec<Vecd> &pos_0_;

SPHINXsys/src/shared/particle_dynamics/solid_dynamics/solid_dynamics.cpp

BaseElasticRelaxation::BaseElasticRelaxation(BaseBodyRelationInner &inner_relation)
    : ParticleDynamics1Level(*inner_relation.sph_body_),
    ElasticSolidDataInner(inner_relation), Vol_(particles_->Vol_),
    rho_n_(particles_->rho_n_), mass_(particles_->mass_),
    pos_n_(particles_->pos_n_), vel_n_(particles_->vel_n_), dvel_dt_(particles_->dvel_dt_),
    pos_0_(particles_->pos_0_),
    B_(particles_->B_), F_(particles_->F_), dF_dt_(particles_->dF_dt_) {}

...

void StressRelaxationFirstHalf::Initialization(size_t index_i, Real dt)
    ...
    // New sigmaPK2 based on the position X
    Real posX = 1.85e-2;
    if(pos_0_[index_i][0] < posX) {
        stress_PK1_[index_i] = F_[index_i] * material_->ConstitutiveRelationPosX(F_[index_i], index_i) * B_[index_i];
    }

...

void KirchhoffStressRelaxationFirstHalf::Initialization(size_t index_i, Real dt)
    ...
    // New sigmaPK2 based on the position X
    Real posX = 1.85e-2;
    if(pos_0_[index_i][0] < posX) {
        stress_PK1_[index_i] = F_[index_i] * material_->ConstitutiveRelationPosX(F_[index_i], index_i);
    }

[Xiangyu-Hu]

Yes. Such simple changes work for your problem.
The limit is that it works only for material properties with small change or smooth variations and need to choose parameters in runtime. Actually, the code is able to do this as those materials classes with "local" in their names. In these classes, the parameters are defined when the material is constructed and there is no need to choose them in runtime. The difficult is to generalized the application for materials with very different properties, such as rubber glued on steel surface. In this case, simple solution with local properties does not work as it leads numerical instabilities. That is the reason why we still need one student specifically work on it.

[rubensamarojr]

Thank you for the detailed explanation Prof. Hu. I understood the challenges related to this topic. Looking forward to the next improvements.

[Xiangyu-Hu]

Some new results from my students suggest simply use different constitute relation in one body is already can handle the cases with quite different stiffness or density. We will soon have test cases in the repository on this.

[Xiangyu-Hu]

Now, the new functionality of composite solid is done.