Two-dimensional finite-difference lattice Boltzmann method for the complete Navier-Stokes equations of binary fluids

0.18 MB·English

by  Aiguo Xu

Europhysics Letters PREPRINT
5
Two-dimensional finite-difference lattice Boltzmann method
0
0
for the complete Navier-Stokes equations of binary fluids
2

n
Aiguo Xu
a
J Department of Physics, Yoshida-South Campus, Kyoto University,

9 Sakyo-ku, Kyoto, 606-8501, Japan
1

]
h PACS.47.11.+j – Computational methods in fluid dynamics.
c
PACS.51.10.+y – Kinetic and transport theory.
e
m PACS.05.20.Db – Kinetic theory.
-
t
a
t
s Abstract. – Based on Sirovich’s two-fluidkinetic theoryand adodecagonal discrete velocity
.
t model, a two-dimensional 61-velocity finite-difference lattice Boltzmann method for the com-
a
m plete Navier-Stokes equations of binary fluids is formulated. Previous constraints, in most
existing lattice Boltzmann methods, on the studied systems, like isothermal and nearly in-
-
d compressible, are released within the present method. This method is designed to simulate
n compressible and thermal binary fluid mixtures. The validity of the proposed method is ver-
o ified by investigating (i) the Couette flow and (ii) the uniform relaxation process of the two
c components.
[

2
v
1
3 Introduction. – Lattice Boltzmann Method (LBM) has become a viable and promising
4 numericalscheme for simulating fluid flows. There are severaloptions to discretize the Boltz-
0
mann equation: (i) Standard LBM (SLBM) [1]; (ii) Finite-Difference LBM (FDLBM) [1–3];
1
(iii) Finite-Volume LBM [1,4]; (iv) Finite-Element LBM [1,5]; etc. These kinds of schemes
4
0 are expected to be complementary in the LBM studies.
/ Even though various LBMs for multicomponent fluids [6–19] have been proposed and
t
a developed,(i)mostexistingmethodsbelongtotheSLBM[6–16],and/orbasedonthesingle-
m
fluidtheory[7–14,16,17,20];(ii)inRef.[6]aSLBMbasedonSirovich’stwo-fluidkinetictheory
- [21]isproposed;(iii)nearlyallthestudiesarefocusedonisothermalandnearlyincompressible
d
systems. In a recent study [22], Sirovich’s kinetic theory is clarified and corresponding two-
n
o fluid FDLBMs for Euler equations and isothermal Navier-Stokes equations are presented. In
c thisletterweproposeatwo-fluidFDLBMforthecompleteNavier-Stokesequations,including
v: the energy equation.
i
X Formulation and verification of the FDLBM. – The formulationof a FDLBM consistsof
r three steps: (i) select or design an appropriate discrete velocity model (DVM), (ii) formulate
a
the discrete local equilibrium distribution function, (iii) choose a finite-difference scheme.
The continuous Boltzmann equation has infinite velocities, so the rotational invariance is
automaticallysatisfied. Recoveringrotationalinvariantmacroscopicequationsfromadiscrete-
finite-velocity microscopic dynamics imposes constraints on the isotropy of DVM and the
(cid:13)c EDPSciences
2
EUROPHYSICSLETTERS
finite-difference scheme used. In this Letter, the proposed FDLBM is based on the following
DVM,
iπ iπ
v =0, v =v cos ,sin ,i=1,2,···,12, (1)
0 ki k
6 6
(cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21)
where k indicates the k-th group of particle velocities and i indicates the direction of the
particle speed. It is easy find that (i) its odd rank tensors are zero, and (ii) its initial four
even rank tensors satisfy
12 v v =6v2δ , 12 v v v v = 3v4∆ ,
i=1 kiα kiβ k αβ i=1 kiα kiβ kiγ kiδ 2 k αβγδ
12 v v v v v v = 1v6∆ , (2)
Pi=1 kiα kiβ kiγ kiδ kiµPkiν 4 k αβγδµν
12 v v v v v v v v = 1 v8∆ ,
Pi=1 kiα kiβ kiγ kiδ kiµ kiν kiλ kiπ 32 k αβγδµνλπ
P
where α, β, ··· indicate x or y component and
∆ =δ δ +δ δ +δ δ , (3)
αβγδ αβ γδ αγ βδ αδ βγ
∆ =δ ∆ +δ ∆ +δ ∆ +δ ∆ +δ ∆ , (4)
αβγδµν αβ γδµν αγ βδµν αδ βγµν αµ βγδν αν βγδµ
∆ = δ ∆ +δ ∆ +δ ∆ +δ ∆
αβγδµνλπ αβ γδµνλπ αγ βδµνλπ αδ βγµνλπ αµ βγδνλπ
+δ ∆ +δ ∆ +δ ∆ . (5)
αν βγδµλπ αλ βγδµνπ απ βγδµνλ
It is clear that this DVM is isotropic up to, at least, its 9th rank tensor.
We consider a binary mixture with two components, A and B, where the masses and
temperaturesofthetwocomponentsarenotsignificantlydifferent. Theinterparticlecollisions
can be divided into two kinds: collisions within the same species (self-collision)and collisions
amongdifferent species (cross-collision)[21]. Basedonthe DVM (1), the 2-dimensionalBGK
[23] kinetic equation for species A reads [22],
∂ vA −uA
∂ fA+vA · fA−aA· ki fA(0) =JAA+JAB (6)
t ki ki ∂r ki ΘA ki ki ki
(cid:0) (cid:1)
where
JAA =− fA−fA(0) /τAA , JAB =− fA−fAB(0) /τAB (7)
ki ki ki ki ki ki
h i h i
nA vA −uA 2 nA vA −uAB 2
fA(0) = exp − ki , fAB(0) = exp − ki (8)
ki 2πΘA 2ΘA ki 2πΘAB 2ΘAB
" (cid:0) (cid:1) # " (cid:0) (cid:1) #
ΘA =k TA/mA, ΘAB =k TAB/mA (9)
B B
fA(0) and fAB(0) are the corresponding Maxwellian distribution functions. nA, uA, TA are
the local density, hydrodynamic velocity and temperature of species A. uAB, TAB are the
hydrodynamic velocity and temperature of the mixture after equilibration process. aA is the
acceleration of species A due to the effective external field.
For species A, we have
1
nA = fkAi , nAuA = vkAifkAi, PA(eAint =nAkBTA)= 2mA(vAki−uA)2fkAi (10)
ki ki ki
X X X
Aiguo Xu: 2-d FDLBM for Navier-Stokes equations 3
where PA (eA ) is the local pressure (internal mean kinetic energy). For species B, we have
int
similar relations. For the mixture, we have
1
uAB = ρAuA+ρBuB /ρ, nk TAB = (vA −uAB)2mAfA+(vB −uBA)2mBfB
B 2 ki ki ki ki
(cid:0) (cid:1) Xki h (1i1)
where ρA = nAmA, n = nA+nB and ρ = ρA +ρB. Three sets of hydrodynamic quantities
(for the two components A, B and for the mixture) are involved, but only two sets of them
are independent. So this is a two-fluid model. Without lossing generality, we focus on hy-
drodynamics of the two individual species. By expanding the local equilibrium distribution
function fAB(0) around fA(0) to the first order in flow velocity and temperature, the BGK
model (6-9) becomes
∂ vA −uA
∂ fA+vA · fA−aA· ki fA(0) =QAA+QAB (12)
t ki ki ∂r ki ΘA ki ki ki
(cid:0) (cid:1)
1 1
QAA =− + fA−fA(0) (13)
ki τAA τAB ki ki
(cid:18) (cid:19)h i
fA(0)
QAB = − ki {µA vA −uA ·(uA−uB)
ki ρAΘA D ki
vA −(cid:0)uA 2 (cid:1) vA −uA 2
+µA ki −1 (TA−TB)−MA ki −1 (uA−uB)2}(14)
T 2ΘA 2ΘA
"(cid:0) (cid:1) # "(cid:0) (cid:1) #
where µA =ρAρB/(τABρ), µA =k nAnB/(τABn), MA =nAρAρB/(2τABnρ).
D T B
Now, we go to the second step: formulate fA(0). The continuous Maxwellian fA(0) pos-
ki
sesses an infinite sequence of moment properties. The Chapman-Enskog analysis [24] shows
that, requiring the discrete fA(0) to follow the initial eight ones is sufficient to describe the
ki
same Navier-Stokes equations,
∂ρA ∂
+ ρAuA =0, (15)
∂t ∂r α
α
(cid:0) (cid:1)
∂∂t ρAuAα + ∂∂rβ ρAuAαuAβ + ∂∂PrαA −ρAaAα − ∂∂rβ ηA ∂∂urβAα + ∂∂urαAβ − ∂∂urAγγ δαβ
(cid:0) (cid:1) (cid:16) (cid:17) +ρAρB uA−uB =0(cid:20), (cid:18) (cid:19)(cid:21) (16)
τABρ α α
(cid:0) (cid:1)
∂∂etA+∂∂rα eA+PA uAα −ρAaA·uA− ∂∂rα kA∂(k∂BrTαA) +ηAuAβ ∂∂urAαβ + ∂∂urαAβ − ∂∂urAγγ δαβ
(cid:20) (cid:18) (cid:19)(cid:21)
+ρ(cid:2)A(cid:0)ρB uA 2(cid:1)−u(cid:3)A·uB + nAnBk TA−TB −nA ρAρB uA−uB 2 =0, (17)
τABρ τABn B 2τABnρ
h(cid:0) (cid:1) i (cid:0) (cid:1) (cid:0) (cid:1)
where
eA =eA +1ρA uA 2, ηA =PAτAAτAB/ τAA+τAB , kA =2nAΘAτAAτAB/ τAA+τAB .
int
2
(18)
(cid:0) (cid:1) (cid:0) (cid:1) (cid:0) (cid:1)
4
EUROPHYSICSLETTERS
Recall that uA (uB) is a small quantity. By using Eq. (15), PA = nAk TA, and neglecting
B
the secondandhigherorderterms inuA, Eq. (16)showsthat the diffusion velocity,uB−uA,
α α
is related to the gradients of nA and TA.
The first three requirements on fA(0) are referred to Eq. (10) with fA replaced by fA(0),
ki ki ki
and the remaining five are
mAvA vA fA(0) =PAδ +ρAuAuA (19)
kiα kiβ ki αβ α β
ki
X
mAvA vA vA fA(0) =PA uAδ +uAδ +uAδ +ρAuAuAuA (20)
kiα kiβ kiγ ki γ αβ α βγ β γα α β γ
ki
X (cid:0) (cid:1)
1mA vA 2vA fA(0) =2nAk TAuA+ 1ρA uA 2uA (21)
2 ki kiα ki B α 2 α
ki
X (cid:0) (cid:1) (cid:0) (cid:1)
1mA vA 2vA vA fA(0) =2PAΘAδ + 1PA uA 2δ
2 ki kiα kiβ ki αβ 2 αβ
ki
X (cid:0) (cid:1) (cid:0) (cid:1)
+3PAuAuA+ 1ρA uA 2uAuA (22)
α β 2 α β
1mA vA 4vA fA(0)(cid:0)= (cid:1)12PAΘA+6PA uA 2+ 1ρA uA 4 uA (23)
2 ki kiα ki 2 α
ki (cid:20) (cid:21)
X (cid:0) (cid:1) (cid:0) (cid:1) (cid:0) (cid:1)
The requirement equation (23) contains the fifth order of the flow velocity uA. So it is
sufficient to expand fA(0) in polynomial up to the fifth order of uA:
ki
(uA)2 (uA)4 vA uA (uA)2 (uA)4
fA(0) = nAFA 1− + + kiξ ξ 1− +
ki k (" 2ΘA 8(ΘA)2# ΘA " 2ΘA 8(ΘA)2#
vA vA uAuA (uA)2 vA vA vA uAuAuA (uA)2
kiξ kiπ ξ π kiξ kiπ kiη ξ π η
+ 1− + 1−
2(ΘA)2 2ΘA 6(ΘA)3 2ΘA
(cid:20) (cid:21) (cid:20) (cid:21)
vA vA vA vA uAuAuAuA vA vA vA vA vA uAuAuAuAuA
kiξ kiπ kiη kiλ ξ π η λ kiξ kiπ kiη kiλ kiδ ξ π η λ δ
+ +
24(ΘA)4 120(ΘA)5 )
+··· (24)
where
1 (vA)2
FA = exp − k . (25)
k 2πΘA 2ΘA
(cid:20) (cid:21)
A(0)
The truncated equilibrium distribution function f (24) contains the fifth rank tensor of
ki
the particle velocity vA and the requirement (20) contains its third rank tensor. Thus, a
DVM being isotropic up to its 8th rank tensors is enough to recover the physical isotropy of
the continuous Boltzmann equations to the Navier-Stokes level. So DVM (1) is an appropri-
ate choice. To calculate the discrete fA(0), one first needs calculate the factor FA. FA is
ki k k
determined by the eight requirements on fA(0) and the isotropic properties of the DVM (1).
ki
Following the same procedure as described in [22], we obtain
FA = 1, FA vA 2 = ΘA, FA vA 4 = 2 ΘA 2,
k k k 6 k k 3
ki k k
X X (cid:0) (cid:1) X (cid:0) (cid:1) (cid:0) (cid:1)
FA vA 6 = 4 ΘA 3, FA vA 8 =32 ΘA 4, FA vA 10 =320 ΘA 5 (26)
k k k k k k
k k k
X (cid:0) (cid:1) (cid:0) (cid:1) X (cid:0) (cid:1) (cid:0) (cid:1) X (cid:0) (cid:1) (cid:0) (cid:1)
Aiguo Xu: 2-d FDLBM for Navier-Stokes equations 5
Once a zero speed, vA =0, and other five nonzero ones, vA (k =1, 2, 3, 4, 5) are chosen, FA
0 k k
(k =0, 1, 2, 3, 4, 5) will be fixed.
Wecometothethirdstep: finite-differenceimplementationofthediscretekineticmethod.
There are several choices [17] available. One possibility is shown below,
vA −uA ∂fA,(n)
fA,(n+1) =fA,(n)+ aA· ki fA(0)+QAA,(n)+QAB,(n)−vA · ki ∆t, (27)
ki ki ΘA ki ki ki ki ∂r
" (cid:0) (cid:1) #
where the second superscripts n, n+1 indicate the consecutive two iteration steps, ∆t the
time step; the spatial derivatives are calculated as
∂fA,(n) (3fA,(n)−4fA,(n) +fA,(n) )/(2∆α) ifvA ≥0
ki = ki,I ki,I−1 ki,I−2 kiα , (28)
∂α ( (3fkAi,,I(n)−4fkAi,,I(n+)1+fkAi,,I(n+)2)/(−2∆α) ifvkAiα <0
where α = x,y, the third subscripts I −2, I −1, I, I +1, I +2 indicate consecutive mesh
nodes in the α direction.
The validity ofthe formulatedFDLBMis verifiedthroughtwotestexamples. (The Boltz-
mann constantk =1.) The firstone is the isothermalandincompressibleCouette flow with
B
a single component. In this case, A=B. The initial state of the fluid is static. The distance
betweenthetwowallsisD. Attimet=0theystarttomoveatvelocitiesU,−U,respectively.
The horizontalvelocityprofilesofspeciesAorB alongaverticallineagreewiththefollowing
analytical solution,
γD 4j2π2η 2jπ
u=γy− (−1)j+1 exp(− t)sin( y), (29)
jπ ρD2 D
j
X
where γ =2U/D is the imposedshearrate,j isaninteger,the twowallslocateaty =±D/2.
(For example, see Fig. 1.)
The secondone is the uniformrelaxationprocess,whichis anidealprocessto indicate the
equilibration behaviorof the mixture [22]. By neglecting the force terms and terms in spatial
derivatives, the Navier-Stokes equations (15)-(17) give
∂
ρA =0, (30)
∂t
∂ 1 ρA ρB
uB −uA =− + uB −uA , (31)
∂t ρ τBA τAB
(cid:18) (cid:19)
(cid:0) (cid:1) (cid:0) (cid:1)
∂ TB−TA =−1 nA + nB TB−TA + ρAρB 1 − 1 uB −uA 2. (32)
∂t n τBA τAB 2k nρ τAB τBA
(cid:0) (cid:1) (cid:18) (cid:19) B (cid:18) (cid:19)
(cid:0) (cid:1) (cid:0) (cid:1)
The flow velocities of the two components equilibrate exponentially with time. (For example,
seeFig. 2(a).) Theequilibrationofflowvelocitiesalsoaffectsthatofthetemperatures. When
the flowvelocitydifferenceiszero,thetemperaturesequilibrateexponentiallywithtime. (For
example, see Fig. 2(b).) The simulation results agree well with Eqs. (31) and (32).
Conclusions and remarks. – The Chapman-Enskog analysis shows what properties the
discrete Maxwellian distribution function fA(0) should follow. Those requirements tell the
ki
lowest order of the flow velocity uA in the Taylor expansion of fA(0). The highest rank of
ki
tensors of the particle velocity vA in the requirements on the truncated fA(0) determines
ki
6
EUROPHYSICSLETTERS
Fig. 1 Fig. 2
Fig.1– Horizontalvelocity profilesalongaverticallinefor thetwospecies, AandB,at timet=8.
Thesymbolsareforsimulation results. Thesolid linecorrespondstothetheoreticalresult,Eq. (29).
Parameters used in the two-fluid FDLBM are mA = mB = 1, T = 1, nA = nB = 1, γ = 0.001,
τAA=τBB =τAB =τBA =0.2. Parameters used in Eq. (29) are η=ηA =0.1, ρ=ρA=1.
Fig. 2 – Uniform relaxation processes. (a) Equilibration of velocities; (b) Equilibration of tem-
peratures. The symbols are for simulation results. The solid lines possess the theoretical slopes.
Common parameters for the simulations in (a) and (b) are nA = 10, nB = 1, mA = 1, mB = 10,
τAA = τBB = 1, τAB = 10, τBA = 1. In (a) the initial conditions are uA(0) = −uB(0) = −0.3,
x x
uA(0) =uB(0) =0, and TA(0) =1.3, TB(0) =0.7. The slope of the solid line in (a) is −11/20, which
y y
is consistent with Eq. (31). In (b) the initial conditions are uA(0) = uB(0) = 0, and TA(0) = 1.3,
TB(0) = 0.7. The slope of the solid line in (b) is −10.1/11, which is consistent with the first term
of right-hand side of Eq.(32). The second superscript “(0)” denotes the corresponding initial value.
This figureshows an examplewhere theparticle masses of thetwospecies are significantly different.
the needed isotropy of the DVM. The incorporation of the force terms makes no additional
requirement on the isotropy of the DVM. The present approach works for binary neutral
fluid mixtures. One possibility to introduce interfacial tension is to modify the pressure
tensors [13], which is implemented by changing the force terms [3]. The specific force terms
or pressure tensors, which are out of the scope of this Letter, depend on the system under
consideration,butcanberesolvedunderthesametwo-dimensional61-velocitymodel(D2V61).
Forbinaryfluidswithdisparate-masscomponents,saymA ≪mB,onlyifthetotalmassesand
temperatures of the two species are not significantly different, does Sirovich’s kinetic theory
works[22],so do the correspondingFDLBMs. (See Fig. 2 for an example.) When the masses
and/or the temperatures of the two components are greatly different, the two-fluid kinetic
theory should be modified. In those cases, the Navier-Stokes equations and the FDLBMs
are not symmetric about the two components, but the FDLBMs can still be resolved under
the D2V61 model. The formulation procedure is straightforward. In practical simulations,
numerical errors from the finite-difference schemes should be quantified.
∗∗∗
The author thanks Prof. G. Gonnella for guiding him into the LBM field and Profs. H.
Hayakawa, V. Sofonea, S. Succi, M. Watari for helpful discussions. This work is partially
supported by Grant-in-Aids for Scientific Research (Grant No. 15540393) and for the 21-
th Century COE “Center for Diversity and Universality in Physics” from the Minstry of
Education, Culture and Sports, Science and Technology (MEXT) of Japan.
Aiguo Xu: 2-d FDLBM for Navier-Stokes equations 7
REFERENCES
[1] Succi S., The Lattice Boltzmann Equation ((Oxford University Press, New York) 2001; Benzi
R., Succi S. and Vergassola M., Phys. Rep., 222 (1992) 3; Higuera F., Succi. S. and
Benzi R.,Europhys. Lett., 9 (1989) 345.
[2] Cao N., Chen S.,Jin S. and Martinez D., Phys. Rev. E, 55(1997)R21; Seta T. and Taka-
hashi R., J. Stat. Phys., 107(2002)557; Sofonea V. and Sekerka R. F. , J. Comp. Phys.,
184(2003)422; Watari M. and Tsutahara M., Phys. Rev. E, 67(2003)36306; Kataoka T.
and Tsutahara M., Phys. Rev. E, 69(2004)56702; Kataoka T. and Tsutahara M., Phys.
Rev. E, 69(2004)R35701.
[3] Sofonea V., Lamura A., Gonnella G., and Cristea A., Phys. Rev. E, 70(2004)046702.
[4] Nannelli F., and Succi S., J. Stat. Phys., 68(401)1992; Xi H., Peng G. and Chou S. H.,
Phys. Rev. E, 60 (1999)3380; Ubertini S., Bella G., and Succi S., Phys. Rev. E, 68 (2003)
16701; See also references in [1].
[5] Li Y., LeBoeufe E. and Basu P.K., Phys. Rev. E, 69 (2004) 065701(R).
[6] Luo L. S. and Girimaji S. S., Phys. Rev. E, 66(2002)R35301; ibid.67(2003)36302.
[7] Guo Z. and Zhao T. S., Phys. Rev. E, 68(2003)R35302.
[8] Xu Aiguo, Gonnella G. and Lamura A., Phys. Rev. E, 67(2003)56105; Physica A,
331(2004)10; Physcia A, (in press); cond-mat/0404205,; Xu Aiguo, Commun. Theor. Phys.,
39(2003)729.
[9] Gunstensen A. K. , Rothman D. H., Zaleski S., and Zanetti G., Phys. Rev. A,
43(1991)4320.
[10] Flekkoy E. G., Phys. Rev. E, 47(1993)4247.
[11] Grunau D., Chen S., and Eggert K., Phys. Fluids A, 5(1993)2557.
[12] Shan X. and Doolen G., J. Stat. Phys., 81(1995)379; Phys. Rev. E, 54(1996)3614.
[13] Orlandini E., Osborn W. R., and Yeomans J. M., Europhys. Lett., 32(1995)463; Osborn
W. R., Orlandini E., Swift M. R., Yeomans J. M., and Banavar J. R., Phys. Rev.
Lett., 75(1995)4031; Swift M. R., Orlandini E., Osborn W. R., and Yeomans J. M.,
Phys. Rev. E,54(1996)5041; Gonnella G., Orlandini E., Yeomans J. M.,Phys. Rev. Lett,
78(1997)1695, Phys. Rev. E, 58(1998) 480; Lamura A., Gonnella G., and Yeomans J. M.,
Europhys. Lett., 45(1999)314.
[14] KendonV.M.,DesplatJ.C.,BladonP.,andCatesM.E.,Phys.Rev.Lett.,83(1999)576;
Kendon V. M., Cates M. E., Pagonabarrage I., DesplatJ. C.,andBladon P.,J. Fluid
Mech., 440(2001)147.
[15] Yu H., Luo L. S., Girimaji S. S., Int. J. Comp. Eng. Sci., 3(2002)73.
[16] Facin P. C., Philippi P. C.,anddos Santos L. O. E.,inICCS2003, LNCS 2657,editedby
Sloot P. M. A. et al. (Spring-Verlag, Heidelberg) 2003, p. 1007-1014.
[17] Cristea A. and Sofonea V., Central European J. Phys., 2(2004)382.
[18] Shan X. and Chen H.,Phys. Rev. E, 47(1993)1815; ibid,49(1994)2941.
[19] Love P. J. and Coveney P.V., Phys. Rev. E, 64(2001)21503; Love P. J.,Maillet J. B.,
and Coveney P. V.,ibid,64(2001)61302; Chin J. andCoveney P. V.,ibid,66(2002) 16303;
Gonzalez-Segredo N, Nekovee M., and Coveney P.V., ibid,67(2003)46304; Gonzalez-
Segredo N. and Coveney P. V., ibid,69(2004)61501.
[20] Sofonea V. and Sekerka R. F., Physica A, 299(2001)494.
[21] Sirovich L., Phys. Fluids, 5(1962)908; ibid.9(1966)2323.
[22] Xu Aiguo, http://arxiv.org/abs/cond-mat/0406012.
[23] Bhatnagar P. L., Gross E. P., and Krook M., Phys. Rev. , 94(1954)511.
[24] ChapmannS.andCowling T.G.,TheMathematical TheoryofNon-UniformGases,3rded.,
(Cambridge UniversityPress, Cambridge) 1970.