Skip to content

Commit 418deeb

Browse files
pabloselesonpabloseleson
committed
Updated HIP cylinder example
1 parent ba446a3 commit 418deeb

File tree

2 files changed

+117
-92
lines changed

2 files changed

+117
-92
lines changed

‎examples/thermomechanics/hip_cylinder.cpp

Lines changed: 110 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -36,15 +36,17 @@ void HIPCylinderExample( const std::string filename )
3636
// Material parameters
3737
// ====================================================
3838
double rho0 = inputs["density"];
39-
double K = inputs["bulk_modulus"];
40-
// double G = inputs["shear_modulus"]; // Only for LPS.
41-
double sc = inputs["critical_stretch"];
39+
double E = inputs["elastic_modulus"];
40+
double nu = 0.25; // Use bond-based model
41+
double K = E / ( 3 * ( 1 - 2 * nu ) );
42+
double G0 = inputs["fracture_energy"];
43+
double sigma_y = inputs["yield_stress"];
4244
double delta = inputs["horizon"];
4345
delta += 1e-10;
44-
// For PMB or LPS with influence_type == 1
45-
double G0 = 9 * K * delta * ( sc * sc ) / 5;
46-
// For LPS with influence_type == 0 (default)
47-
// double G0 = 15 * K * delta * ( sc * sc ) / 8;
46+
47+
// Problem parameters
48+
double sigma0 = inputs["traction"];
49+
double temp0 = inputs["reference_temperature"];
4850

4951
// ====================================================
5052
// Discretization
@@ -57,10 +59,10 @@ void HIPCylinderExample( const std::string filename )
5759
int halo_width = m + 1; // Just to be safe.
5860

5961
// ====================================================
60-
// Force model
62+
// Force model type
6163
// ====================================================
6264
using model_type = CabanaPD::PMB;
63-
CabanaPD::ForceModel force_model( model_type{}, delta, K, G0 );
65+
using mechanics_type = CabanaPD::ElasticPerfectlyPlastic;
6466

6567
// ====================================================
6668
// Custom particle generation and initialization
@@ -83,96 +85,120 @@ void HIPCylinderExample( const std::string filename )
8385
return true;
8486
};
8587

86-
// ====================================================
87-
// Simulation run with contact physics
88-
// ====================================================
89-
if ( inputs["use_contact"] )
88+
CabanaPD::Particles particles(
89+
memory_space{}, model_type{}, low_corner, high_corner, num_cells,
90+
halo_width, Cabana::InitUniform{}, init_op, exec_space{} );
91+
92+
auto rho = particles.sliceDensity();
93+
auto x = particles.sliceReferencePosition();
94+
double W = inputs["wall_thickness"];
95+
96+
auto init_functor = KOKKOS_LAMBDA( const int pid )
9097
{
91-
using contact_type = CabanaPD::NormalRepulsionModel;
92-
CabanaPD::Particles particles(
93-
memory_space{}, contact_type{}, low_corner, high_corner, num_cells,
94-
halo_width, Cabana::InitRandom{}, init_op, exec_space{} );
95-
96-
auto rho = particles.sliceDensity();
97-
auto x = particles.sliceReferencePosition();
98-
auto v = particles.sliceVelocity();
99-
auto f = particles.sliceForce();
100-
auto dx = particles.dx;
101-
102-
double vrmax = inputs["max_radial_velocity"];
103-
double vrmin = inputs["min_radial_velocity"];
104-
double vzmax = inputs["max_vertical_velocity"];
105-
double zmin = z_center - 0.5 * H;
106-
107-
auto init_functor = KOKKOS_LAMBDA( const int pid )
108-
{
109-
// Density
98+
double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) +
99+
( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center );
100+
if ( rsq > ( Rin + W ) * ( Rin + W ) &&
101+
rsq > ( Rout - W ) * ( Rout - W ) &&
102+
x( pid, 2 ) > low_corner[2] + W )
103+
{ // Powder density
104+
rho( pid ) = 0.7 * rho0;
105+
}
106+
else
107+
{ // Container density
110108
rho( pid ) = rho0;
111-
112-
// Velocity
113-
double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1;
114-
double vr = vrmax - vrmin * zfactor * zfactor;
115-
v( pid, 0 ) =
116-
vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
117-
v( pid, 1 ) =
118-
vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
119-
v( pid, 2 ) = vzmax * zfactor;
120109
};
121-
particles.updateParticles( exec_space{}, init_functor );
110+
};
111+
particles.updateParticles( exec_space{}, init_functor );
122112

123-
// Use contact radius and extension relative to particle spacing.
124-
double r_c = inputs["contact_horizon_factor"];
125-
double r_extend = inputs["contact_horizon_extend_factor"];
126-
// NOTE: dx/2 is when particles first touch.
127-
r_c *= dx[0] / 2.0;
128-
r_extend *= dx[0];
113+
// ====================================================
114+
// Boundary conditions planes
115+
// ====================================================
116+
CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane(
117+
low_corner[0], high_corner[0], low_corner[1], high_corner[1],
118+
low_corner[2], high_corner[2] );
129119

130-
contact_type contact_model( delta, r_c, r_extend, K );
120+
// ====================================================
121+
// Force model
122+
// ====================================================
123+
CabanaPD::ForceModel force_model( model_type{}, delta, K, G0 );
124+
// CabanaPD::ForceModel force_model( model_type{}, mechanics_type{}, rho,
125+
// delta, K, G0, sigma_y, rho0 );
131126

