. . . . (- simpleFoam, -- pUCoupledFoam) 10 0 Turbulent energy - - PowerPoint PPT Presentation

simplefoam pucoupledfoam 10 0 turbulent energy 10 1
SMART_READER_LITE
LIVE PREVIEW

. . . . (- simpleFoam, -- pUCoupledFoam) 10 0 Turbulent energy - - PowerPoint PPT Presentation

. Block coupled calculations in OpenFOAM Chalmers University of Technology Klas Jareteg . Implementing a pressure-velocity coupled solver October 22, 2012 . . . . . . (- simpleFoam, -- pUCoupledFoam) 10 0 Turbulent energy 10 -1


slide-1
SLIDE 1

. .

.

.

.

.

.

.

Block coupled calculations in OpenFOAM

Implementing a pressure-velocity coupled solver

200 400 600 800 1000 Iteration 10

  • 10

10

  • 9

10

  • 8

10

  • 7

10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 Residuals (- simpleFoam, -- pUCoupledFoam) Turbulent energy Velocity (x) Pressure

Klas Jareteg

Chalmers University of Technology

October 22, 2012

slide-2
SLIDE 2

.

.

Overview

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

. .Overview . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 2

slide-3
SLIDE 3

.

.

.

.

Presentation overview

  • Introduction
  • Block coupling formulation
  • OpenFOAM background
  • Pressure-velocity coupled solver
  • Tests and results
  • Conclusion

.Overview . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 3

slide-4
SLIDE 4

.

.

.

.

Learning outcomes

  • General understanding of the block matrix structure
  • Overview of the BlockLduMatrix class with utilities
  • Understand implementation structure of block coupled solvers
  • Increased general knowledge on discretization and assembling

.Overview . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 4

slide-5
SLIDE 5

.

.

Introduction

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

. .Introduction . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 5

slide-6
SLIDE 6

.

.

.

.

Introduction

  • Many problems coupled PDE. Example Navier-Stokes:

∇ · (U) = 0 (1) ∇ · (UU) − ∇(ν∇U) = −1 ρ∇p (2)

  • Traditionally segragated solvers
  • Alternatively: coupled solvers
  • OpenFOAM 1.6-ext: BlockLduMatrix

. .Introduction . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 6

slide-7
SLIDE 7

.

.

.

.

Introduction

  • Many problems coupled PDE. Example Navier-Stokes:

∇ · (U) = 0 (1) ∇ · (UU) − ∇(ν∇U) = −1 ρ∇p (2)

  • Traditionally segragated solvers
  • Alternatively: coupled solvers
  • OpenFOAM 1.6-ext: BlockLduMatrix

. .Introduction . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 6

slide-8
SLIDE 8

.

.

.

.

Introduction

  • Many problems coupled PDE. Example Navier-Stokes:

∇ · (U) = 0 (1) ∇ · (UU) − ∇(ν∇U) = −1 ρ∇p (2)

  • Traditionally segragated solvers
  • Alternatively: coupled solvers
  • OpenFOAM 1.6-ext: BlockLduMatrix

. .Introduction . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 6

slide-9
SLIDE 9

.

.

.

.

Introduction

  • Many problems coupled PDE. Example Navier-Stokes:

∇ · (U) = 0 (1) ∇ · (UU) − ∇(ν∇U) = −1 ρ∇p (2)

  • Traditionally segragated solvers
  • Alternatively: coupled solvers
  • OpenFOAM 1.6-ext: BlockLduMatrix

. .Introduction . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 6

slide-10
SLIDE 10

.

.

.

.

Introduction

  • Many problems coupled PDE. Example Navier-Stokes:

∇ · (U) = 0 (1) ∇ · (UU) − ∇(ν∇U) = −1 ρ∇p (2)

  • Traditionally segragated solvers
  • Alternatively: coupled solvers
  • OpenFOAM 1.6-ext: BlockLduMatrix

. .Introduction . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 6

slide-11
SLIDE 11

.

.

Formulation

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

. .Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 7

slide-12
SLIDE 12

.

.

.

.

Formulation

  • General coupled matrix system:

A(y)x = a (3) B(x)y = b (4)

  • Solved together (still segregated):

[A(y) B(x) ] [x y ] = [a b ] (5)

  • Coupled solution:

[A′ Ay Bx B′ ] [x y ] = [a b ] (6)

. .Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 8

slide-13
SLIDE 13

.

.

.

.

