Duncan-Selig soil model #195
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## dev #195 +/- ##
================================================
- Coverage 90.7760% 90.7600% -0.0160%
================================================
Files 661 662 +1
Lines 29087 29253 +166
Branches 3314 3330 +16
================================================
+ Hits 26404 26550 +146
- Misses 2683 2703 +20 ☔ View full report in Codecov by Sentry. |
| @@ -0,0 +1,39 @@ | |||
| # A TEST MODEL FOR DUNCAN-SELIG MATERIAL | |||
There was a problem hiding this comment.
This is a two-element uniaxial compression test.
| node 5 1 1 | ||
| node 6 1 2 | ||
|
|
||
| material DuncanSelig 1 14.7 6000 .6 4000 .2 .7 .1 .7 .5 |
There was a problem hiding this comment.
Nine parameters are required to define the model, see below for details.
| double p_atm = 14.7; | ||
| double ref_elastic = 400. * p_atm, n = .6; | ||
| double ref_bulk = 300. * p_atm, m = .2; | ||
| double ini_phi = .7, ten_fold_phi_diff = .1, r_f = .7, cohesion = .5; |
There was a problem hiding this comment.
Nine parameters shall be defined in such an order. Definitions are compatible with the ones presented in the paper.
| void new_duncanselig(unique_ptr<Material>& return_obj, istringstream& command) { | ||
| unsigned tag; | ||
| if(!get_input(command, tag)) { | ||
| suanpan_error("A valid tag is required.\n"); | ||
| return; | ||
| } | ||
|
|
||
| auto pool = vec{14.7, 400. * 14.7, .6, 300 * 14.7, .2, .7, .1, .7, .5}; | ||
| if(!get_input(command, pool)) { | ||
| suanpan_error("A valid parameter is required.\n"); | ||
| return; | ||
| } | ||
|
|
||
| auto density = 0.; | ||
| if(command.eof()) | ||
| suanpan_debug("Zero density assumed.\n"); | ||
| else if(!get_input(command, density)) { | ||
| suanpan_error("A valid density is required.\n"); | ||
| return; | ||
| } | ||
|
|
||
| return_obj = make_unique<DuncanSelig>(tag, pool, density); | ||
| } |
There was a problem hiding this comment.
This function handles the parsing of the text input file.
Material/Material2D/DuncanSelig.cpp
Outdated
| while(true) { | ||
| if(max_iteration == ++counter) { | ||
| suanpan_error("Cannot converge within {} iterations.\n", max_iteration); | ||
| return SUANPAN_FAIL; | ||
| } | ||
|
|
||
| const auto [elastic, bulk, deds, dkds] = compute_moduli(); | ||
|
|
||
| const auto factor_a = 3. * bulk * (3. * bulk + elastic); | ||
| const auto factor_b = 3. * bulk * (3. * bulk - elastic); | ||
|
|
||
| mat33 right(fill::zeros); | ||
| right(0, 0) = right(1, 1) = factor_a; | ||
| right(0, 1) = right(1, 0) = factor_b; | ||
| right(2, 2) = 3. * bulk * elastic; | ||
|
|
||
| const vec3 residual = (9. * bulk - elastic) * (trial_stress - net_stress) - right * net_strain; | ||
|
|
||
| const rowvec3 dfads = (18. * bulk + 3. * elastic) * dkds + 3. * bulk * deds; | ||
| const rowvec3 dfbds = (18. * bulk - 3. * elastic) * dkds - 3. * bulk * deds; | ||
|
|
||
| jacobian = (9. * bulk - elastic) * eye(3, 3) + (trial_stress - net_stress) * (9. * dkds - deds); | ||
| jacobian.row(0) -= net_strain(0) * dfads + net_strain(1) * dfbds; | ||
| jacobian.row(1) -= net_strain(0) * dfbds + net_strain(1) * dfads; | ||
| jacobian.row(2) -= net_strain(2) * 3. * (bulk * deds + elastic * dkds); | ||
|
|
||
| if(!solve(incre, jacobian, residual)) return SUANPAN_FAIL; | ||
|
|
||
| const auto error = inf_norm(incre); | ||
| if(1u == counter) ref_error = error; | ||
| suanpan_debug("Local iteration error: {:.5E}.\n", error / ref_error); | ||
| if(error < tolerance * ref_error || ((error < tolerance || inf_norm(residual) < tolerance) && counter > 5u)) { | ||
| if(!solve(trial_stiffness, jacobian, right)) return SUANPAN_FAIL; | ||
|
|
||
| max_dev_stress = dev(trial_stress); | ||
|
|
||
| return SUANPAN_SUCCESS; | ||
| } | ||
|
|
||
| trial_stress -= incre; | ||
| } |
There was a problem hiding this comment.
This loop effectively implements the incremental constitutive relation Eq. (3a). To obtain the consistent tangent matrix, we need to take full derivative w.r.t. stress as now the elastic and bulk moduli are functions of stress. The reference implementation provided by Prof. Katona does not follow this approach, as a result, it cannot achieve the quadratic convergence rate.
No description provided.