132-
CabanaPD::Solver solver( inputs, particles, force_model,
133-
contact_model );
134-
solver.init();
135-
solver.run();
136-
}
137127
// ====================================================
138-
// Simulation run without contact
128+
// Create solver
139129
// ====================================================
140-
else
141-
{
142-
CabanaPD::Particles particles(
143-
memory_space{}, model_type{}, low_corner, high_corner, num_cells,
144-
halo_width, Cabana::InitRandom{}, init_op, exec_space{} );
130+
CabanaPD::Solver solver( inputs, particles, force_model );
145131

146-
auto rho = particles.sliceDensity();
147-
auto x = particles.sliceReferencePosition();
148-
auto v = particles.sliceVelocity();
149-
auto f = particles.sliceForce();
132+
// ====================================================
133+
// Boundary conditions
134+
// ====================================================
150135

151-
double vrmax = inputs["max_radial_velocity"];
152-
double vrmin = inputs["min_radial_velocity"];
153-
double vzmax = inputs["max_vertical_velocity"];
154-
double zmin = z_center - 0.5 * H;
136+
// Create BC last to ensure ghost particles are included.
137+
double dx = particles.dx[0];
138+
// double dy = particles.dx[1];
139+
double dz = particles.dx[2];
140+
// double sigma0 = inputs["traction"];
141+
double b0 = sigma0 / dx;
142+
auto f = particles.sliceForce();
143+
x = particles.sliceReferencePosition();
144+
// Create an isostatic pressure BC.
145+
auto bc_op = KOKKOS_LAMBDA( const int pid, const double )
146+
{
147+
double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) +
148+
( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center );
149+
double theta =
150+
Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center );
155151

156-
auto init_functor = KOKKOS_LAMBDA( const int pid )
152+
// BC on outer boundary
153+
if ( rsq > ( Rout - dx ) * ( Rout - dx ) )
157154
{
158-
// Density
159-
rho( pid ) = rho0;
155+
f( pid, 0 ) += -b0 * Kokkos::cos( theta );
156+
f( pid, 1 ) += -b0 * Kokkos::sin( theta );
157+
}
158+
// BC on inner boundary
159+
else if ( rsq < ( Rin + dx ) * ( Rin + dx ) )
160+
{
161+
f( pid, 0 ) += b0 * Kokkos::cos( theta );
162+
f( pid, 1 ) += b0 * Kokkos::sin( theta );
163+
};
160164

161-
// Velocity
162-
double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1;
163-
double vr = vrmax - vrmin * zfactor * zfactor;
164-
v( pid, 0 ) =
165-
vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
166-
v( pid, 1 ) =
167-
vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
168-
v( pid, 2 ) = vzmax * zfactor;
165+
// BC on top boundary
166+
if ( x( pid, 2 ) > z_center + 0.5 * H - dz )
167+
{
168+
f( pid, 2 ) += -b0;
169+
}
170+
// BC on bottom boundary
171+
else if ( x( pid, 2 ) < z_center - 0.5 * H + dz )
172+
{
173+
f( pid, 2 ) += b0;
169174
};
170-
particles.updateParticles( exec_space{}, init_functor );
175+
};
176+
auto bc =
177+
createBoundaryCondition( bc_op, exec_space{}, particles, true, plane );
178+
179+
// ====================================================
180+
// Imposed field
181+
// ====================================================
182+
/*
183+
x = particles.sliceReferencePosition();
184+
auto temp = particles.sliceTemperature();
185+
const double low_corner_y = low_corner[1];
186+
// This is purposely delayed until after solver init so that ghosted
187+
// particles are correctly taken into account for lambda capture here.
188+
auto temp_func = KOKKOS_LAMBDA( const int pid, const double t )
189+
{
190+
temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t;
191+
};
192+
CabanaPD::BodyTerm body_term( temp_func, particles.size(), false );
193+
*/
171194

172-
CabanaPD::Solver solver( inputs, particles, force_model );
173-
solver.init();
174-
solver.run();
175-
}
195+
// ====================================================
196+
// Simulation run
197+
// ====================================================
198+
solver.init( bc );
199+
solver.run( bc );
200+
// solver.init( body_term );
201+
// solver.run( body_term );
176202
}
177203

178204
// Initialize MPI+Kokkos.

‎examples/thermomechanics/inputs/hip_cylinder.json

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,13 @@
11
{
22
"num_cells" : {"value": [100, 100, 140]},
33
"system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"},
4-
"density" : {"value": 7800, "unit": "kg/m^3"},
5-
"bulk_modulus" : {"value": 130e+9, "unit": "Pa"},
6-
"shear_modulus" : {"value": 78e+9, "unit": "Pa"},
7-
"critical_stretch" : {"value": 0.02},
8-
"horizon" : {"value": 0.003, "unit": "m"},
9-
"use_contact" : {"value": false},
10-
"contact_horizon_factor" : {"value": 0.9},
11-
"contact_horizon_extend_factor" : {"value": 0.01},
4+
"density" : {"value": 1200, "unit": "kg/m^3"},
5+
"elastic_modulus" : {"value": 3e+9, "unit": "Pa"},
6+
"fracture_energy" : {"value": 3000, "unit": "J/m^2"},
7+
"horizon" : {"value": 0.003, "unit": "m"},
8+
"yield_stress" : {"value": 30e6, "unit": "Pa"},
9+
"traction" : {"value": 2e6, "unit": "Pa"},
10+
"reference_temperature" : {"value": 300.0, "unit": "oC"},
1211
"cylinder_outer_radius_old" : {"value": 0.025, "unit": "m"},
1312
"cylinder_inner_radius_old" : {"value": 0.02, "unit": "m"},
1413
"cylinder_height_old" : {"value": 0.1, "unit": "m"},

0 commit comments

Comments
 (0)