Formulation

  • General coupled matrix system:

A(y)x = a (3) B(x)y = b (4)

  • Solved together (still segregated):

[A(y) B(x) ] [x y ] = [a b ] (5)

  • Coupled solution:

[A′ Ay Bx B′ ] [x y ] = [a b ] (6)

. .Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 8

slide-14
SLIDE 14

.

.

.

.

Formulation

  • General coupled matrix system:

A(y)x = a (3) B(x)y = b (4)

  • Solved together (still segregated):

[A(y) B(x) ] [x y ] = [a b ] (5)

  • Coupled solution:

[A′ Ay Bx B′ ] [x y ] = [a b ] (6)

. .Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 8

slide-15
SLIDE 15

.

.

.

.

Block coupled approach

  • Alternative formulation:

Cz = c (7) C = Ci,j = [ ca,a

i,j

ca,b

i,j

cb,a

i,j

cb,b

i,j

]

i,j

(8) c = ci = [ ai bi ]⊤ (9) z = zi = [ xi yi ]⊤ (10)

  • Element in vectors and matrices: vectors and tensors

.Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 9

slide-16
SLIDE 16

.

.

.

.

Block coupled approach

  • Alternative formulation:

Cz = c (7) C = Ci,j = [ ca,a

i,j

ca,b

i,j

cb,a

i,j

cb,b

i,j

]

i,j

(8) c = ci = [ ai bi ]⊤ (9) z = zi = [ xi yi ]⊤ (10)

  • Element in vectors and matrices: vectors and tensors

.Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 9

slide-17
SLIDE 17

.

.

.

.

Block coupled approach

  • Alternative formulation:

Cz = c (7) C = Ci,j = [ ca,a

i,j

ca,b

i,j

cb,a

i,j

cb,b

i,j

]

i,j

(8) c = ci = [ ai bi ]⊤ (9) z = zi = [ xi yi ]⊤ (10)

  • Element in vectors and matrices: vectors and tensors

.Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 9

slide-18
SLIDE 18

.

.

.

.

Block coupled approach

  • Coupled system is larger
  • Implicit coupling → faster convergence
  • Convergence properties changed

Non-linearities in equations still left!

. .Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 10

slide-19
SLIDE 19

.

.

.

.

Block coupled approach

  • Coupled system is larger
  • Implicit coupling → faster convergence
  • Convergence properties changed

Non-linearities in equations still left!

. .Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 10

slide-20
SLIDE 20

.

.

.

.

Block coupled approach

  • Coupled system is larger
  • Implicit coupling → faster convergence
  • Convergence properties changed

Non-linearities in equations still left!

.Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 10

slide-21
SLIDE 21

.

.

.

.

Block coupled approach

  • Coupled system is larger
  • Implicit coupling → faster convergence
  • Convergence properties changed

Non-linearities in equations still left!

. .Formulation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 10

slide-22
SLIDE 22

.

.

OpenFOAM background

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

. OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 11

slide-23
SLIDE 23

.

.

.

.

Discretization, assembling and matrices I

polyMesh and fvMesh Face based computational cell:

OWNER

NEIGHBOUR

Sf

face f cell i cell j

fvMesh is the finite volume specialization V() Volumes of the cells. Numbered according to cell numbering. Sf() Surface normals with magnitude equal to the

  • area. Numbered according to face numbers.

. . OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 12

slide-24
SLIDE 24

.

.

.

.

Discretization, assembling and matrices II

lduMatrix Basic square sparse matrix. Stored in three arrays: the diagonal, the upper and the lower part:

85 // − C o e f f i c i e n t s ( not including i n t e r f a c e s ) 86 scalarField *lowerPtr_ , *diagPtr_ , *upperPtr_ ;

Listing 1: lduMatrix.H

  • Diagonal elements: numbered as cell numbers
  • Off-diagonal elements: are numbered according to faces.

1 const surfaceVectorField& Sf = p . mesh ( ) . Sf ( ) ; 2 const unallocLabelList& owner = mesh . owner ( ) ; 3 const unallocLabelList& neighbour = mesh . neighbour ( ) ;

Listing 2: Surface normal, owner and neighbour for each face

. OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 13

slide-25
SLIDE 25

.

.

.

.

Discretization, assembling and matrices III

Sparsity of matrix: A = Ai,j (11)

  • i, j: contribution from cell j on cell i
  • j, i: contribution from cell i on cell j
  • i > j: upper elements
  • i < j: lower elements
  • i = j: diagonal elements

. . OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 14

slide-26
SLIDE 26

.

.

.

.

Discretization, assembling and matrices IV

fvMatrix

  • Specialization for finite volume
  • Adds source and reference to field
  • Helper functions:

3 volScalarField AU = UEqn ( ) . A ( ) ; 4 U = UEqn ( ) . H()/ AU ;

Listing 3: Part of pEqn.H in simpleFoam

. . OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 15

slide-27
SLIDE 27

.

.

.

.

Discretization, assembling and matrices V

fvm and fvc

  • In general implicit and explicit operators:

12 fvScalarMatrix pEqn 13 ( 14 fvm : : laplacian (1.0/ AU , p) == fvc : : div(phi) 15 ) ;

Listing 4: Part of pEqn.H in simpleFoam

  • Explicit→source, Implicit→source,matrix
  • Mathematical operators defined in matrix class

. OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 16

slide-28
SLIDE 28

.

.

.

.

Discretization, assembling and matrices VI

Boundary conditions

  • Working on faces in a patch
  • Calculates contribution to diagonal and source.
  • E.g. fixedValue:

139 // − Return the matrix diagonal c o e f f i c i e n t s corresponding to the 140 // evaluation

  • f

the value

  • f

t h i s patchField with given weights 141 v i r t u a l tmp<Field<Type> > valueInternalCoeffs 142 ( 143 const tmp<scalarField>& 144 ) const ; 145 146 // − Return the matrix source c o e f f i c i e n t s corresponding to the 147 // evaluation

  • f

the value

  • f

t h i s patchField with given weights 148 v i r t u a l tmp<Field<Type> > valueBoundaryCoeffs 149 ( 150 const tmp<scalarField>& 151 ) const ; 152 153 // − Return the matrix diagonal c o e f f i c i e n t s corresponding to the 154 // evaluation

  • f

the gradient

  • f

t h i s patchField 155 v i r t u a l tmp<Field<Type> > gradientInternalCoeffs () const ; 156 157 // − Return the matrix source c o e f f i c i e n t s corresponding to the 158 // evaluation

  • f

the gradient

  • f

t h i s patchField 159 v i r t u a l tmp<Field<Type> > gradientBoundaryCoeffs () const ;

Listing 5: Part of fixedValueFvPatchField.H

. . OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 17

slide-29
SLIDE 29

.

.

.

.

Discretization, assembling and matrices VII

  • Usage examples:

92 forAll (fvm . psi ( ) . boundaryField () , patchI ) 93 { 94 const fvPatchField<Type>& psf = fvm . psi ( ) . boundaryField ( ) [ patchI ] ; 95 const fvsPatchScalarField& patchFlux = faceFlux . boundaryField ( ) [ patchI ] ; 96 const fvsPatchScalarField& pw = weights . boundaryField ( ) [ patchI ] ; 97 98 fvm . internalCoeffs ( ) [ patchI ] = patchFlux*psf . valueInternalCoeffs (pw ) ; 99 fvm . boundaryCoeffs ( ) [ patchI ] = −patchFlux*psf . valueBoundaryCoeffs (pw ) ; 100 }

Listing 6: Part of gaussConvectionScheme.C

70 forAll (fvm . psi ( ) . boundaryField () , patchI ) 71 { 72 const fvPatchField<Type>& psf = fvm . psi ( ) . boundaryField ( ) [ patchI ] ; 73 const fvsPatchScalarField& patchGamma = 74 gammaMagSf . boundaryField ( ) [ patchI ] ; 75 76 fvm . internalCoeffs ( ) [ patchI ] = patchGamma*psf . gradientInternalCoeffs ( ) ; 77 fvm . boundaryCoeffs ( ) [ patchI ] = −patchGamma*psf . gradientBoundaryCoeffs ( ) ; 78 }

Listing 7: Part of gaussianLaplacianScheme.C

. . OpenFOAM background . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 18

slide-30
SLIDE 30

.

.

Pressure-velocity coupling model

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

. Pressure-velocity coupling model . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 19

slide-31
SLIDE 31

.

.

.

.

Pressure-velocity coupling I

  • Navier-Stokes, incompressible, steady-state:

∇ · (U) = 0 (12) ∇ · (UU) − ∇(ν∇U) = −1 ρ∇p (13)

  • Semi-discretized form:

faces

Uf · Sf = 0 (14) ∑

faces

