Skip to content

Commit 2742c05

Browse files
pabloselesonstreeve
pabloseleson
authored and
streeve
committed
Updated HIP cylinder example;
new model form
1 parent b0fbe08 commit 2742c05

File tree

2 files changed

+112
-34
lines changed

2 files changed

+112
-34
lines changed

‎examples/thermomechanics/hip_cylinder.cpp

Lines changed: 105 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -36,18 +36,20 @@ 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;
4446
// For PMB or LPS with influence_type == 1
45-
double G0 = 9 * K * delta * ( sc * sc ) / 5;
47+
// double G0 = 9 * K * delta * ( sc * sc ) / 5;
4648
// For LPS with influence_type == 0 (default)
4749
// double G0 = 15 * K * delta * ( sc * sc ) / 8;
48-
double sigma_y = 10.0;
4950
double alpha = inputs["thermal_expansion_coeff"];
5051
double temp0 = inputs["reference_temperature"];
52+
double sigma0 = inputs["traction"];
5153

5254
// ====================================================
5355
// Discretization
@@ -60,7 +62,7 @@ void HIPCylinderExample( const std::string filename )
6062
int halo_width = m + 1; // Just to be safe.
6163

6264
// ====================================================
63-
// Force model
65+
// Force model type
6466
// ====================================================
6567
using model_type = CabanaPD::PMB;
6668
using mechanics_type = CabanaPD::ElasticPerfectlyPlastic;
@@ -94,39 +96,116 @@ void HIPCylinderExample( const std::string filename )
9496

9597
auto rho = particles.sliceDensity();
9698
auto x = particles.sliceReferencePosition();
97-
auto v = particles.sliceVelocity();
98-
auto f = particles.sliceForce();
99-
100-
double vrmax = inputs["max_radial_velocity"];
101-
double vrmin = inputs["min_radial_velocity"];
102-
double vzmax = inputs["max_vertical_velocity"];
103-
double zmin = z_center - 0.5 * H;
99+
double W = inputs["wall_thickness"];
104100

105101
auto init_functor = KOKKOS_LAMBDA( const int pid )
106102
{
107-
// Density
108-
rho( pid ) = rho0;
109-
110-
// Velocity
111-
double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1;
112-
double vr = vrmax - vrmin * zfactor * zfactor;
113-
v( pid, 0 ) =
114-
vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
115-
v( pid, 1 ) =
116-
vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
117-
v( pid, 2 ) = vzmax * zfactor;
103+
double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) +
104+
( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center );
105+
if ( rsq > ( Rin + W ) * ( Rin + W ) &&
106+
rsq > ( Rout - W ) * ( Rout - W ) &&
107+
x( pid, 2 ) > low_corner[2] + W )
108+
{ // Powder density
109+
rho( pid ) = 0.7 * rho0;
110+
}
111+
else
112+
{ // Container density
113+
rho( pid ) = rho0;
114+
};
118115
};
119116
particles.updateParticles( exec_space{}, init_functor );
120117

118+
// ====================================================
119+
// Boundary conditions planes
120+
// ====================================================
121+
CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane(
122+
low_corner[0], high_corner[0], low_corner[1], high_corner[1],
123+
low_corner[2], high_corner[2] );
124+
125+
// ====================================================
126+
// Force model
127+
// ====================================================
121128
rho = particles.sliceDensity();
122129
auto temp = particles.sliceTemperature();
123130
CabanaPD::ForceDensityModel force_model( model_type{}, mechanics_type{},
124131
rho, delta, K, G0, sigma_y, rho0,
125132
temp, alpha, temp0 );
126133

134+
// ====================================================
135+
// Create solver
136+
// ====================================================
127137
CabanaPD::Solver solver( inputs, particles, force_model );
128-
solver.init();
129-
solver.run();
138+
139+
// ====================================================
140+
// Boundary conditions
141+
// ====================================================
142+
143+
// Create BC last to ensure ghost particles are included.
144+
double dx = particles.dx[0];
145+
// double dy = particles.dx[1];
146+
double dz = particles.dx[2];
147+
// double sigma0 = inputs["traction"];
148+
double b0 = sigma0 / dx;
149+
auto f = particles.sliceForce();
150+
x = particles.sliceReferencePosition();
151+
// Create an isostatic pressure BC.
152+
auto bc_op = KOKKOS_LAMBDA( const int pid, const double )
153+
{
154+
double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) +
155+
( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center );
156+
double theta =
157+
Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center );
158+
159+
// BC on outer boundary
160+
if ( rsq > ( Rout - dx ) * ( Rout - dx ) )
161+
{
162+
f( pid, 0 ) += -b0 * Kokkos::cos( theta );
163+
f( pid, 1 ) += -b0 * Kokkos::sin( theta );
164+
}
165+
// BC on inner boundary
166+
else if ( rsq < ( Rin + dx ) * ( Rin + dx ) )
167+
{
168+
f( pid, 0 ) += b0 * Kokkos::cos( theta );
169+
f( pid, 1 ) += b0 * Kokkos::sin( theta );
170+
};
171+
172+
// BC on top boundary
173+
if ( x( pid, 2 ) > z_center + 0.5 * H - dz )
174+
{
175+
f( pid, 2 ) += -b0;
176+
}
177+
// BC on bottom boundary
178+
else if ( x( pid, 2 ) < z_center - 0.5 * H + dz )
179+
{
180+
f( pid, 2 ) += b0;
181+
};
182+
};
183+
auto bc =
184+
createBoundaryCondition( bc_op, exec_space{}, particles, true, plane );
185+
186+
// ====================================================
187+
// Imposed field
188+
// ====================================================
189+
/*
190+
x = particles.sliceReferencePosition();
191+
auto temp = particles.sliceTemperature();
192+
const double low_corner_y = low_corner[1];
193+
// This is purposely delayed until after solver init so that ghosted
194+
// particles are correctly taken into account for lambda capture here.
195+
auto temp_func = KOKKOS_LAMBDA( const int pid, const double t )
196+
{
197+
temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t;
198+
};
199+
CabanaPD::BodyTerm body_term( temp_func, particles.size(), false );
200+
*/
201+
202+
// ====================================================
203+
// Simulation run
204+
// ====================================================
205+
solver.init( bc );
206+
solver.run( bc );
207+
// solver.init( body_term );
208+
// solver.run( body_term );
130209
}
131210

132211
// 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)