Skip to content

Commit 0f63379

Browse files
mrp089mrp089
authored
Add loss models, test STENOSIS (#76)
1 parent 525e7cd commit 0f63379

File tree

4 files changed

+292
-3
lines changed

4 files changed

+292
-3
lines changed

‎Code/Source/cvOneDModelManager.cxx‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -97,6 +97,18 @@ int cvOneDModelManager::CreateSegment(char *segName,long segID, double segLen
9797
// convert char string to boundary condition type
9898
if (!strcmp(lossType, "NONE")){
9999
loss = MinorLossScope::NONE;
100+
}else if(!strcmp(lossType, "STENOSIS")){
101+
loss = MinorLossScope::STENOSIS;
102+
}else if(!strcmp(lossType, "BRANCH_THROUGH_DIVIDING")){
103+
loss = MinorLossScope::BRANCH_THROUGH_DIVIDING;
104+
}else if(!strcmp(lossType, "BRANCH_SIDE_DIVIDING")){
105+
loss = MinorLossScope::BRANCH_SIDE_DIVIDING;
106+
}else if(!strcmp(lossType, "BRANCH_THROUGH_CONVERGING")){
107+
loss = MinorLossScope::BRANCH_THROUGH_CONVERGING;
108+
}else if(!strcmp(lossType, "BRANCH_SIDE_CONVERGING")){
109+
loss = MinorLossScope::BRANCH_SIDE_CONVERGING;
110+
}else if(!strcmp(lossType, "BIFURCATION_BRANCH")){
111+
loss = MinorLossScope::BIFURCATION_BRANCH;
100112
}else{
101113
return CV_ERROR;
102114
}

‎Code/Source/cvOneDMthSegmentModel.cxx‎

Lines changed: 210 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -110,15 +110,222 @@ double cvOneDMthSegmentModel::N_MinorLoss(long ith){
110110
S[1] = currSolution->Get( eqNumbers[0]);
111111
Q[1] = currSolution->Get( eqNumbers[1]);
112112

113+
if(minorLoss == MinorLossScope::NONE ){//|| minorLoss != MinorLossScope::STENOSIS ){
114+
return sub->GetMaterial()->GetN(S[1]);
115+
}
113116

114-
if(minorLoss == MinorLossScope::NONE ){
117+
// In case of branch, previous seg might not be adjacent seg
118+
// Subdomain *adjacent =subdomainList[ith-1];
119+
cvOneDSubdomain *upstream = subdomainList[sub->GetUpstreamSeg()];
120+
cvOneDSubdomain *branch = subdomainList[sub->GetBranchSeg()];
115121

116-
return sub->GetMaterial()->GetN(S[1]);
122+
double L = sub->GetLength();
123+
double Kv, Kt, Re0, La, D[2], func, N, theta, min;
124+
strcpy(propName,"kinematic viscosity");
125+
double kinViscosity0 = upstream->GetMaterial()->GetProperty(propName);
126+
127+
min=0.6; // if a for branches is too small, troubles abound
128+
129+
// The last node of segment (3) at the upstream of the stenosis (minor loss) segment // not true for converging flow
130+
GetEquationNumbers(upstream->GetNumberOfElements()-1, eqNumbers, sub->GetUpstreamSeg());
131+
S[0] = currSolution->Get( eqNumbers[2]);
132+
Q[0] = currSolution->Get( eqNumbers[3]);
133+
134+
// first node of the branch (minor loss) segment, only used for through(2) calcs
135+
GetEquationNumbers(0, eqNumbers, sub->GetBranchSeg());
136+
S[2] = currSolution->Get( eqNumbers[0]);
137+
Q[2] = currSolution->Get( eqNumbers[1]);
138+
139+
// ratio of flow rate between branch and combined flow
140+
double q = Q[1]/Q[0];
141+
// ratio of area between branch and combined flow
142+
double a = S[1]/S[0];
117143

144+
// why are we doing this? Sometimes flow is negative.
145+
if(Q[0] < 0.000001 || Q[1] < 0.000001){
146+
return sub->GetMaterial()->GetN(S[1]);
118147
}
119148

149+
switch(minorLoss){
150+
case MinorLossScope::STENOSIS:
151+
Kt = 1.52;
152+
// modify to lessen stenosis
153+
D[0] = sqrt(4*S[0]/3.14159265358979323846);
154+
D[1] = sqrt(4*S[1]/3.14159265358979323846);
155+
La = 0.83*L + 1.64*D[1];
156+
// La = L;
157+
Kv = 32. * La / D[0] * 1.0/a*1.0/a;
158+
Re0 = D[0] * Q[0] / (S[0] * kinViscosity0);
159+
func= Kv/Re0 + Kt / 2. * (1.0/a - 1.) * (1.0/a - 1.);
160+
// func= Kv/Re0 + Kt / 2. * (1/a - 1.) * (1/a - 1.)*abs(Q[0])/Q[0]+Ku*1.06*L*sub->getDQDt()/Q; // for non-steady flow
161+
func = 2.0 * func*a*a/q/q; // a^2/q^2 switch from upstream segment to stenosed segment
162+
break;
163+
164+
case MinorLossScope::BRANCH_THROUGH_DIVIDING:
165+
theta = sub -> GetBranchAngle();
166+
167+
q = Q[2]/Q[0]; // branch(1)/combined(0)
168+
a = S[2]/S[0];
169+
170+
// having problems with a too small, try hokey fix for now 02-08-02
171+
if(a<min){
172+
// a=min;
173+
// return sub->GetMaterial()->GetN(S[1]);
174+
}
175+
// Wood K32, modified, times combined flow div by leg flow
176+
func = (-0.03*(1-q)*(1-q) - 0.35 * q*q + 0.2*q*(1-q)); // no difference
177+
// modified coeff -- too strong
178+
// func=func * Q[0]*Q[0]/(Q[1]*Q[1]);
179+
// cout<<"through dividing q "<<q<< " func "<<func<<endl;
180+
break;
181+
182+
// problem with this case.
183+
case MinorLossScope::BRANCH_SIDE_DIVIDING:
184+
theta = sub -> GetBranchAngle();
185+
186+
// having problems with a too small, try hokey fix for now 02-08-02
187+
if(a<min){
188+
// a=min;
189+
// return sub->GetMaterial()->GetN(S[1]);
190+
}
191+
// q= Q[1]/Q[0]; // branch(this=1)/combined(0)
192+
193+
// Wood K31
194+
func = (-0.95*(1-q)*(1-q) + q*q*(1.3*cot((3.14159265358979323846-theta)/2) - 0.3 + (.4 - .1*a)/(a*a)) +
195+
-0.4*q*(1-q)*(1+(1/a))*cot((3.14159265358979323846 - theta)/2));// one sign change
196+
// func = 0.95*(1-q)*(1-q) + q*q*(1.3*cot((3.14159265358979323846-theta)/2) - 0.3 + (.4 - .1*a)/(a*a)) +
197+
// 0.4*q*(1-q)*(1+(1/a))*cot((3.14159265358979323846 - theta)/2);
198+
199+
// modified coefficient -- too strong
200+
// func=func/(q*q);
201+
// cout<<"side dividing q "<<q <<" func "<<func<< endl;
202+
break;
203+
204+
case MinorLossScope::BRANCH_THROUGH_CONVERGING:
205+
206+
// The segment at the upstream of the stenosis (minor loss) segment // not true for converging flow
207+
GetEquationNumbers(0, eqNumbers, sub->GetUpstreamSeg());
208+
S[0] = currSolution->Get( eqNumbers[0]);
209+
Q[0] = currSolution->Get( eqNumbers[1]);
210+
211+
// the stenosis (minor loss) segment // not true for converging flow
212+
GetEquationNumbers(subdomainList[ith]->GetNumberOfElements()-1, eqNumbers, ith);
213+
S[1] = currSolution->Get( eqNumbers[2]);
214+
Q[1] = currSolution->Get( eqNumbers[3]);
215+
216+
theta = sub -> GetBranchAngle();
217+
218+
q = Q[2]/Q[0]; // branch/combined
219+
a = S[2]/S[0];
220+
// having problems with a too small, try hokey fix for now 02-08-02
221+
if(a<min){
222+
a = min;
223+
return sub->GetMaterial()->GetN(S[1]);
224+
}
225+
226+
// this method from Marc Serre Using upstream area instead of downstreram area in divisor a
227+
// it also wants branch flow instead of through flow, so modify Q for 90 degree
228+
// func = 2*(1-0.2)*q-q*q+.04*q*q*(1/sqrt(a*a*a));
229+
230+
// Serre for arbitrary angle
231+
//func = (1.61*q+(-1+ .04*(1/sqrt(a*a*a)) - 1.74*(1/a)*cos(theta))*q*q);
232+
//cout<<"Serre method conv-through "<<func<<endl;
233+
234+
// Wood K23
235+
func = (-0.03*(1+q)*(1+q) + q*q*(1 + 1.62*(cos(theta)/a - 1) - 0.38*(1-a))+(2-a)*q*(1+q)); // sign change in front of last term and in last term
236+
// func =(0.03*(1-q)*(1-q) - q*q*(1 + 1.62*(cos(theta)/a - 1) - 0.38*(1-a))+(2-a)*q*(1-q));
237+
// modified value, divide by flow ratio squared. I found that this was too strong.
238+
// func= func * Q[0]*Q[0]/(Q[1]*Q[1]);
239+
// cout <<" branch converging q:"<<q<<"func:"<< func ;
240+
break;
241+
242+
case MinorLossScope::BRANCH_SIDE_CONVERGING:
243+
// The segment of combined flow (downstream), want first segment
244+
GetEquationNumbers(0, eqNumbers, sub->GetUpstreamSeg());
245+
S[0] = currSolution->Get(eqNumbers[0]);
246+
Q[0] = currSolution->Get(eqNumbers[1]);
247+
248+
// the branch segment. want last element
249+
GetEquationNumbers(subdomainList[ith]->GetNumberOfElements()-1, eqNumbers, ith);
250+
S[1] = currSolution->Get(eqNumbers[2]);
251+
Q[1] = currSolution->Get(eqNumbers[3]);
252+
253+
q = Q[1]/Q[0]; // branch(1)/combined(3)
254+
a = S[1]/S[0];
255+
// having problems with a too small, try hokey fix for now 02-08-02
256+
if(a < min){
257+
a = min;
258+
return sub->GetMaterial()->GetN(S[1]);
259+
}
260+
theta = sub -> GetBranchAngle();
261+
// Marc Serre eq29 q=branch/upstream
262+
// func=(1-1.8*a)*((q*q)/(a*a) - 1);
263+
// cout<< "Serre method - conv-branch "<< func <<endl;
264+
// Wood K13 // try changing first and last terms to (1+q)
265+
func =(0.92*(1+q)*(1+q)+(q*q)*(1.2*(cos(theta)/a-1)+0.8*(1-1/(a*a))-(1-a)*cos(theta)/a) +(2-a)*q*(1+q) );
266+
// func =-0.92*(1-q)*(1-q)-(q*q)*(1.2*(cos(theta)/a-1)+0.8*(1-1/(a*a))-(1-a)*cos(theta)/a) +(2-a)*q*(1-q) ;
267+
// use modified value..
268+
//func=func/(q*q);
269+
// cout<<"side converging " <<func<< " q "<<q ;
270+
break;
271+
case MinorLossScope::BIFURCATION_BRANCH:
272+
273+
theta = sub ->GetBranchAngle();
274+
275+
/* JNIEnv* env=BFSolver::jenv;
276+
jobject obj=BFSolver::jobj;// OneDModel object
277+
// cout<< "in bif-branch " << endl;
278+
jclass cls = env->GetObjectClass(obj);
279+
// cout<< "jclass " << cls << endl;
280+
281+
// jmethodID mid = env->GetMethodID(cls, "callback", "(DDD)D");
282+
jmethodID mid = env->GetMethodID(cls, "findBifurcationK", "(DDD)D");
283+
if (mid == 0) {
284+
cout<< "method id is zero " << endl;
285+
return sub->GetMaterial()->GetProperty( "N");
286+
}
287+
jdouble ans = env->CallDoubleMethod(obj, mid, theta,Q[0],S[0]);
288+
// cout<< ">>found method id ans = "<<ans << endl;
289+
func = ans;
290+
*/
291+
//this is turned off because we don't use java in the debug version
292+
strcpy(propName,"N");
293+
return sub->GetMaterial()->GetProperty(propName);
294+
break;
295+
296+
/*
297+
case MinorLossScope::BEND_90:
298+
func = .45;
299+
break;
300+
case MinorLossScope::BEND_45:
301+
func = .2;
302+
break;
303+
case MinorLossScope::BEND_180: // don't think I can really model this, case
304+
func = 0.4;
305+
break;
306+
default:
307+
return sub->GetMaterial()->GetProperty( "N");
308+
*/
309+
}// end of switch
310+
311+
// integrate (1) over z to get formula for dP
312+
// jing's
313+
// N = func * Q[0] * Q[0] * S[1] * S[1] / (S[0] * S[0] * Q[1] * L);
314+
// mine.. simple
315+
316+
N = func * Q[1] / (2 * L);
317+
318+
sub->SaveK(N,(int)(cvOneDBFSolver::currentTime / cvOneDBFSolver::deltaTime) );
319+
320+
// don't want to have less than the default Puoseille
321+
strcpy(propName,"N");
322+
double std= sub->GetMaterial()->GetProperty(propName);
323+
if( -N > std){
324+
return std;
325+
}
120326

121-
return 0;
327+
// cout << " -N " << -N << " Q(1) :"<< Q[1] << endl;
328+
return -N;
122329

123330
}
124331

‎test/run_tests.py‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,12 @@ def get_tests():
6161
Test('flow', 0, -1, -1, 100.0, 1.0e-16, 'point'),
6262
Test('area', 0, -1, -1, 1.0, 1.0e-5, 'point')]
6363

64+
tests['tube_stenosis_r'] = [Test('pressure', 0, 0, -1, 10150.68211, 1.0e-6, 'point'),
65+
Test('pressure', 2, -1, -1, 10000.0, 1.0e-8, 'point'),
66+
Test('flow', 0, -1, -1, 100.0, 1.0e-16, 'point'),
67+
Test('area', 0, -1, -1, 10.0, 1.0e-4, 'point')]
68+
69+
6470
tests['bifurcation_P'] = [Test('pressure', 0, 0, -1, 4039.45953118937, 1e-5, 'point'),
6571
Test('pressure', 0, -1, -1, 4026.67220709878, 1e-5, 'point'),
6672
Test('pressure', 1, 0, -1, 4026.67220709878, 1e-5, 'point'),

‎test/tube_stenosis_r.in‎

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
MODEL results_tube_stenosis_r_
2+
3+
NODE 0 0.0 0.0 0.0
4+
NODE 1 0.0 0.0 10.0
5+
NODE 2 0.0 0.0 20.0
6+
NODE 3 0.0 0.0 30.0
7+
8+
JOINT JOINT0 1 IN0 OUT0
9+
JOINTINLET IN0 1 0
10+
JOINTOUTLET OUT0 1 1
11+
12+
JOINT JOINT1 2 IN1 OUT1
13+
JOINTINLET IN1 1 1
14+
JOINTOUTLET OUT1 1 2
15+
16+
SEGMENT seg0 0 10.0 50 0 1 10.0 10.0 0.0 MAT1 NONE 0.0 0 0 NOBOUND NONE
17+
SEGMENT seg1 1 10.0 50 1 2 5.0 5.0 0.0 MAT1 STENOSIS 0.0 0 0 NOBOUND NONE
18+
SEGMENT seg2 2 10.0 50 2 3 10.0 10.0 0.0 MAT1 NONE 0.0 0 0 RESISTANCE OUTLETDATA
19+
20+
DATATABLE INLETDATA LIST
21+
0.0 100.0
22+
10.0 100.0
23+
ENDDATATABLE
24+
25+
DATATABLE OUTLETDATA LIST
26+
0.0 100.0
27+
ENDDATATABLE
28+
29+
SOLVEROPTIONS 0.001 1000 1000 2 INLETDATA FLOW 1.0e-8 1 1
30+
31+
MATERIAL MAT1 OLUFSEN 1.06 0.04 0 2.0 1.0e15 -20 1e9
32+
33+
OUTPUT TEXT
34+
35+
# analytical solution
36+
37+
# parameters
38+
# viscosity mu 0.04
39+
# vessel length L 10
40+
# vessel cross-section normal S0 10
41+
# vessel cross-section stenosis S1 5
42+
# resistance BC R_BC 100
43+
# prescribed inflow Q 100
44+
# loss factor Kt 1.52
45+
46+
# derived parameters
47+
# vessel diameter normal D0 3.568248232
48+
# vessel diameter stenosis D1 2.523132522
49+
# kinematic visocity nu 0.03773584906
50+
51+
# stenosis loss model
52+
# corrected length La = 0.83*L + 1.64*D1 = 12.43793734
53+
# loss factor Kv = 32*La/D0*(S0/S1)^2 = 446.1729889
54+
# Reynolds number Re0 = D0*Q/S0/nu
55+
# loss factor K = 2*(Kv/Re0 + Kt/2 * (S0/S1 - 1)^2) * (S1/S0)^2
56+
57+
# reference solution
58+
# resistance normal vessels R1 = R3 = 8*mu*L*PI/A^2 = 10.05309649
59+
# pressure drop stenosis delta P = K/2*rho * (Q/S1)^2 = 130.5759137
60+
61+
62+
# results to be checked
63+
# pressure inlet Pin = Q*(R1+R3+R_BC) + delta P = 10150.68211
64+
# pressure outlet Pout = Q*R_BC = 10000

0 commit comments

Comments
 (0)