[UU − µ∇U]f · Sf = − ∑

faces

PfSf (15)

  • Modified pressure:

p ρ = P (16)

  • Rhie-Chow in continuity equation:

faces

[ Uf − Df ( ∇Pf − ∇Pf )] · Sf = 0 (17)

. . Pressure-velocity coupling model . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 20

slide-32
SLIDE 32

.

.

OpenFOAM implementation

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

. OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 21

slide-33
SLIDE 33

.

.

.

.

Implementation I

Solution vector: xP = xP

l =

    uP vP wP PP    

l

(18)

117 // Block vector f i e l d for the pressure and v e l o c i t y f i e l d to be solved for 118 volVector4Field pU 119 ( 120 IOobject 121 ( 122 ”pU” , 123 runTime . timeName () , 124 mesh , 125 IOobject : : NO_READ , 126 IOobject : : NO_WRITE 127 ) , 128 mesh , 129 dimensionedVector4 (word () , dimless , vector4 : : zero) 130 ) ; . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 22

slide-34
SLIDE 34

.

.

.

.

Implementation II

132 // I n s e r t the pressure and v e l o c i t y i n t e r n a l f i e l d s in to the volVector2Field 133 { 134 vector4Field blockX = pU . internalField ( ) ; 135 136 // Separately add the three v e l o c i t y components 137 for ( i n t i=0; i<3;i++) 138 { 139 tmp<scalarField> tf = U . internalField ( ) . component (i ) ; 140 scalarField& f = tf ( ) ; 141 blockMatrixTools : : blockInsert (i , f , blockX ) ; 142 } 143 144 // Pressure i s the 2nd component 145 scalarField& f = p . internalField ( ) ; 146 blockMatrixTools : : blockInsert (3 ,f , blockX ) ; 147 } . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 23

slide-35
SLIDE 35

.

.

.

.

Implementation III

Equation system to be formed: APxP + ∑

F∈{N}

AFxF = bP (19) AX = [ aX

k,l

]

i

k, l ∈ {u, v, w, p}, X ∈ {P, F} (20) Construct block matrix:

188 // Matrix block 189 BlockLduMatrix<vector4> B(mesh ) ;

Retrieve fields:

191 // Diagonal i s set separately 192 Field<tensor4>& d = B . diag ( ) . asSquare ( ) ; 193 194 // Off−diagonal also as square 195 Field<tensor4>& u = B . upper ( ) . asSquare ( ) ; 196 Field<tensor4>& l = B . lower ( ) . asSquare ( ) ;

Source:

198 // Source term for the block matrix 199 Field<vector4> s(mesh . nCells () , vector4 : : zero ) ; . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 24

slide-36
SLIDE 36

.

.

.

.

Discretizing the momentum equation I

LHS: Turbulence is introduced by calling the divDivReff(U)

182 tmp<fvVectorMatrix> UEqnLHS 183 ( 184 fvm : : div(phi , U) 185 + turbulence− >divDevReff (U) 186 ) ;

Retrieve matrix coefficients:

202 tmp<scalarField> tdiag = UEqnLHS ( ) . D ( ) ; 203 scalarField& diag = tdiag ( ) ; 204 scalarField& upper = UEqnLHS ( ) . upper ( ) ; 205 scalarField& lower = UEqnLHS ( ) . lower ( ) ;

Add boundary contribution:

211 // Add source boundary contribution 212 vectorField& source = UEqnLHS ( ) . source ( ) ; 213 UEqnLHS ( ) . addBoundarySource (source , f a l s e ) ; . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 25

slide-37
SLIDE 37

.

.

.

.

Discretizing the momentum equation II

Considering RHS as separate problem: ∑

faces

PfSf = 0 (21) Interpolation weights:

218 // I n t e r p o l a t i o n scheme for the pressure weights 219 tmp<surfaceInterpolationScheme<scalar> > 220 tinterpScheme_ 221 ( 222 surfaceInterpolationScheme<scalar >::New 223 ( 224 p . mesh () , 225 p . mesh ( ) . interpolationScheme (”grad (p)”) 226 ) 227 ) ; 218 tmp<surfaceScalarField> tweights = tinterpScheme_ ( ) . weights (p ) ; 219 const surfaceScalarField& weights = tweights ( ) ;

wN = 1 − wP (22)

. . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 26

slide-38
SLIDE 38

.

.

.

.

Discretizing the momentum equation III

Equivalent to matrix fields:

