ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixing bug in link pressure

Browse Source
This commit is contained in:
jtclemm
2023-09-19 13:37:47 +02:00
parent ec65fc48ad
commit 77a5fd16dd
1 changed files with 7 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/fix_deform.cpp
View File

@ -1302,6 +1302,7 @@ void FixDeform::set_volume()
double p2 = pressure->vector[fixed]; double p2 = pressure->vector[fixed];
double p1i = set[i].prior_pressure; double p1i = set[i].prior_pressure;
double p2i = set[fixed].prior_pressure; double p2i = set[fixed].prior_pressure;
double denominator;
if (e3 == 0) { if (e3 == 0) {
e1 = 0.0; e1 = 0.0;
@ -1314,7 +1315,12 @@ void FixDeform::set_volume()
if (!linked_pressure) { if (!linked_pressure) {
// Calculate first strain rate by expanding stress to linear order, p1(t+dt) = p2(t+dt) // Calculate first strain rate by expanding stress to linear order, p1(t+dt) = p2(t+dt)
// Calculate second strain rate to preserve volume // Calculate second strain rate to preserve volume
e1 = -e3 / (1 + e3 * dt) * (p2 - p2i) - e2i * (p1 - p2) / (p2 - p2i + (p1 - p1i) / e1i * e2i); denominator = (p2 - p2i + (p1 - p1i) / e1i * e2i);
if (denominator != 0.0 && e1i != 0.0) {
e1 = (-e3 / (1 + e3 * dt) * (p2 - p2i) - e2i * (p1 - p2)) / denominator;
} else {
e1 = e2i;
}
e2 = (1.0 - (1 + e3 * dt) * (1 + e1 * dt)) / ((1 + e3 * dt) * (1 + e1 * dt) * dt); e2 = (1.0 - (1 + e3 * dt) * (1 + e1 * dt)) / ((1 + e3 * dt) * (1 + e1 * dt) * dt);
shift = 0.5 * L1i * (1.0 + e1 * dt); shift = 0.5 * L1i * (1.0 + e1 * dt);