Skip to content

Hgebouw improved #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
arnoldkoopman wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Hgebouw improved #9

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
092c7b6
Hgebouw
Jun 24, 2022
b7fbd16
Hgebouw
Jun 26, 2022
85e9dbb
Hgebouw
Jun 26, 2022
6f21484
Merge branch 'HgebouwImproved' of https://github.com/rivm-syso/OURS i…
Jul 16, 2022
a53c98a
Hgebouw
Jul 17, 2022
b160ad0
Hgebouw debugje
Sep 19, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
117 changes: 89 additions & 28 deletions Naverwerking/Hgebouw.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,12 @@
Created on Thu Sep 12 12:16:59 2019

@author: Arnold Koopman

version tracker
1.01 24 juni 2022 Arnold Koopman
- Hfundering improved, for large buildings
- Impedance jump limited, based on consideration of low building stiffness
- Tau effect reduced, based on limitation of foundation stiffness
"""

import argparse
Expand All @@ -19,7 +25,7 @@
ondergrenzen = octaafbanden/math.sqrt(2);
omega = 2*np.pi*octaafbanden;

MCgrootte = 10; #3*333; # streven: 333
MCgrootte = 3*333; # streven: 333
np.random.seed(1234); # om te zorgen voor steeds dezelfde output, laatste cijfer kan nl. wat zwabberen

# rekenwaardes
Expand All @@ -33,17 +39,19 @@
var_zetaBeton = .1 ;

# defaultwaardes, voor als gebruiker niets invoert
StandaardWandlengte = 10; # meter, diepte van de woning
StandaardGevellengte = 6; # meter
StandaardAantalBouwlagen = 2; # 2, dus BG en eerste verdieping (zonder dak)
StandaardGebouwHoogte = 5.6;# meter
StandaardVloerHoogte = 2.8; # eerste verdieping
StandaardVar_wandlengte = 1; # in meters, 2*std
StandaardVar_gevellengte = 1; # in meters, 2*std
StandaardVar_aantalBouwlagen = 1; # in aantalverdiepingen, 2*std
StandaardVar_frequenties = .1; # voor zowel Mid- als Quarterspan
StandaardWandlengte = 10; # meter, diepte van de woning
StandaardGevellengte = 6; # meter
StandaardAantalBouwlagen = 2; # 2, dus BG en eerste verdieping (zonder dak)
StandaardGebouwHoogte = 5.6;# meter
StandaardGebouwC1Hz = 180; # m/s, buiggolfsnelheid, over de hoogte van het gebouw, bij 1Hz
StandaardVloerHoogte = 2.8; # eerste verdieping
StandaardVar_wandlengte = 1; # in meters, 2*std
StandaardVar_gevellengte = 1; # in meters, 2*std
StandaardVar_aantalBouwlagen = 1; # in aantalverdiepingen, 2*std
StandaardVar_frequenties = .1; # voor zowel Mid- als Quarterspan
StandaardVar_gebouwHoogte = 2.8;# in meters, 2*std
StandaardVar_vloerHoogte = 2.8;# in meters, 2*std
StandaardVar_gebouwC1Hz = 50; # in meters, 2*std
StandaardVar_vloerHoogte = 2.8; # in meters, 2*std
StandaardVloeroverspanningHout = 4;
StandaardVloeroverspanningBeton = 6;
# onzekerheden over vloerfrequenties, indien die freqs nog berekend moeten worden
Expand Down Expand Up @@ -117,6 +125,28 @@ def Hgebouw(Bodem,Gebouw,Vloer):
if var_gebouwHoogte==0:
var_gebouwHoogte=.01;

if "gebouwBuigfrequentie" in Gebouw: # in Hz
gebouwBuigfrequentie = np.array(Gebouw["gebouwBuigfrequentie"]);
else:
gebouwBuigfrequentie = [];
if len(gebouwBuigfrequentie)==0:
gebouwC1Hz = StandaardGebouwC1Hz;
var_gebouwC1Hz = StandaardVar_gebouwC1Hz;
else:
gebouwC1Hz = 4*gebouwHoogte*np.sqrt(gebouwBuigfrequentie);
if "var_gebouwBuigfrequentie" in Gebouw: # in Hz, 2*std
var_gebouwBuigfrequentie = np.array(Gebouw["var_gebouwBuigfrequentie"]);
else:
var_gebouwBuigfrequentie = [];
if len(var_gebouwBuigfrequentie)==0: # maw: leeg
var_gebouwC1Hz = StandaardVar_gebouwC1Hz;
else:
var_gebouwC1Hz = np.sqrt((var_gebouwBuigfrequentie * gebouwC1Hz/(2*gebouwBuigfrequentie))**2 +
(4*var_gebouwHoogte*np.sqrt(gebouwBuigfrequentie))**2);
if var_gebouwC1Hz==0:
var_gebouwC1Hz = StandaardVar_gebouwC1Hz;


# vloerhoogte bepalen
# vloerhoogte kan expliciet worden gegeven of via VerdiepingNr
# evenzo voor de onzekerheid ervan
Expand Down Expand Up @@ -228,7 +258,8 @@ def Hgebouw(Bodem,Gebouw,Vloer):
HfxxmcArray = np.zeros([MCgrootte,aantalbanden]);
HfzzmcArray = np.zeros([MCgrootte,aantalbanden]);
ZbMC = np.zeros([aantalbanden,MCgrootte], dtype=complex);
ZgMC = np.zeros([aantalbanden,MCgrootte], dtype=complex);
ZgXMC = np.zeros([aantalbanden,MCgrootte], dtype=complex);
ZgZMC = np.zeros([aantalbanden,MCgrootte], dtype=complex);
cZMC = np.zeros([aantalbanden,MCgrootte]);
cXMC = np.zeros([aantalbanden,MCgrootte]);

Expand All @@ -239,6 +270,7 @@ def Hgebouw(Bodem,Gebouw,Vloer):
gebouwhoogteMC = stats.lognorm.ppf(kansenreeks,var_gebouwHoogte/gebouwHoogte/2) * gebouwHoogte;
vloerhoogteMC = stats.norm.ppf(kansenreeks,vloerHoogte,var_vloerHoogte/2);
vloerhoogteMC = np.where(vloerhoogteMC<0,0,vloerhoogteMC);
gebouwC1HzMC = stats.lognorm.ppf(kansenreeks,var_gebouwC1Hz/gebouwC1Hz/2) * gebouwC1Hz;
wandlengteMC = stats.lognorm.ppf(kansenreeks,var_wandlengte/wandlengte/2) * wandlengte;
gevellengteMC = stats.lognorm.ppf(kansenreeks,var_gevellengte/gevellengte/2) * gevellengte;
zetaMC = stats.lognorm.ppf(kansenreeks,var_zeta) * zeta;
Expand All @@ -248,6 +280,7 @@ def Hgebouw(Bodem,Gebouw,Vloer):
np.random.shuffle(wandlengteMC);
np.random.shuffle(gevellengteMC);
np.random.shuffle(zetaMC);
np.random.shuffle(gebouwC1HzMC);
# fase, Y en c zijn spectra, dus covariantie tussen banden meenemen
# dat kan met np.random.multivariate, maar helaas alleen voor normale verdelingen
# En dus helaas geen pseudorandom.
Expand All @@ -272,11 +305,21 @@ def Hgebouw(Bodem,Gebouw,Vloer):
#cZMC[band] = stats.lognorm.ppf(kansenreeks,var_c[band]) * c[band];
#np.random.shuffle(cZMC[band]);
cXMC[band] = c_ratioMC[band] * cZMC[band];
ZgMC[band] = 1j * gebouwdichtheidMC * gebouwhoogteMC * wandlengteMC * gevellengteMC * omega[band]; # complexe impedantie van het gebouw
Zb = np.swapaxes(ZbMC,0,1);
Zg = np.swapaxes(ZgMC,0,1);
cZ = np.swapaxes(cZMC,0,1);
cX = np.swapaxes(cXMC,0,1);
lamdbaX = cXMC[band] / octaafbanden[band];
lamdbaZ = cZMC[band] / octaafbanden[band];
# ZgMC[band] = 1j * gebouwdichtheidMC * gebouwhoogteMC * wandlengteMC * gevellengteMC * omega[band]; # complexe impedantie van het gebouw
VolumeTotaal = gebouwhoogteMC * wandlengteMC * gevellengteMC;
VolumeActiefX = np.pi * lamdbaX**3 / 100;
VolumeActiefX = np.minimum(VolumeTotaal,VolumeActiefX);
VolumeActiefZ = np.pi * lamdbaZ**3 / 100;
VolumeActiefZ = np.minimum(VolumeTotaal,VolumeActiefZ);
ZgXMC[band] = 1j * gebouwdichtheidMC * VolumeActiefX * omega[band]; # complexe impedantie van het gebouw
ZgZMC[band] = 1j * gebouwdichtheidMC * VolumeActiefZ * omega[band]; # complexe impedantie van het gebouw
Zb = np.swapaxes(ZbMC ,0,1);
ZgX = np.swapaxes(ZgXMC,0,1);
ZgZ = np.swapaxes(ZgZMC,0,1);
cZ = np.swapaxes(cZMC ,0,1);
cX = np.swapaxes(cXMC ,0,1);

for i1 in range(MCgrootte): # Monte Carlo over variaties, om covariantiematrix te kunnen maken
# Hfundering
Expand All @@ -289,15 +332,19 @@ def Hgebouw(Bodem,Gebouw,Vloer):
# cZ = np.random.lognormal(0,var_c) *c;
# cX = np.random.lognormal(0,var_c_ratio)*cZ*c_ratio;
# Zg = 1j * gebouwdichtheidMC * gebouwhoogteMC * wandlengteMC * gevellengteMC * omega; # complexe impedantie van het gebouw
DictOut = Hfundering(Zb[i1],cX[i1],Zg[i1],wandlengteMC[i1]); # in x richting kloppen Z's eigenlijk niet
DictOut = Hfundering(Zb[i1],cX[i1],ZgX[i1],wandlengteMC[i1]); # in x richting kloppen Z's eigenlijk niet
Hfxx = DictOut['Htranslatie'];
DictOut = Hfundering(Zb[i1],cZ[i1],Zg[i1],wandlengteMC[i1]);
DictOut = Hfundering(Zb[i1],cZ[i1],ZgZ[i1],wandlengteMC[i1]);
Hfzz = DictOut['Htranslatie'];
Hfr = DictOut['Hrocking'];

# Hconstructie
Hcxx = 1.0; # doen we even niets mee, gaat over schuif en buig, beiden xx
Hczx = vloerhoogteMC[i1]*omega/cZ[i1]; # zwaaien, geometrische versterkin

maxvloerhoogte = gebouwC1HzMC[i1]/np.sqrt(octaafbanden)/4;
vloerhoogte = np.minimum(vloerhoogteMC[i1],maxvloerhoogte);

Hczx = vloerhoogte*omega/cZ[i1]; # zwaaien, geometrische versterkin
Hczz = 1.0; # hoogteverzwakking

# eerste resultaten: maaiveld naar x-richting bovenste verdieping
Expand Down Expand Up @@ -393,8 +440,16 @@ def Hgebouw(Bodem,Gebouw,Vloer):
Hfunvloer[2] = np.mean(Hv2mcArray,axis=0); # opslingering midspan
Hgebouw = np.max(Hfunvloer,axis=0); # volgens Level memo

# varianties uit de MC halen
# varianties uit de MC halen, eerst wat voorwerk
legeCov = np.zeros([6,6]);
# delen door nul is flauwekul,
Hxx[Hxx==0] = .001;
Hzx[Hzx==0] = .001;
Hzz1[Hzz1==0] = .001;
Hzz2[Hzz2==0] = .001;
Hfxx[Hfxx==0] = .001;
Hfzz[Hfzz==0] = .001;
# nu dus de covs maken
cov_Hxx = np.cov(HxxmcArray, rowvar=False) / (Hxx *np.transpose([Hxx ]));
cov_Hxz = legeCov;
cov_Hzx = np.cov(HzxmcArray, rowvar=False) / (Hzx *np.transpose([Hzx ]));
Expand Down Expand Up @@ -453,22 +508,28 @@ def Hfundering(Zbodem, cbodem, Zgebouw, Lgebouw):
aantallijnen = 10;
tauRockingff = np.ones(aantallijnen);
tauTransff = np.ones(aantallijnen);
# length of bending waves in foundation, assuming concrete or brick
# assuming 1 meter height; walls (up to underside window) can take part
# higher walls will "buckle"
# NB: double heigth shifts wave lengths one octave
lambdafundering = np.array([52,37,26,19,13,9.5]); # nb: vervangen door functie golflengte.py
for octaafnr in range(aantalbanden):
if cbodem[octaafnr]<1: # als van een band geen c is bepaald (dus =0)
buren = range(octaafnr-1,octaafnr+2,2); # info bij de buurbanden halen
octnr = range(aantalbanden); # bestaande buurbanden
buren = list(set(buren)&set(octnr)); # doorsnede
cbodem[octaafnr]=np.mean(cbodem[buren]);
ff = 2**(np.log2(octaafnr+1)+np.linspace(-.45,.45,aantallijnen));
LgebouwEff = min(Lgebouw,lambdafundering[octaafnr]);
labda = cbodem[octaafnr]/ff;
coeff = np.pi*Lgebouw/labda;
coeff = np.pi*LgebouwEff/labda;
for lijnnr in range(aantallijnen):
tauRockingff[lijnnr] = math.sin(coeff[lijnnr])/coeff[lijnnr];
tauTransff[lijnnr] = math.sqrt((1+math.cos(coeff[lijnnr]))/2);
tauRockingff[labda<Lgebouw] = 0;
tauTransff[labda<Lgebouw] = 0;
tauRockingOct[octaafnr] = np.mean(tauRockingff,axis=0);
tauTransOct[octaafnr] = np.mean(tauTransff,axis=0);
tauRockingff[lijnnr] = math.sin(coeff[lijnnr])/coeff[lijnnr];
tauTransff[lijnnr] = math.sqrt((1+math.cos(coeff[lijnnr]))/2);
tauRockingff[labda<Lgebouw] = 0;
tauTransff[labda<LgebouwEff/2] = 0;
tauRockingOct[octaafnr] = np.mean(tauRockingff,axis=0);
tauTransOct[octaafnr] = np.mean(tauTransff,axis=0);
# dan impedantiesprong nog in rekening brengen:
H = np.abs(Zbodem/(Zbodem+Zgebouw));
# effecten bijelkaar:
Expand Down