Mercurial > repos > public > sbplib

comparison +scheme/ViscoElastic2d.m @ 1314:58df4a35fe43 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add traction operators that include the viscous terms and use them to clean up the imposition of traction BC.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 25 Jul 2020 15:56:31 -0700
parents a4c8bf2371f3
children 6650b111742d
comparison
equal deleted inserted replaced
1313:a4c8bf2371f3 1314:58df4a35fe43
39 % Bundary restriction ops 39 % Bundary restriction ops
40 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n 40 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n
41 41
42 n_w, n_e, n_s, n_n % Physical normals 42 n_w, n_e, n_s, n_n % Physical normals
43 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents 43 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
44
45 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
46 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
47 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
48 tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
44 49
45 elasticObj 50 elasticObj
46 51
47 end 52 end
48 53
221 226
222 obj.e_scalar_w = elasticObj.e_scalar_w; 227 obj.e_scalar_w = elasticObj.e_scalar_w;
223 obj.e_scalar_e = elasticObj.e_scalar_e; 228 obj.e_scalar_e = elasticObj.e_scalar_e;
224 obj.e_scalar_s = elasticObj.e_scalar_s; 229 obj.e_scalar_s = elasticObj.e_scalar_s;
225 obj.e_scalar_n = elasticObj.e_scalar_n; 230 obj.e_scalar_n = elasticObj.e_scalar_n;
231
232 % -- Create new traction operators including viscous strain contribution --
233 tau1 = struct;
234 tau2 = struct;
235 tau_n = struct;
236 tau_t = struct;
237 boundaries = {'w', 'e', 's', 'n'};
238 for bNumber = 1:numel(boundaries)
239 b = boundaries{bNumber};
240
241 n = elasticObj.getNormal(b);
242 t = elasticObj.getTangent(b);
243 e = elasticObj.getBoundaryOperatorForScalarField('e', b);
244 tau1.(b) = (elasticObj.getBoundaryOperator('tau1', b)'*Eu')';
245 tau2.(b) = (elasticObj.getBoundaryOperator('tau2', b)'*Eu')';
246
247 % Add viscous contributions
248 for i = 1:dim
249 for k = 1:dim
250 for l = 1:dim
251 tau1.(b) = tau1.(b) - (n{i}*e'*C{i,1,k,l}*eGamma{k,l}')';
252 tau2.(b) = tau2.(b) - (n{i}*e'*C{i,2,k,l}*eGamma{k,l}')';
253 end
254 end
255 end
256
257 % Compute normal and tangential components
258 tau_n.(b) = tau1.(b)*n{1} + tau2.(b)*n{2};
259 tau_t.(b) = tau1.(b)*t{1} + tau2.(b)*t{2};
260 end
261
262 obj.tau1_w = tau1.w;
263 obj.tau1_e = tau1.e;
264 obj.tau1_s = tau1.s;
265 obj.tau1_n = tau1.n;
266
267 obj.tau2_w = tau2.w;
268 obj.tau2_e = tau2.e;
269 obj.tau2_s = tau2.s;
270 obj.tau2_n = tau2.n;
271
272 obj.tau_n_w = tau_n.w;
273 obj.tau_n_e = tau_n.e;
274 obj.tau_n_s = tau_n.s;
275 obj.tau_n_n = tau_n.n;
276
277 obj.tau_t_w = tau_t.w;
278 obj.tau_t_e = tau_t.e;
279 obj.tau_t_s = tau_t.s;
280 obj.tau_t_n = tau_t.n;
281 %----------------------------------------
226 282
227 % Misc. 283 % Misc.
228 obj.elasticObj = elasticObj; 284 obj.elasticObj = elasticObj;
229 obj.m = elasticObj.m; 285 obj.m = elasticObj.m;
230 obj.h = elasticObj.h; 286 obj.h = elasticObj.h;
278 penalty = Eu*penalty; 334 penalty = Eu*penalty;
279 335
280 switch component 336 switch component
281 case 't' 337 case 't'
282 dir = obj.getTangent(boundary); 338 dir = obj.getTangent(boundary);
339 tau = obj.getBoundaryOperator('tau_t', boundary);
283 case 'n' 340 case 'n'
284 dir = obj.getNormal(boundary); 341 dir = obj.getNormal(boundary);
342 tau = obj.getBoundaryOperator('tau_n', boundary);
285 case 1 343 case 1
286 dir = {1, 0}; 344 dir = {1, 0};
345 tau = obj.getBoundaryOperator('tau1', boundary);
287 case 2 346 case 2
288 dir = {0, 1}; 347 dir = {0, 1};
348 tau = obj.getBoundaryOperator('tau2', boundary);
289 end 349 end
290 350
291 switch type 351 switch type
292 case {'F','f','Free','free','traction','Traction','t','T'} 352 case {'F','f','Free','free','traction','Traction','t','T'}
293 353
294 % Add viscous part of closure 354 % Set elastic closure to zero
295 for i = 1:dim 355 closure = 0*closure;
296 for j = 1:dim 356
297 for k = 1:dim 357 for m = 1:dim
298 for l = 1:dim 358 closure = closure - eU{m}*( (RHO*H)\(e*dir{m}*H_gamma*tau') );
299 for m = 1:dim
300 closure = closure + eU{m}*( (RHO*H)\(C{i,j,k,l}*e*dir{m}*dir{j}*H_gamma*n{i}*e'*eGamma{k,l}') );
301 end
302 end
303 end
304 end
305 end 359 end
306 360
307 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} 361 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
308 362
309 % Add penalty to strain rate eq 363 % Add penalty to strain rate eq
512 % op -- string 566 % op -- string
513 function o = getBoundaryOperator(obj, op, boundary) 567 function o = getBoundaryOperator(obj, op, boundary)
514 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 568 assertIsMember(boundary, {'w', 'e', 's', 'n'})
515 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'}) 569 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
516 570
517 o = obj.elasticObj.([op, '_', boundary]); 571 o = obj.([op, '_', boundary]);
518 572
519 end 573 end
520 574
521 % Returns the boundary operator op for the boundary specified by the string boundary. 575 % Returns the boundary operator op for the boundary specified by the string boundary.
522 % op -- string 576 % op -- string