229 // Pressure gradient contributions − corresponds to an i m p l i c i t 230 // gradient

  • perator

231 tmp<vectorField> tpUv = tmp<vectorField> 232 ( 233 new vectorField ( upper . size () , pTraits<vector >:: zero) 234 ) ; 235 vectorField& pUv = tpUv ( ) ; 236 tmp<vectorField> tpLv = tmp<vectorField> 237 ( 238 new vectorField ( lower . size () , pTraits<vector >:: zero) 239 ) ; 240 vectorField& pLv = tpLv ( ) ; 241 tmp<vectorField> tpSv = tmp<vectorField> 242 ( 243 new vectorField ( source . size () , pTraits<vector >:: zero) 244 ) ; 245 vectorField& pSv = tpSv ( ) ; 246 tmp<vectorField> tpDv = tmp<vectorField> 247 ( 248 new vectorField (diag . size () , pTraits<vector >:: zero) 249 ) ; 250 vectorField& pDv = tpDv ( ) ; . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 27

slide-39
SLIDE 39

.

.

.

.

Discretizing the momentum equation IV

Calcualte elements:

256 for ( i n t i=0;i<owner . size ( ) ; i++) 257 { 258 i n t o = owner [ i ] ; 259 i n t n = neighbour [ i ] ; 260 scalar w = weights . internalField ( ) [ i ] ; 261 vector s = Sf [ i ] ; 262 263 pDv [ o]+=s*w ; 264 pDv [ n] −=s*(1−w ) ; 265 pLv [ i]= −s*w ; 266 pUv [ i]=s*(1−w ) ; 267 268 } . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 28

slide-40
SLIDE 40

.

.

.

.

Discretizing the momentum equation V

Boundary contribution:

271 p . boundaryField ( ) . updateCoeffs ( ) ; 272 forAll (p . boundaryField () , patchI ) 273 { 274 // Present fvPatchField 275 fvPatchField<scalar> & fv = p . boundaryField ( ) [ patchI ] ; 276 277 // Retrieve the weights for the boundary 278 const fvsPatchScalarField& pw = weights . boundaryField ( ) [ patchI ] ; 279 280 // Contributions from the boundary c o e f f i c i e n t s 281 tmp<Field<scalar> > tic = fv . valueInternalCoeffs (pw ) ; 282 Field<scalar>& ic = tic ( ) ; 283 tmp<Field<scalar> > tbc = fv . valueBoundaryCoeffs (pw ) ; 284 Field<scalar>& bc = tbc ( ) ; 285 286 // Get the fvPatch

  • nly

287 const fvPatch& patch = fv . patch ( ) ; 288 289 // Surface normals for t h i s patch 290 tmp<Field<vector> > tsn = patch . Sf ( ) ; 291 Field<vector> sn = tsn ( ) ; 292 293 // Manually add the contributions from the boundary 294 // This what happens with addBoundaryDiag , addBoundarySource 295 forAll (fv , facei ) 296 { 297 label c = patch . faceCells ( ) [ facei ] ; 298 299 pDv [ c]+=ic [ facei ]* sn [ facei ] ; 300 pSv [ c] −=bc [ facei ]* sn [ facei ] ; 301 } 302 } . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 29

slide-41
SLIDE 41

.

.

.

.

Discretizing the momentum equation VI

au,u, av,v, aw,w, ap,u, ap,v, ap,w:

317 forAll (d , i) 318 { 319 d [ i ](0 ,0) = diag [ i ] ; 320 d [ i ](1 ,1) = diag [ i ] ; 321 d [ i ](2 ,2) = diag [ i ] ; 322 323 d [ i ](0 ,3) = pDv [ i ] . x ( ) ; 324 d [ i ](1 ,3) = pDv [ i ] . y ( ) ; 325 d [ i ](2 ,3) = pDv [ i ] . z ( ) ; 326 } 327 forAll (l , i) 328 { 329 l [ i ](0 ,0) = lower [ i ] ; 330 l [ i ](1 ,1) = lower [ i ] ; 331 l [ i ](2 ,2) = lower [ i ] ; 332 333 l [ i ](0 ,3) = pLv [ i ] . x ( ) ; 334 l [ i ](1 ,3) = pLv [ i ] . y ( ) ; 335 l [ i ](2 ,3) = pLv [ i ] . z ( ) ; 336 } 337 forAll (u , i) 338 { 339 u [ i ](0 ,0) = upper [ i ] ; 340 u [ i ](1 ,1) = upper [ i ] ; 341 u [ i ](2 ,2) = upper [ i ] ; 342 343 u [ i ](0 ,3) = pUv [ i ] . x ( ) ; 344 u [ i ](1 ,3) = pUv [ i ] . y ( ) ; 345 u [ i ](2 ,3) = pUv [ i ] . z ( ) ; 346 } 347 forAll (s , i) 348 { 349 s [ i ](0) = source [ i ] . x()+pSv [ i ] . x ( ) ; . . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 30

slide-42
SLIDE 42

.

.

.

.

Discretizing the momentum equation VII

350 s [ i ](1) = source [ i ] . y()+pSv [ i ] . y ( ) ; 351 s [ i ](2) = source [ i ] . z()+pSv [ i ] . z ( ) ; 352 } . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 31

slide-43
SLIDE 43

.

.

.

.

Discretizing the continuity equation I

One implicit and one explicitcontribution:

439 tmp<volScalarField> tA = UEqnLHS ( ) . A ( ) ; 440 volScalarField& A = tA ( ) ; 442 tmp<volVectorField> texp = fvc : : grad (p ) ; 443 volVectorField& exp = texp ( ) ; 444 tmp<volVectorField> texp2 = exp/A ; 445 volVectorField exp2 = texp2 ( ) ; 446 447 tmp<fvScalarMatrix> MEqnLHSp 448 ( 449 −fvm : : laplacian(1/A , p) 450 = = 451 −fvc : : div(exp2 ) 452 ) ; 454 // Add the boundary contributions 455 scalarField& pMdiag = MEqnLHSp ( ) . diag ( ) ; 456 scalarField& pMupper = MEqnLHSp ( ) . upper ( ) ; 457 scalarField& pMlower = MEqnLHSp ( ) . lower ( ) ; 458 459 // Add diagonal boundary contribution 460 MEqnLHSp ( ) . addBoundaryDiag (pMdiag , 0 ) ; 461 462 // Add source boundary contribution 463 scalarField& pMsource = MEqnLHSp ( ) . source ( ) ; 464 MEqnLHSp ( ) . addBoundarySource (pMsource , f a l s e ) ;

Need implicit divergence scheme:

. OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 32

slide-44
SLIDE 44

.

.

.

.

Discretizing the continuity equation II

348 // Again an i m p l i c i t version not existing , now the div

  • perator

349 tmp<surfaceInterpolationScheme<scalar> > 350 UtinterpScheme_ 351 ( 352 surfaceInterpolationScheme<scalar >::New 353 ( 354 U . mesh () , 355 U . mesh ( ) . interpolationScheme (” div (U)( i m p l i c i t )”) 356 ) 357 ) ; 358 359 360 // 1) Setup diagonal , source , upper and lower 361 tmp<vectorField> tMUpper = tmp<vectorField> 362 (new vectorField ( upper . size () , pTraits<vector >:: zero ) ) ; 363 vectorField& MUpper = tMUpper ( ) ; 364 365 tmp<vectorField> tMLower = tmp<vectorField> 366 (new vectorField ( lower . size () , pTraits<vector >:: zero ) ) ; 367 vectorField& MLower = tMLower ( ) ; 368 369 tmp<vectorField> tMDiag = tmp<vectorField> 370 (new vectorField (diag . size () , pTraits<vector >:: zero ) ) ; 371 vectorField& MDiag = tMDiag ( ) ; 372 373 tmp<vectorField> tMSource = tmp<vectorField> 374 ( 375 new vectorField 376 ( 377 source . component ( 0 ) ( ) . size () , pTraits<vector >:: zero 378 ) 379 ) ; 380 vectorField& MSource = tMSource ( ) ; 381 . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 33

slide-45
SLIDE 45

.

.

.

.

Discretizing the continuity equation III

382 // 2) Use i n t e r p o l a t i o n weights to assemble the contributions 383 tmp<surfaceScalarField> tMweights = 384 UtinterpScheme_ ( ) . weights (mag(U ) ) ; 385 const surfaceScalarField& Mweights = tMweights ( ) ; 386 387 for ( i n t i=0;i<owner . size ( ) ; i++) 388 { 389 i n t o = owner [ i ] ; 390 i n t n = neighbour [ i ] ; 391 scalar w = Mweights . internalField ( ) [ i ] ; 392 vector s = Sf [ i ] ; 393 394 MDiag [ o]+=s*w ; 395 MDiag [ n] −=s*(1−w ) ; 396 MLower [ i]= −s*w ; 397 MUpper [ i]=s*(1−w ) ; 398 } 399 400 // Get boundary condition contributions for the pressure grad (P) 401 U . boundaryField ( ) . updateCoeffs ( ) ; 402 forAll (U . boundaryField () , patchI ) 403 { 404 // Present fvPatchField 405 fvPatchField<vector> & fv = U . boundaryField ( ) [ patchI ] ; 406 407 // Retrieve the weights for the boundary 408 const fvsPatchScalarField& Mw = 409 Mweights . boundaryField ( ) [ patchI ] ; 410 411 // Contributions from the boundary c o e f f i c i e n t s 412 tmp<Field<vector> > tic = fv . valueInternalCoeffs (Mw ) ; 413 Field<vector>& ic = tic ( ) ; 414 tmp<Field<vector> > tbc = fv . valueBoundaryCoeffs (Mw ) ; 415 Field<vector>& bc = tbc ( ) ; 416 . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 34

slide-46
SLIDE 46

.

.

.

.

Discretizing the continuity equation IV

417 // Get the fvPatch

  • nly

418 const fvPatch& patch = fv . patch ( ) ; 419 420 // Surface normals for t h i s patch 421 tmp<Field<vector> > tsn = patch . Sf ( ) ; 422 Field<vector> sn = tsn ( ) ; 423 424 // Manually add the contributions from the boundary 425 // This what happens with addBoundaryDiag , addBoundarySource 426 forAll (fv , facei ) 427 { 428 label c = patch . faceCells ( ) [ facei ] ; 429 430 MDiag [ c]+=cmptMultiply (ic [ facei ] , sn [ facei ] ) ; 431 MSource [ c] −=cmptMultiply (bc [ facei ] , sn [ facei ] ) ; 432 } 433 } . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 35

slide-47
SLIDE 47

.

.

.

.

Discretizing the continuity equation V

au,p, av,p, aw,p, ap,p:

469 forAll (d , i) 470 { 471 d [ i ](3 ,0) = MDiag [ i ] . x ( ) ; 472 d [ i ](3 ,1) = MDiag [ i ] . y ( ) ; 473 d [ i ](3 ,2) = MDiag [ i ] . z ( ) ; 474 d [ i ](3 ,3) = pMdiag [ i ] ; 475 } 476 forAll (l , i) 477 { 478 l [ i ](3 ,0) = MLower [ i ] . x ( ) ; 479 l [ i ](3 ,1) = MLower [ i ] . y ( ) ; 480 l [ i ](3 ,2) = MLower [ i ] . z ( ) ; 481 l [ i ](3 ,3) = pMlower [ i ] ; 482 } 483 forAll (u , i) 484 { 485 u [ i ](3 ,0) = MUpper [ i ] . x ( ) ; 486 u [ i ](3 ,1) = MUpper [ i ] . y ( ) ; 487 u [ i ](3 ,2) = MUpper [ i ] . z ( ) ; 488 u [ i ](3 ,3) = pMupper [ i ] ; 489 } 490 forAll (s , i) 491 { 492 s [ i ](3) = MSource [ i ] . x () 493 +MSource [ i ] . y () 494 +MSource [ i ] . z () 495 +pMsource [ i ] ; 496 } . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 36

slide-48
SLIDE 48

.

.

.

.

Solving the block matrix and retrieving the solution I

Solve the system:

501 BlockSolverPerformance<vector4> solverPerf = 502 BlockLduSolver<vector4 >::New 503 ( 504 word(”blockVar” ) , 505 B , 506 mesh . solver (”blockVar”) 507 ) − >solve (pU , s ) ; 508 509 solverPerf . print ( ) ; . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 37

slide-49
SLIDE 49

.

.

.

.

Solving the block matrix and retrieving the solution II

Retrieve the solution:

515 tmp<scalarField> tUx = U . internalField ( ) . component (0); 516 scalarField& Ux = tUx ( ) ; 517 blockMatrixTools : : blockRetrieve (0 , Ux , pU ) ; 518 U . internalField ( ) . replace (0 , Ux ) ; 519 520 tmp<scalarField> tUy = U . internalField ( ) . component (1); 521 scalarField& Uy = tUy ( ) ; 522 blockMatrixTools : : blockRetrieve (1 , Uy , pU ) ; 523 U . internalField ( ) . replace (1 , Uy ) ; 524 525 tmp<scalarField> tUz = U . internalField ( ) . component (2); 526 scalarField& Uz = tUz ( ) ; 527 blockMatrixTools : : blockRetrieve (2 , Uz , pU ) ; 528 U . internalField ( ) . replace (2 , Uz ) ; 529 530 blockMatrixTools : : blockRetrieve (3 , p . internalField () , pU ) ; 531 532 UEqnLHS . clear ( ) ; 533 534 p . relax ( ) ; 535 536 U . correctBoundaryConditions ( ) ; 537 p . correctBoundaryConditions ( ) ; . OpenFOAM implementation . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 38

slide-50
SLIDE 50

.

.

Results

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

.Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 39

slide-51
SLIDE 51

.

.

.

.

Results - 2D I

Figure: Pressure profile for solution of the dailyPitz using pUCoupledFoam.

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 40

slide-52
SLIDE 52

.

.

.

.

Results - 2D II

Figure: Pressure comparison for simpleFoam and pUCoupledFoam at line indicated in figure 1

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 41

slide-53
SLIDE 53

.

.

.

.

Results - 2D III

200 400 600 800 1000 Iteration 10

  • 10

10

  • 9

10

  • 8

10

  • 7

10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 Residuals

(- simpleFoam, -- pUCoupledFoam) Turbulent energy Velocity (x) Pressure

(a) Iterations

10 20 30 40 50 60 70 80 90 ['Elapsed time [s]'] 10

  • 10

10

  • 9

10

  • 8

10

  • 7

10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 Residuals

(- simpleFoam, -- pUCoupledFoam) Turbulent energy Velocity (x) Pressure

(b) Elapsed time

Figure: Comparison of convergence for simpleFoam and pUCoupledFoam. Turbulent case.

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 42

slide-54
SLIDE 54

.

.

.

.

Results - 2D IV

200 400 600 800 1000 Iteration 10

  • 10

10

  • 9

10

  • 8

10

  • 7

10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 Residuals

(- simpleFoam, -- pUCoupledFoam) Velocity (x) Pressure

(a) Iterations

5 10 15 20 25 30 35 40 ['Elapsed time [s]'] 10

  • 10

10

  • 9

10

  • 8

10

  • 7

10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 Residuals

(- simpleFoam, -- pUCoupledFoam) Velocity (x) Pressure

(b) Elapsed time

Figure: Comparison of convergence for simpleFoam and pUCoupledFoam. Laminar case.

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 43

slide-55
SLIDE 55

.

.

.

.

Results - 3D I

Figure: Pressure profile comparison for (y, z) = (0.075, 0.075).

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 44

slide-56
SLIDE 56

.

.

.

.

Results - 3D II

Figure: Velocity profile comparison for (y, z) = (0.075, 0.075).

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 45

slide-57
SLIDE 57

.

.

.

.

Results - 3D III

10 20 30 40 50 60 70 80 Iteration 10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 Residuals

(- pUCoupledFoam3D, -- simpleFoam) Turbulent energy Velocity (x) Pressure

(a) Iterations

200 400 600 800 1000 1200 1400 1600 1800 ['Elapsed time [s]'] 10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 Residuals

(- pUCoupledFoam3D, -- simpleFoam) Turbulent energy Velocity (x) Pressure

(b) Elapsed time

Figure: Comparison of convergence of simpleFoam and pUCoupledFoam3D (50 iterations).

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 46

slide-58
SLIDE 58

.

.

.

.

Results - 3D IV

Figure: Parts of the mesh generated by snappyHexMesh in the motorBike case.

. .Results . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 47

slide-59
SLIDE 59

.

.

Conclusions

Overview Introduction Formulation OpenFOAM background Pressure-velocity coupling model OpenFOAM implementation Results Conclusions

.Conclusions . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 48

slide-60
SLIDE 60

.

.

.

.

Conclusions

  • Working in 2D and 3D
  • Significant performance gain (primarily 2D and laminar)
  • Turbulence is limiting the performance gain from coupling pressure and

velocity

  • Much left to do (parallelizing, coupled patches, operator implementation

...)

. .Conclusions . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 49

slide-61
SLIDE 61

.

.

Questions?

.Conclusions . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 50

slide-62
SLIDE 62

.

.

.Conclusions . Klas Jareteg, Project in CFD with OpenSource software, 2012-10-22 . 51